Thermal-fluid analysis of an RTG-powered ice penetrator for europa subsurface access: CFD modeling with phase-regime analysis

Abstract

Europa, Jupiter’s ice-covered moon, harbors a subsurface ocean representing a prime target for astrobiological investigation. This study presents computational fluid dynamics (CFD) analysis of a radioisotope-powered thermal ice penetrator (cryobot) designed for Europa subsurface access, with explicit treatment of the melting-sublimation phase regime transition. A dual-regime thermal model is developed: a sublimation-dominated surface regime (depth <0.5 m, pressure <611 Pa) and a pressure-assisted melting regime (depth >0.5 m) where the penetrator’s weight and meltwater column raise local interface pressure above the triple point. Two-dimensional axisymmetric CFD simulations using the enthalpy-porosity method predict steady-state descent rates of 0.19 mm/s in the sublimation regime and 0.92 mm/s in the melting regime under baseline conditions (4.4 kW thermal power, temperature-dependent ice conductivity k(T) = 651/T W/(mK)). A multi-angle parametric study (cone half-angles 30°, 45°, 60°, 75°, 90°) independently establishes geometry-dependent flux concentration factors ranging from 1.0 to 14.8, revealing that the 120° apex (60° half-angle) achieves η_geom,eff = 12.4 through combined geometric focusing and thermal conduction concentration. Richardson extrapolation on three systematically refined meshes confirms spatial convergence with Grid Convergence Index GCI_fine = 0.32%. Sensitivity analyses demonstrate descent rate varies linearly with power (±10% power yields ±9.1% rate change) and inversely with ambient ice temperature (4.7% improvement per 10 K). Critical refreezing analysis confirms safe operation (L_cr < 4 cm) throughout the 125–180 K range. Including the sublimation-limited surface transient and temperature-dependent properties, a 10 km penetration depth is achievable about 120 days (97–148 days). This work establishes validated thermal-fluid design parameters for kilometer-scale ice penetration missions to Europa and other ocean worlds.

1 Introduction

Europa, one of Jupiter’s four Galilean moons, has emerged as a premier target for astrobiological exploration due to compelling evidence of a global subsurface ocean beneath its ice shell (Pappalardo et al., 1999; Carr et al., 1998). Geophysical models constrained by Galileo spacecraft gravity data indicate a water-ice crust several kilometers to tens of kilometers thick overlying a briny ocean in contact with a silicate mantle (Wakita et al., 2024; Khurana et al., 1998), with independent ice-shell thickness estimates from elastically supported topography supporting this range (Nimmo et al., 2003). This configuration, comprising liquid water, chemical energy sources, and essential elements, meets the basic requirements for habitability as understood from terrestrial biology (McKay, 2004). The successful launch of the Europa Clipper from the space agency, i.e., the National Aeronautics and Space Administration (NASA) mission in 2024 designed for orbital reconnaissance, is an important milestone and forerunner of future in-situ surface and subsurface investigations (Korth et al., 2025; Hand et al., 2022).

1.1 Prior Europa lander concepts

Preliminary investigations into Europa surface access were done by Gershman et al. (Gershman et al., 2003), who put forward a baseline lander for shallow surface sampling. Subsequently, studies of radioisotope-powered landers were considered by the Jet Propulsion Laboratory of the National Aeronautics and Space Administration of the United States as part of the Jupiter Icy Moons Orbiter (JIMO) project (Abelson and Shirley, 2005). By 2011, joint plans of the joint missions by the National Aeronautics and Space Administration (NASA) and Russian Academy of Sciences under the Europa Jupiter System Mission (EJSM) were focused on modular lander designs with relay orbiters (Zelenyi et al., 2011). The most detailed mission architecture to date, presented by Hand et al. (Hand et al., 2022), proposed the use of a sky crane descent system for the delivery of a near-surface probe that could perform spectroscopic analysis. However, these concepts were limited to surface or near-surface sampling (less than 10 cm depth) and did not cover penetration through the multi-kilometer ice shell of Europa.

1.2 Thermal ice penetration technology

Mechanical drilling through kilometers of cryogenic ice poses formidable challenges such as torque reaction management, debris removal in vacuum and drill string logistics (Zacny et al., 2013). Thermal ice penetration (melt probes or “cryobots”) is an alternative technique where the thermal energy is used to melt or sublimate the ice and the probe can descend under gravity without cutting through (Philberth, 1962; Hansen et al., 1984). This methodology eliminates the problems associated with torque and allows constant operation without the use of depth limiting drill strings.

Several thermal penetrator concepts have been devised for terrestrial and planetary applications. The project VALKYRIE (Very-deep Autonomous Laser-powered Kilowatt-class Yo-yoing Robotic Ice Explorer) demonstrated the performance of laser powered ice penetration at rates approaching 1 mm/s in Antarctic field trials using a hemispherical tip (Stone et al., 2018). The THOR (Thermal High-voltage Ocean-penetrating Research) system investigated RTG power for Europa (We et al., 2011). The SLUSH (Subsurface Life Universal Search for Habitability) concept suggested the use of hybrid thermal-mechanical methods (Zacny et al., 2018). Schuller and Kowalski (Schüller and Kowalski, 2019) made basic analysis of critical refreezing length in the melt probes, which provide theoretical limitations for extraterrestrial ice penetration. Scully et al. (Scully et al., 2025) presented a comprehensive review of melting probe concepts for planetary ice exploration, which includes the IceMole technology demonstration. Dachwald et al. (2020) made an analysis of key technologies and system requirements for in-situ exploration of ice-covered ocean worlds. Tarau et al. (2022) reported the development status of an RTG powered thermal probe for icy planet exploration, addressing power and structural design considerations directly relevant to the present penetrator. Weiss et al. (2008) studied a thermal drill head design for subsurface planetary ice layer exploration, informing tip geometry considerations adopted in this work.

1.3 Limitations of prior work and present contributions

We do not claim that CFD, thermal-probe modeling, or cone-geometry studies are individually new; each has precedent in terrestrial and planetary ice-penetration research. The contribution of this work is their integration under Europa-specific conditions. The primary innovations are: (i) an explicit dual-regime model of the sublimation-to-melting transition at Europa surface pressure; (ii) a systematic parametric CFD study of cone half-angle that separates geometric concentration from thermal conduction focusing; and (iii) the use of temperature-dependent ice conductivity and specific heat as the baseline rather than constant properties. To the best of the authors knowledge, the combination of these three elements for a kilometer-scale Europa cryobot has not previously been reported.

Table 1 gives a quantitative comparison of this work to previous thermal penetrator concepts and points out the specific methodological improvements in the present work.

Taken together, these limitations define the research gap addressed in this study: no existing Europa-specific thermal-penetration model simultaneously resolves the sublimation-to-melting regime transition at Europa surface pressure, implements temperature-dependent ice properties, and optimizes tip geometry through systematic CFD with formal numerical verification. The present work is directed at closing this gap.

2 Phase regime analysis for Europa conditions

2.1 Triple point constraint and regime identification

A basic consideration of thermal penetration of Europa’s ice shell is the phase behavior of water ice at the prevailing thermodynamic conditions. The surface pressure of Europa is about 10⁻⁷ Pa (Spencer et al., 1999) which is far less than the triple point pressure of water (P_tp = 611.73 Pa, T_tp = 273.16 K). Under these conditions, ice heated above ∼150 K sublimates, passing directly from the solid phase to the vapor phase without going through an intermediate liquid phase (Bramson et al., 2015).

This places two different sets of operational regimes on the thermal penetrator:

Regime I—Sublimation-dominated (surface to depth z_tr): At depths where the local pressure is below 611 Pa, ice at the penetrator interface is sublimating. The relevant latent heat is , approximately 8.5 times larger than the latent heat of fusion . The descent rate in this regime is correspondingly slower.

Regime II—Pressure-assisted melting (): As the penetrator descends, two mechanisms raise the local interface pressure above the triple point: (a) the hydrostatic pressure of accumulated meltwater/refrozen ice above the penetrator, and (b) the contact pressure from the penetrator’s weight concentrated at the conical tip.

2.2 Transition depth estimation

The hydrostatic pressure at depth z in refrozen ice above the penetrator is:where

Additionally, the contact pressure at the conical tip provides supplementary pressurization. For a penetrator of mass m = 80 kg with effective contact area :

In addition to the hydrostatic pressure, the penetrator weight produces a localized mechanical contact stress at the conical tip. This contact stress is distinct from the uniform hydrostatic (thermodynamic) pressure of Equations 1, 2: it acts only over the small effective contact area at the tip and represents an average mechanical stress rather than a bulk fluid pressure. For a penetrator of mass m = 80 kg distributed over an effective contact area , the average contact stress is given by Equation 3.

This average contact stress (about 13.2 kPa) exceeds the triple-point pressure by a factor of about 21.5 within the immediate tip-contact region. Because it is a localized mechanical stress rather than a uniform hydrostatic pressure, it raises the local thermodynamic pressure only inside the small contact zone, where it helps sustain a thin liquid film; away from the tip the hydrostatic pressure of Equation 1 governs the phase state. To remain conservative, and because effective tip contact may be reduced during initial descent through porous regolith, we adopt approximately 0.5 m based on the hydrostatic criterion alone.

For the effective contact area derived in Section 4 ():

2.3 Descent rate in the sublimation regime

In Regime I, the energy balance at the interface requires sublimation rather than fusion, consistent with one-dimensional Stefan-problem formulations accounting for volume change and sensible heat (Huynh et al., 2024). The descent rate is:where replaces . Compared to the melting-regime rate:

For T₀ = 125 K:

As shown by Equation 7, the sublimation-regime rate of descent is around a fifth of the melting-regime rate, or about 20.5%. For the first 0.5 m of penetration, this is equivalent to a time penalty of:

The sublimation-phase time penalty in Equation 8 is evaluated using the same effective interface heat flux as in the melting regime. This is justified by the nature of the heat source: the RTG delivers a fixed thermal power to the tip through solid conduction in the high-conductivity cone, and this delivered power is independent of whether the interface ice melts or sublimates. The phase process changes only the energy required per unit mass removed (the latent-heat term in the denominator of Equation 5), not the energy supplied to the interface. The principal second-order effect neglected is the additional sensible heating of the escaping vapor; because the sublimation regime is confined to the first 0.5 m and adds under 1 hour to a mission lasting more than 120 days, this simplification has no practical effect on the duration estimate and is conservative for the limited surface transient.

2.4 Vapor transport considerations

In the sublimation regime, water vapor produced at the interface must be carried away to avoid back-pressure build-up. Under the conditions of near vacuum of Europa, the molecular mean free path is:where . The molecular mean free path is given by Equation 9, where (kinetic diameter of H₂O), and P = 10⁻⁷ Pa. At T = 200 K, λ ≈ 10⁶ m—vastly exceeding the probe dimensions. Thus, vapor transport takes place in the free-molecular flow regime (Knudsen number K_n >> 1) and sublimated molecules escape ballistically to the vacuum environment without developing significant back-pressure. The vapor does not interfere with the descent.

As the penetrator moves down into the melting regime (depth >0.5 m), the annular gap between the penetrator and borehole wall is filled with liquid water. The gap is sealed off from above with refrozen ice so that this results in a pressurized melt zone that keeps the liquid phase stable.

3 Thermal penetrator design

3.1 Design requirements and configuration

The extreme conditions of Europa’s surface (low gravity (1.315 m/s²), intense radiation (5.4 Sv/day), near-vacuum atmosphere (approximately 10⁻⁷ Pa), and cryogenic temperatures (86–132 K at the surface, with increasing depth)) have imposed stringent requirements on the design of the penetrator (Hand et al., 2022; Spencer et al., 1999). The following thermal penetrator is proposed: cylindrical thermal penetrator with diameter D = 0.6 m, total length L = 1.2 m, mass m ≈ 80 kg. The lower section ends in a conical tip of 120-degree apex angle (60-degree half-angle). The important design parameters are summarized in Table 2.

Figure 1 is a dimensioned schematic cross section of the thermal penetrator, illustrating the conical tip geometry, internal GPHS-RTG fuel modules, radial heat pipe array for thermal energy distribution and multi-layer insulation (MLI) with aerogel shielding.

FIGURE 1

3.2 Power source and thermal management

The power source is a General-Purpose Heat Source Radioisotope Thermoelectric Generator (GPHS-RTG) which produces 4.4 kW of thermal power and roughly 300 W of electrical power (Bennet et al., 2008). This power level is consistent with heritage RTG systems flown on Cassini, New Horizons and Mars Science Laboratory missions. The thermal to electrical conversion efficiency of about 7% is typical for thermoelectric systems operating at the relevant temperature differentials.

Thermal energy is conducted to the penetrator tip via an array of high conductivity copper heat pipes radially located around the bottom barrel. The non-contact surfaces are insulated with multi-layer insulation (MLI) which consists of 20 layers of Kapton with aluminized film separated by Dacron netting and supplemented with 10 mm thickness silica aerogel panels (Colgan et al., 2023; Durka et al., 2022). This configuration has an effective emissivity . This value corresponds to a 20-layer aluminized-Kapton MLI blanket supplemented with silica aerogel; measured effective emissivities for spacecraft MLI of this layer count typically lie in the range 0.01–0.03, increasing toward the upper end of the range when seams, attachment points, and edge losses are included (Gilmore, 2002; Colgan et al., 2023; Durka et al., 2022). We adopt the conservative upper bound for the radiative-loss estimate.

Radiative heat loss from the insulated surfaces is estimated by the Stefan Boltzmann law:As expressed in Equation 10, where ε ≈ 0.03, σ = 5.67 × 10⁻⁸ W/(m²·K⁴), A ≈ 2.4 m² (lateral surface area), T_s ≈ 500 K (insulated surface temperature), and T_env ≈ 125 K. This yields:Evaluating Equation 11 yields a radiative loss representing 5.7% of the 4.4 kW source power. Combined with axial conduction losses (∼8%) about 86% of thermal power is available for ice phase change.

3.3 Planetary protection considerations

A detailed planetary protection compliance assessment would be required as part of mission-level review and is beyond the scope of this thermal-fluid analysis.

4 CFD simulation methodology

4.1 Geometric model and computational domain

A two-dimensional axisymmetric model represents the thermal penetrator and surrounding ice domain. The penetrator is modeled as a conical section (baseline: 120° apex angle) connected to a cylindrical shaft with external diameter 0.6 m and total height 1.2 m. The computational domain extends 3 m radially (10× probe diameter) and 10 m axially (8.3× probe length) to ensure far-field boundary independence, verified by confirming that thermal perturbations at the outer boundaries remain below 0.1 K.

Figure 2 shows the computational domain with boundary condition annotations.

FIGURE 2

The domain represents the melting-regime conditions (depth >0.5 m) where the local pressure exceeds the triple point and liquid water exists at the interface. The sublimation regime is treated analytically (Section 2.3) and does not require separate CFD simulation because the sublimation front velocity follows directly from the energy balance (Equation 5) without convective transport in the vapor phase (free-molecular flow regime, Section 2.4).

4.2 Governing equations

The simulation solves the coupled energy, continuity, and momentum equations for conjugate heat transfer with phase change.

Energy equation:

Msomentum equation:

Continuity equation:

Phase change is captured using the enthalpy-porosity method (Voller and Prakash, 1987), consistent with coupled thermomechanical approaches used elsewhere to model ice removal (Song et al., 2020), wherein the liquid fraction β varies between 0 (solid) and 1 (liquid) within the mushy zone:where and are solidus and liquidus temperatures respectively. For pure water ice, ; a numerical mushy zone width of is adopted for computational stability, consistent with standard practice (Voller and Prakash, 1987).

Applicability of the enthalpy-porosity method: This model is applied exclusively in Regime II (depth >0.5 m) where local pressure exceeds 611 Pa and liquid water exists at the interface. The pressure in the sealed melt zone is maintained by the combination of hydrostatic head and contact pressure (Equations 1–4). The enthalpy-porosity method is well-validated for this configuration (Voller and Prakash, 1987; Brent et al., 1988).

The source term incorporates latent heat:

The momentum source term implements the Carman-Kozeny equation to suppress velocity in solid regions:where C = 10⁵ is the mushy zone constant and ε = 10⁻³ prevents division by zero. The sensitivity of results to the mushy zone constant was verified: varying C from 10⁴ to 10⁶ produced less than 2% change in predicted descent rate, consistent with literature findings (Voller and Prakash, 1987).

4.3 Temperature-dependent material properties

Unlike prior studies using constant ice properties, the present work employs temperature-dependent thermal conductivity as the baseline:

The conductivity correlation of Equation 18 follows the established relationship for crystalline ice (Slack, 1980), (Ohno et al., 2006). Alternative published correlations for polycrystalline ice differ from this form by less than about 10 percent over 125–273 K. Because the descent rate scales sub-linearly with conductivity (normalized sensitivity minus 0.6, Table 7), a 10 percent change in the correlation alters the predicted rate by about 6 percent, within the plus or minus 15 percent conductivity uncertainty already in the budget. The specific-heat correlation is the standard linear fit for ice (Ohno et al., 2006) and affects the rate only through the sensible-heating term (normalized sensitivity minus 0.35); alternative Cp correlations change the predicted rate by less than 4 percent. The sublimation-to-melting transition depth depends on ice density and gravity rather than on these correlations and is therefore insensitive to the choice of property model. Table 3 lists the complete material property specifications.

Specific heat capacity is also temperature-dependent (Ohno et al., 2006):yielding at 125 K and at 273 K. Evaluating Equation 19 yields the specific heat at 125 K and at 273 K. Ice density is treated as constant (ρ = 920 kg/m³) as the variation with temperature is less than 1% over the relevant range.

4.4 Boundary conditions

The gravity vector is set to g = 1.315 m/s² (Europa) directed axially downward to capture buoyancy-driven convection in the melt layer. The Rayleigh number based on the melt film thickness δ ≈ 2 mm is:where is the thermal expansion coefficient and m²/s is the thermal diffusivity. Since Ra <1708 (critical value for onset of convection), the melt film is in the conduction-dominated regime and natural convection is negligible. This justifies the enthalpy-porosity approach without turbulence modeling.

Note: Equations 12–14 express conservation of energy, momentum, and mass for the coupled solid-liquid system. In the energy Equation 12, the first left-hand term is the local rate of change of sensible enthalpy and the second is enthalpy advection by the melt flow; on the right, the first term is conduction with temperature-dependent conductivity and is the latent-heat source associated with phase change (Equation 16). In the momentum Equation 13, the terms represent, in order, transient inertia, advection, the pressure gradient, viscous diffusion, the gravitational body force, and the porous-media damping source (Equation 17) that drives velocity to zero in fully solid cells. The continuity Equation 14 enforces mass conservation. The set is closed by the enthalpy-porosity relations (Equations 15–17), in which the liquid fraction beta maps temperature to the local phase state and progressively suppresses momentum as cells solidify. Because the melt film is thin and conduction-dominated (Rayleigh number approximately 86, well below the convective threshold of 1708; Equation 20), advective terms are small and the interface energy balance reduces essentially to conduction balanced by latent heat.

4.5 Mesh generation and grid convergence study

The computational domain is discretized using structured quadrilateral elements in the near-probe region and unstructured triangular elements in the far field. Mesh refinement is concentrated at the ice-penetrator interface where the largest thermal gradients occur, with minimum element size of 0.125 mm (fine mesh) resolving the phase-transition boundary.

Grid independence was verified using three systematically refined meshes with refinement ratio r = 2.0, as summarized in Table 4:

The observed order of convergence is calculated following Roache (Roache, 1994).

The Richardson-extrapolated value is:

The Grid Convergence Index for the fine mesh is:where is the safety factor for three-grid studies. The fine mesh achieves less than 0.4% numerical uncertainty, confirming spatial convergence. The Richardson-extrapolated value of 0.915 mm/s is adopted as the grid-independent reference; the fine-mesh solution (0.917 mm/s) deviates by only 0.2%.

The fine mesh (270,000 cells) is used for all subsequent analyses. For reporting purposes, the baseline descent rate is quoted as 0.92 mm/s (rounded from the extrapolated 0.915 mm/s).

In summary, the grid-convergence study yields an observed order of convergence (consistent with the second-order discretization), a Richardson-extrapolated grid-independent descent rate of 0.915 mm/s, and a fine-grid Grid Convergence Index , corresponding to an estimated numerical uncertainty below 0.4 percent. The fine-mesh solution (0.917 mm/s) lies 0.2% above the extrapolated value. The grid-independent reference rate of 0.915 mm/s is used consistently throughout the manuscript and is reported as 0.92 mm/s where rounded to two significant figures.

4.6 Solver configuration

Simulations were performed using ANSYS Fluent 2024 R² with a pressure-based coupled algorithm. Temporal discretization employed a second-order implicit scheme with adaptive time stepping: initial time step during the transient phase-transition period, gradually increasing to 0.1 s as steady-state conditions developed. Convergence criteria were set to scaled residuals below 10⁻⁶ for energy and 10⁻⁴ for continuity and momentum equations, combined with monitoring of interface heat flux variation (<0.1% between consecutive iterations). Each simulation was run for a physical time of 100,000 s (∼27.8 h) to ensure fully developed quasi-steady-state conditions. Total wall-clock computation time was approximately 48 h per case on a workstation with dual Intel Xeon Gold 6248 R processors.

5 Results and analysis

5.1 Thermal field evolution and descent rate

Figure 3 presents the steady-state temperature contours from the baseline simulation (4.4 kW, 120° apex, T₀ = 125 K, k(T) = 651/T). Three regions are labeled. Region A, directly beneath the conical tip, is the liquid melt zone where T > 273.15 K and the liquid fraction beta = 1; it forms a compact cap about 10 cm across that conforms to the cone. Region B is the ice-water interface, the 273.15 K isotherm, which appears as a thin curved front (melt-film thickness 1.5–2.5 mm) wrapping the tip. Region C is the far-field conduction zone, where the isotherms form nested, elliptically stacked surfaces characteristic of axisymmetric conduction-dominated heat transfer. The tight clustering of isotherms at the apex relative to the cone flank is the visual signature of the flux concentration quantified in Section 5.2.

FIGURE 3

Figure 4 presents the temporal evolution of the descent rate. Three stages are visible: an initial start-up transient (0 to about 5,000 s) during which the melt cap forms and the rate rises steeply from zero; an approach stage (about 5,000–25000 s) in which the rate asymptotes as the thermal field around the tip matures; and a quasi-steady stage (beyond about 25000 s, roughly 7 h of simulated time) in which the rate settles at 0.92 mm/s with variation below 0.1 percent between iterations. The steady value reached in Figure 4 is the rate reported throughout the remainder of the paper.

FIGURE 4

5.2 Multi-angle parametric study and geometry factor derivation

To independently establish the relationship between cone geometry and effective heat flux concentration, CFD simulations were performed for five cone half-angles: 30°, 45°, 60°, 75°, and 90° (flat), holding all other parameters constant (4.4 kW, 125 K, k(T)). This parametric study is the primary means of characterizing the geometry effect, replacing the circular single-point calibration of prior work.

For each geometry, the conical surface area varies while the projected circular area remains constant. The geometric area ratio provides a first-order flux concentration estimate:where θ is the cone half-angle. However, the effective flux concentration η_geom,eff is defined from the CFD results as:where is the one-dimensional Stefan descent rate for uniform flux over the full cross-sectional area:

Note: Two complementary measures of flux concentration are used. Equation 24 is a first order, purely geometric estimate equal to the ratio of the fixed projected area A_proj to the cone surface area Acone, which for an ideal cone reduces to sin (theta) and bounds the concentration available from area reduction alone. Equation 25 defines the effective concentration factor geom,eff extracted from the CFD results, namely, the ratio of the CFD descent rate for a given cone to the one-dimensional flat-tip Stefan rate of Equation 26. Equations 24, 25 are therefore related as a lower-bound geometric estimate Equation 24 and the full CFD-derived value Equation 25; the gap between them isolates the contribution of thermal conduction focusing within the metallic tip. In Equation 26, Q_net is the net thermal power delivered to the ice interface after subtracting radiative and axial-conduction losses from the RTG output, Q_net = 0.86 times 4.4 kW = 3.784 kW, the same value applied as the tip boundary condition in Section 4.4.

Table 5 presents the results of the multi-angle parametric study. The effective concentration factor combines three distinct effects, which the CFD results allow us to separate. First, geometric concentration: the conical surface projects the delivered power onto a smaller cross-sectional area, giving a first-order factor 1/sin (theta) that rises only modestly as the cone sharpens (1.0 at 90° to 2.0 at 30° half-angle). Second, contact-area variation: the effective high-flux contact zone at the apex, where the flux exceeds 50 percent of maximum, shrinks as the cone sharpens, concentrating power over a smaller footprint. Third, thermal conduction focusing: heat conducted through the high-conductivity metallic cone converges toward the apex; this is the dominant effect and accounts for the order-of-magnitude enhancement over the geometric factor alone.

These effects explain why the 30-degree cone does not outperform the 60-degree cone despite its higher geometric concentration factor. As the cone sharpens below 45°, three penalties grow and offset the geometric gain: the longer cone flank increases lateral heat loss; the extended conduction path from the heat-pipe array to the apex raises the internal temperature drop and reduces the flux delivered to the contact zone; and the smaller apex footprint limits the rate at which melted volume can be cleared. The descent rate increases by only 18 percent from the 45° to the 30-degree half-angle even though the cone surface area increases by 41 percent, a clear sign of saturation. The 60-degree half-angle captures most of the available rate benefit while avoiding the buckling risk, reduced internal volume, and manufacturing difficulty of sharper cones, and is therefore selected as the practical optimum.

Figure 5 plots η_geom,eff versus cone half-angle, revealing several important features.

FIGURE 5

Physical interpretation of η_geom,eff = 12.4:

Thermal-conduction focusing factor is the part of the effective flux concentration that arises from heat being conducted through the high-conductivity metallic cone and concentrated at the apex contact region, as distinct from the purely geometric area ratio 1 = 1/sin (theta). For the 120-degree apex, 1 = 1/sin (60°) = 1.155 and geom,eff = 12.4, so 2 = 12.4/1.155 = 10.7. A value of 2 greater than one indicates that conduction focusing, not area reduction, is the dominant concentration mechanism.

Heat conducting via the high conductivity copper alloy cone (k_tip = 320 W/(mK)) is focused on the apex contact region. This effect is quantified by analyzing the distribution of the heat flux on the surface of the cone from the results of the CFD calculations.

Figure 6 plots the local heat flux q along the cone surface as a function of normalized arc length s/s_max, from the apex (s = 0) to the base (s = s_max). The curve is annotated to show (i) the peak apex flux, (ii) the 50 percent of maximum threshold that defines the effective contact area A_eff, and (iii) the base flux, which is about one-11th of the apex value. The shaded region beneath the curve up to the 50 percent threshold marks the effective contact zone (A_eff approximately 0.008 square meters, about 10 cm in diameter). The steep decay from apex to base is direct evidence of the conduction focusing quantified by 2.

FIGURE 6

The heat flux at the apex is approximately 11× higher than at the cone base, confirming that thermal conduction within the metallic cone strongly concentrates energy at the contact point. The effective contact area (where heat flux exceeds 50% of maximum) is , corresponding to a circular contact diameter of ∼10 cm. This is consistent with the localized high-temperature zone observed in the temperature contours.

Comparison with analytical conical heat source solutions: Carslaw and Jaeger (1959) provide solutions for steady heat conduction from conical surfaces showing flux concentration proportional to (1/r) near the apex, where r is the distance from the tip. For a cone of thermal conductivity k_tip embedded in a medium of conductivity k_ice, the concentration effect scales approximately as:where f(θ) is a geometry-dependent function. For and θ = 60°, analytical estimates yield depending on contact assumptions (Carslaw and Jaeger, 1959), bracketing the CFD-derived value of 10.7.

Where f(θ) is a dimensionless geometry function that captures how the apex half-angle theta modifies the conduction concentration; it increases as the cone sharpens (smaller theta). It is not evaluated in closed form here. Equation 27 is used only as an order-of-magnitude scaling check: with the conductivity ratio and θ = 60°, analytical estimates yield depending on contact assumptions (Carslaw and Jaeger, 1959), bracketing the CFD-derived value of 10.7. A precise functional form for f (theta) would require the full series solution in (Carslaw and Jaeger, 1959) and is beyond the present scope. Here denotes the radial distance measured from the cone apex (the tip), following the conical-source solution of Carslaw and Jaeger (1959), in which the conduction flux near a conical tip scales as It is a spatial coordinate associated with the conical geometry, not the cone radius or a material property; the dependence is precisely what produces the apex flux concentration

5.3 Validation against published experimental data

The CFD model is validated against published experimental data by simulating the specific conditions of each experiment (geometry, power, temperature) rather than applying a single universal correction factor. Three distinct levels of credibility should be distinguished. Numerical verification, established in Section 4.5 through Richardson extrapolation and the Grid Convergence Index, confirms that the discrete solution converges to the continuum solution of the model. Benchmarking, presented here, shows agreement within 4–7 percent against terrestrial and Mars-chamber experiments and builds confidence in the physical model. True validation under Europa surface conditions (about 125 K and 10 to the minus 7 Pa) is not currently possible and remains an objective for future cryogenic-vacuum testing. The 4 to 7 percent agreement therefore constitutes benchmarking under terrestrial and Mars-relevant conditions, not validation under Europa conditions. Table 6 presents the direct comparison.

For each case in Table 6, validation was performed by reproducing the reported experiment directly in the CFD model rather than by applying a single universal scaling factor. For VALKYRIE (Stone et al., 2018), a 0.25 m hemispherical tip at 5.0 kW and 263 K was modeled and the steady descent rate compared with the Antarctic field measurement. For Kömle et al. (2018), a 0.06 m parabolic tip at 0.5 kW and 253 K was modeled to match their Mars-pressure laboratory ice-penetration tests. For Schüller and Kowalski (2019), a 0.12 m 90-degree conical tip at 1.0 kW and 263 K was modeled against their reported rate. In every case the domain size, boundary conditions, and temperature-dependent properties were set to the reported experimental values, so the deviation column reflects a like-for-like prediction. Ice properties were set to terrestrial values at the reported temperatures (k(T) = 651/T, ρ = 917 kg/m³).

All predictions fall within 4%–7% of measured values, which is within the reported experimental uncertainties (±8%–14%). This provides confidence in the numerical methodology despite the inability to directly validate under Europa conditions.

Validation limitations: Direct experimental confirmation under Europa surface conditions (125 K, 10⁻⁷ Pa) is not currently possible. The validation cases were conducted at terrestrial temperatures (253–263 K) and at pressures ranging from atmospheric (terrestrial field tests) to about 600 Pa (the Mars-chamber conditions of Kömle et al. (2018)). These pressures pertain to the validation experiments and are unrelated to the Europa surface pressure of about 10⁻⁷ Pa.

These uncertainties are partially captured in the sensitivity analysis (Section 5.4) and propagated into the mission duration estimate.

5.4 Sensitivity analysis

5.4.1 Influence of RTG thermal power

Figure 7 plots the descent rate against RTG thermal power for both the CFD points and the one-dimensional Stefan model carrying the geometry factor (). The two data sets coincide closely, and the linear least-squares fit (Equation 28, R squared = 0.998) confirms the expected proportionality between descent rate and thermal power. The figure links the system input (thermal power) to the predicted response (descent rate): each 10 percent change in delivered power maps to a 9.1 percent change in descent rate, and the nominal 4.4 kW operating point yields 0.92 mm/s, rising to 1.04 mm/s at 5.0 kW. The close agreement between the CFD points and the analytical line also confirms that the geometry factor extracted in Section 5.2 transfers consistently across the full power range

FIGURE 7

A linear fit to the CFD data points yields:confirming the analytical proportionality . A ±10% power variation yields ±9.1% descent rate change. At the nominal 4.4 kW, the predicted rate is 0.92 mm/s; increasing power to 5.0 kW increases the rate to 1.04 mm/s.

5.4.2 Influence of ambient ice temperature

Figure 8 shows the descent-rate sensitivity to ambient ice temperature over 100–200 K. The dependence is mildly nonlinear at the cold end (100–130 K) and becomes nearly linear above about 140 K. The nonlinearity is not a true polynomial law but arises from two temperature-dependent properties acting in opposition: the ice thermal conductivity k(T) = 651/T, which falls with temperature and reduces far-field heat loss (aiding descent), and the specific heat Cp(T) = 7.49 T + 90, which rises with temperature and increases the sensible-heating penalty. Because the rate depends on the combination of these properties, the response is a smooth rational function of temperature rather than a polynomial; at cryogenic temperatures the rapid 1/T rise in conductivity dominates and produces the visible curvature, while at higher temperatures the two effects nearly offset to give the near-linear trend.

FIGURE 8

The rate increases with ambient temperature due to reduced sensible heating requirements (smaller ) and reduced ice thermal conductivity (less heat lost to far-field conduction). About 4.7% gain in descent rate is obtained for each 10 K rise in ambient temperature. This favorable trend shows an increasing depth of good penetration as the probe enters warmer and warmer ice layers.

The nonlinearity at low temperatures (100–130 K) represents the combination of increased value for k(T) and higher C_p penalty at cryogenic conditions. At 100 K, k = 6.51 W/(mK) and ; at 200 K, and . These two effects act in opposition: the decrease in thermal conductivity with temperature reduces far-field heat loss and increases the descent rate, whereas the increase in specific heat raises the sensible-heating requirement and decreases it. The two contributions partially cancel, producing the near-linear dependence observed at higher temperatures

5.4.3 Comparison: temperature-dependent vs. constant properties

In order to quantify the significance of the implementation of the k(T), prediction of descent rates with temperature-dependent and constant (k = 2.2 W/(m*K)) properties were compared in the ambient temperature range in Figure 9.

FIGURE 9

At (Europa surface) the constant property model over-predicts descent rate by 22% due to underestimated thermal conductivity of ice (), which results in underestimated heat losses. At the discrepancy is lowered to 8%. This confirms that temperature-dependent properties are crucial to the accurate modeling of Europa ice penetration, especially in the cold near-surface region.

The constant value k = 2.2 W/(m K) is used as the comparison baseline because it is the conductivity of polycrystalline water ice near the melting point (about 273 K), the value most commonly adopted in prior constant-property thermal-penetrator studies. It is consistent with the standard tabulated conductivity of ice I_h at 273 K (Slack, 1980) and with the value returned by the present temperature-dependent relation at the melting point (k (273 K) = 651/273 = 2.38 W/(m K)). Using this near-melting-point constant is the most favorable case for the constant-property assumption, which is why it still over-predicts the descent rate at the cold surface, where the true conductivity is more than twice as large.

5.5 Uncertainty propagation

A structured uncertainty analysis was performed by propagating parameter uncertainties through the validated CFD-based descent rate correlation (Equation 28) combined with the analytical sensitivity relationships. Table 7 shows the uncertainty budget for baseline descent rate.

The combined uncertainty (root-sum-square) yields a baseline descent rate of 0.92 ± 0.18 mm/s (±19%, 95% confidence assuming 2σ). The geometry factor is the dominant contributor for two compounding reasons. First, it carries the largest normalized sensitivity in Table 7 (plus 1.0): the descent rate is directly proportional to geom,eff, so a given fractional error transfers one-for-one into the rate, whereas the other parameters have sub-unity normalized sensitivities (for example, −0.6 for conductivity and −0.35 for specific heat). Second, it carries one of the largest input uncertainties (±15 percent), because geom,eff is obtained from CFD calibrated against terrestrial experiments and must be extrapolated to Europa cryogenic, low-pressure conditions, where contact area and conduction focusing are less certain. The product of unity sensitivity and 15 percent input uncertainty (±0.138 mm/s) exceeds every other line in the budget and sets the overall ±0.175 mm/s root-sum-square value. Cryogenic validation of the geometry factor would therefore give the largest single improvement in predictive confidence.

5.6 Critical refreezing length analysis

During descent, meltwater in the annular gap between the penetrator and borehole wall refreezes due to heat conduction into the surrounding cold ice. If the refreezing length exceeds the penetrator length, borehole closure could trap the probe.

Following Schüller and Kowalski (2019), the critical refreezing length is estimated from Equation 29:where is the lateral heat flux from the penetrator’s cylindrical surface. For the insulated design (ε = 0.03), consists primarily of radiative losses, as given by Equation 30:

At the cone tip, the relevant heat flux is much higher (concentrated q'' from RTG). The refreezing length behind the tip is given by Equation 31:

For T₀ = 125 K and the effective tip flux = 3784/0.008 = 4.73 × 10⁵ W/m²:

The analytical estimate (Equation 32) and the CFD prediction are complementary. The analytical formula gives a lower-bound steady-state refreezing length of 0.16 cm at the tip. The CFD simulations capture the transient, multi-dimensional refreezing behind the moving probe and predict a larger zone of 2–4 cm, because they include latent-heat redistribution and the finite time before meltwater is carried below the freezing front. The CFD value is the more conservative and physically complete estimate and is used to assess closure risk. A transient refreezing zone of 2–4 cm does not threaten borehole stability, because it is far shorter than the 120 cm penetrator and forms behind the descending probe rather than ahead of it. The ratio of penetrator length to worst-case refreezing length (120/4 = 30) is the relevant closure-margin metric and remains above 30 throughout the 10 km descent (Figure 10), confirming that borehole closure does not constrain mission feasibility.

FIGURE 10

Spencer et al. (1999); the length stays below 4 cm throughout, well under the 120 cm penetrator length5.7 Mission Duration Estimate.

The total penetration time includes the sublimation-regime surface transient and the melting-regime bulk penetration, is expressed in Equation 33:where is the depth-averaged melting-regime descent rate accounting for the temperature profile. Here is the target penetration depth, that is the full ice-shell thickness to be traversed, taken as 10 km (10,000 m) in this paper.

Using the thermal gradient from Spencer et al. (1999) (approximately 10 K/km for the conductive ice shell):

Numerical integration using the sensitivity relationship from Section 5.4.2 yields ≈ 0.97 mm/s (accounting for the favorable temperature increase with depth). The depth-averaged melting-regime rate in Equation 34 is evaluated as follows. The conductive ice shell has an approximately linear temperature rise with depth of about 10 K/km (Spencer et al., 1999), so . At each depth the local descent rate at temperature is obtained from the ambient-temperature sensitivity established in Section 5.4.2 (Figure 8), which expresses the rate as a function of local ice temperature. Substituting into that relation gives the rate as a function of depth, which is integrated numerically (trapezoidal rule, 100 m steps) from 0 to and divided by to obtain the depth-averaged rate. Because the ice warms with depth and a warmer ambient raises the rate, the depth-averaged value (0.97 mm/s) is slightly higher than the 125 K surface value (0.92 mm/s). This averaged rate is then used in Equation 36.

Equations 35–37 show that the sublimation-regime surface transient ( = 0.74 h) is negligible relative to the melting-regime bulk penetration ( = 119.3 days), contributing less than 0.001 percent of the total mission time. This result is significant: although sublimation is energetically far more demanding than melting (latent heat larger by a factor of about 8.5), it is confined to only the first 0.5 m of a 10 km descent, so the slow sublimation regime imposes essentially no penalty on the overall mission duration. The phase-regime transition therefore matters for the physical correctness of the interface model but not for the mission timeline, which is governed almost entirely by the melting regime.

Incorporating the ±19% uncertainty on descent rate:or approximately 97–148 days (95% confidence).

Figure 11 shows the mission timeline with uncertainty bands.

FIGURE 11

Pu-238 decay correction: Over a 6-year Earth-Jupiter transit, Pu-238 power decreases by approximately 5.2% (half-life 87.7 years). This is accommodated by the +10% power margin assumed in the uncertainty analysis. Over the 120-day surface mission, decay is negligible (<0.3%).

This mission-duration figure should be interpreted as a model-based approximation rather than a fixed prediction. It is conditional on the principal modeling assumptions: axisymmetric geometry, quasi-steady descent, homogeneous ice, and the geometry factor extrapolated from terrestrial validation. Propagating the ±19 percent descent-rate uncertainty gives a range of approximately 97–148 days at 95% confidence (Equation 38). The assumptions that most strongly influence the estimate are the geometry factor and the ice thermal conductivity, followed by the assumed thermal gradient with depth. Ice heterogeneity (Section 6.2) could shift the estimate by a further 5%–15%, so about 120 days should be regarded as a central estimate within a plausible range of roughly 90–160 days once model-form uncertainty is included.

6 Discussion

6.1 Comparison with prior Europa penetrator analyses

The predicted baseline descent rate of 0.92 mm/s is consistent with the range reported by prior Europa penetrator studies (0.3–1.0 mm/s for 1–5 kW systems, Table 1) and, to the best of the authors knowledge, represents one of the first values derived from CFD analysis with temperature-dependent properties and formal grid convergence documentation.

The multi-angle parametric study provides, to the best of the authors knowledge, one of the first systematic CFD-based quantifications of cone-angle effects on thermal penetrator performance.

The multi-angle parametric study (Section 5.2) provides, to the authors’ knowledge, the first systematic CFD-based quantification of cone angle effects on thermal penetrator performance. The result that effective flux concentration is much higher than geometric area ratios - by an order of magnitude for the 120deg apex - has important implications for penetrator design optimization. It demonstrates that tip material thermal conductivity is as important as tip geometry in determining descent rate.

6.2 Effects of ice heterogeneity

A comprehensive treatment of heterogeneous ice is beyond the scope of this work but represents an important direction for future 3D simulations.

6.3 Limitations

6.4 Design implications

Three engineering conclusions follow from the analysis. First, the parameters with the greatest influence on descent rate are the delivered thermal power, to which the rate is directly proportional (Equation 28), and the tip geometry factor; ambient ice temperature is a secondary influence (about 4.7 percent per 10 K). Second, the assumptions contributing most to predictive uncertainty are the geometry factor (±15 percent) and the ice thermal conductivity (±15 percent), which dominate the ±19 percent rate uncertainty, with the axisymmetric assumption adding a further estimated 10 to 15 percent for heterogeneous ice. Third, the principal engineering compromise is between tip sharpness and robustness: sharper cones marginally increase the rate but reduce structural margin and internal volume, so the 60-degree half-angle is preferred. Future designs should prioritize maximizing delivered thermal power and validating the geometry factor under cryogenic conditions, while retaining sufficient structural and planetary-protection thermal margins.

7 Conclusion

Future work should address: (1) cryogenic-vacuum chamber experiments at T < 150 K to validate the sublimation-to-melting transition model; (2) 3D CFD with heterogeneous ice (salt inclusions, clathrates, cracks) to quantify off-axis effects; (3) coupled melt-film lubrication analysis for high-fidelity descent rate prediction; and (4) ice-ocean interface transition dynamics including buoyancy control.

References

• 1

  AbelsonR. D.ShirleyJ. H. (2005). Exploring Europa with a surface lander powered by a small radioisotope power system (RPS). AIP Conf. Proc.746 (1), 303–317. 10.1063/1.1867147

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  BennettG. L.LombardoJ. J.HemlerR. J.SilvermanG.WhitmoreC. W.AmosW. R.et al (2008). The general-purpose heat source radioisotope thermoelectric generator: a truly general-purpose space RTG. Proc. AIP Conf.969 (1), 663–671. 10.1063/1.2845028

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BramsonA. M.ByrneS.PutzigN. E.SuttonS.PlautJ. J.BrothersT. C.et al (2015). Widespread excess ice in Arcadia planitia, Mars. Geophys. Res. Lett.42 (16), 6566–6574. 10.1002/2015gl064844

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BrentA. D.VollerV. R.ReidK. J. (1988). Enthalpy-porosity technique for modeling convection-diffusion phase change: application to the melting of a pure metal. Numer. Heat. Transf.13 (3), 297–318. 10.1080/10407788808913615

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  CarrM. H.BeltonM. J.ChapmanC. R.DaviesM. E.GeisslerP.GreenbergR.et al (1998). Evidence for a subsurface ocean on europa. Nature391 (6665), 363–365. 10.1038/34857

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  CarslawH. S.JaegerJ. C. (1959). Conduction of Heat in Solids. 2nd ed.Oxford, UK: Clarendon Press.

  • Google Scholar
  • View reference in article

• 7

  ColganW.ShieldsC.TalalayP.FanX.LinesA. P.ElliottJ.et al (2023). Design and performance of the hotrod melt-tip ice-drilling system. Geosci. Instrum. Methods Data Syst.12 (1), 121–140. 10.5194/gi-12-121-2023

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  DachwaldB.UlamecS.PostbergF.SohlF.de VeraJ. P.WaldmannC.et al (2020). Key technologies and instrumentation for subsurface exploration of ocean worlds. Space Sci. Rev.216 (5), 83. 10.1007/s11214-020-00707-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  DurkaM. J.SmithM. W.UllmanM. J.CasslerB.OtisR.CwikT. A.et al (2022). Modeling a class of thermal ice probes for accessing the solar system's ocean worlds. Acta Astronaut.193, 483–495. 10.1016/j.actaastro.2021.12.018

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  GershmanR.NilsenE.ObertoR. (2003). Europa lander. Acta Astronaut.52 (2–6), 253–258. 10.1016/s0094-5765(02)00164-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  GilmoreD. G. (2002). Spacecraft Thermal Control Handbook, Volume I: Fundamental Technologies. 2nd ed.El Segundo, CA, USA: The Aerospace Press.

  • Google Scholar
  • View reference in article

• 12

  HandK. P.PhillipsC. B.MurrayA.GarvinJ. B.MaizeE. H.GibbsR. G.et al (2022). Science goals and mission architecture of the europa lander mission concept. Planet. Sci. J.3 (1), 22. 10.3847/psj/ac4493

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  HansenB. L.KerstenL. (1984). “An in-situ sampling thermal probe,” in Ice-Core Drilling. Editor SplettstoesserJ. F. (Lincoln, NE, USA: Univ. Nebraska Press), 119–129.

  • Google Scholar
  • View reference in article

• 14

  HuynhM.-T. T.NguyenH. T.GjerløwE. (2024). An alternative two-phase formulation and analytical solution for a stefan problem in one-dimension accounting for volume change and sensible heat. Int. J. Thermofluids22, 100692. 10.1016/j.ijft.2024.100692

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  KhuranaK. K.KivelsonM. G.StevensonD. J.SchubertG.RussellC. T.WalkerR. J.et al (1998). Induced magnetic fields as evidence for subsurface oceans in europa and callisto. Nature395 (6704), 777–780. 10.1038/27394

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  KömleN. I.TiefenbacherP.WeissP.BendiukovaA. (2018). Melting probes revisited—Ice penetration experiments under Mars surface pressure conditions. Icarus308, 117–127. 10.1016/j.icarus.2017.08.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  KorthH.BurattiB. J.SenskeD. A. (2025). Europa clipper: a mission to explore ocean world habitability (editorial). Space Sci. Rev.221 (2), 1–3. 10.1007/s11214-025-01150-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  McKayC. P. (2004). What is Life—And how do we search for it in other worlds?PLoS Biol.2 (9), e302. 10.1371/journal.pbio.0020302

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  National Research Council (2012). Assessment of Planetary Protection Requirements for Spacecraft Missions to Icy Solar System Bodies. Washington, DC: The National Academies Press.

  • Google Scholar
  • View reference in article

• 20

  NimmoF.GieseB.PappalardoR. T. (2003). Estimates of Europa's ice shell thickness from elastically-supported topography. Geophys. Res. Lett.30 (5), 1233. 10.1029/2002gl016660

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  OhnoH.IgarashiM.HondohT. (2006). Characteristics of salt inclusions in polar ice from dome fuji, east Antarctica. Geophys. Res. Lett.33 (8), L08501. 10.1029/2006gl025774

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  PappalardoR. T.BeltonM. J. S.BrenemanH. H.CarrM. H.ChapmanC. R.CollinsG. C.et al (1999). Does europa have a subsurface ocean? Evaluation of the geological evidence. J. Geophys. Res. Planets104 (E10), 24015–24055. 10.1029/1998je000628

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  PhilberthK. (1962). Une méthode pour mesurer les températures à l'intérieur d'un inlandsis. C. R. Hebd. Seances Acad. Sci.254, 3881–3883.

  • Google Scholar
  • View reference in article

• 24

  RoacheP. J. (1994). Perspective: a method for uniform reporting of grid refinement studies. J. Fluids Eng.116 (3), 405–413. 10.1115/1.2910291

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  SafiE.ThompsonS. P.EvansA.DayS. J.MurrayC. A.ParkerJ. E.et al (2017). Properties of CO₂ clathrate hydrates formed in the presence of MgSO₄ solutions with implications for icy moons. Astron. Astrophys.600, A73. 10.1051/0004-6361/201629791

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  SchüllerK.KowalskiJ. (2019). Melting probe technology for subsurface exploration of extraterrestrial ice—Critical refreezing length and the role of gravity. Icarus317, 1–9. 10.1016/j.icarus.2018.05.022

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  ScullyJ. E. C.BelgacemI.ParekhR.GrimaC.PhillipsC. B.CraftK. L.et al (2025). Potential landing sites: a comprehensive reconnaissance assessment of the europa clipper trajectory, planet. Sci. J.6 (6), 142. 10.3847/psj/add469

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  SlackG. A. (1980). Thermal conductivity of ice. Phys. Rev. B22 (6), 3065–3071. 10.1103/physrevb.22.3065

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  SongY.LiuR.LiS.KangZ.ZhangF. (2020). Peridynamic modeling and simulation of coupled thermomechanical removal of ice from frozen structures. Meccanica55 (4), 961–976. 10.1007/s11012-019-01106-z

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  SpencerJ. R.TamppariL. K.MartinT. Z.TravisL. D. (1999). Temperatures on europa from galileo photopolarimeter-radiometer: nighttime thermal anomalies. Science284 (5419), 1514–1516. 10.1126/science.284.5419.1514

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  StoneW.HoganB.SiegelV.HarmanJ.FlesherC.ClarkE.et al (2018). “Project VALKYRIE: laser-Powered cryobots and other methods for penetrating deep ice on ocean worlds,” in Outer Solar System: Prospective Energy and Material Resources. Editors BadescuV.ZacnyK. (Cham, Switzerland: Springer), 47–165.

  • Google Scholar
  • View reference in article

• 32

  TarauC.LeeK.-L.LeithererB.ChettyK.LutzA.RokkamS. (2022). “Status of development of a thermal probe for icy planet exploration - II,” in Presented at the Nuclear Emerging Technol. Space (NETS) Conf. Cleveland, OH, USA.

  • Google Scholar
  • View reference in article

• 33

  TonerJ. D.CatlingD. C.LightB. (2015). A revised pitzer model for low-temperature soluble salt assemblages at the phoenix site, Mars. Geochim. Cosmochim. Acta166, 327–343. 10.1016/j.gca.2015.06.011

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  VollerV. R.PrakashC. (1987). A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems. Int. J. Heat. Mass Transf.30 (8), 1709–1719. 10.1016/0017-9310(87)90317-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  WakitaS.JohnsonB. C.SilberE. A.SingerK. N. (2024). Multiring basin formation constrains Europa's ice shell thickness. Sci. Adv.10 (8), eadj8455. 10.1126/sciadv.adj8455

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  WeissP.YungK. L.KömleN.KoS. M.KaufmannE.KarglG. (2011). Thermal drill sampling system onboard high-velocity impactors for exploring the subsurface of europa. Adv. Space Res.48 (4), 743–754. 10.1016/j.asr.2010.01.015

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  WeissP.YungK. L.NgT. C.KömleN.KarglG.KaufmannE. (2008). Study of a thermal drill head for the exploration of subsurface planetary ice layers. Planet. Space Sci.56 (9), 1280–1292. 10.1016/j.pss.2008.04.004

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  ZacnyK.PaulsenG.McKayC. P.GlassB.DavéA.DavilaA. F.et al (2013). Reaching 1 m deep on Mars: the icebreaker drill. Astrobiology13 (12), 1166–1198. 10.1089/ast.2013.1038

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  ZacnyK.MuellerJ.CostaT.CwikT.GrayA.ZimmermanW. (2018). “SLUSH: europa hybrid deep drill,” in Proc. 2018 IEEE Aerosp. Conf. Big Sky (MT, USA), 1–14.

  • Google Scholar
  • View reference in article

• 40

  ZacnyK.StolovL.MellerowiczB.PalmowskiJ.BywatersK.WareL.et al (2019). SLUSH: search for life using submersible heated drill. LPI Contrib.2100, 3013.

  • Google Scholar
  • View reference in article

• 41

  ZelenyiL.KorablevO.MartynovM.PopovG.BlancM.LebretonJ.et al (2011). Europa lander mission and the context of international cooperation. Adv. Space Res.48 (4), 615–628. 10.1016/j.asr.2010.11.027

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  ZolotovM. Y.ShockE. L. (2001). Composition and stability of salts on the surface of europa and their Oceanic origin. J. Geophys. Res.106 (E12), 32815–32827. 10.1029/2000je001413

  • CrossRef
  • Google Scholar
  • View reference in article