EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 711 (July 2026) A&A, 711 (2026) A178 Full HTML
Open Access
Issue
A&A
Volume 711, July 2026
Article Number A178
Number of page(s) 15
Section Planets, planetary systems, and small bodies
DOI https://doi.org/10.1051/0004-6361/202554392
Published online 13 July 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Since the discovery of the TRAPPIST-1 planet system (Gillon et al. 2017), the system has provided an intriguing opportunity to the exoplanet community for observational and theoretical studies (see e.g., Krissansen-Totton & Fortney 2022; Zieba et al. 2023; Pichierri et al. 2024; Schoonenberg, Djoeke et al. 2019). TRAPPIST-1 hosts seven known, transiting exoplanets in orbital resonance, and in almost edge-on orbits (Teyssandier et al. 2022; Agol et al. 2021). It is the first known system of potentially rocky, terrestrial exoplanets with incident fluxes in the range of the solar system rocky planets, three of which are located within the star’s habitable zone.

The small radius of the star aids the planets’ characterization. Agol et al. (2021) performed an analysis of transit timing variations (TTVs) that includes all transit data from Spitzer since the discovery of the system, and provided precise mass and radius measurements for the planets. In addition, they observed a slight downward trend of density with orbital distance. One possible explanation is a different volatile or water content for the seven planets. If confirmed, the three outer planets, f, g, and h, might have between 4.5 and 6.4 wt% water assuming an Earth-like core (Agol et al. 2021). The latter, together with their low equilibrium temperatures, makes them good candidates for what Quick et al. (2020) define as “cold ocean planets”.

Cold ocean planets are ocean-bearing worlds with ice covered surfaces (Quick et al. 2020). Their interior structure and properties are expected to be similar to those of the solar system icy moons (Ehrenreich et al. 2006). As such, they might harbor subsurface oceans fueled by tidal or radiogenic heating, underneath an icy shell (Nimmo & Pappalardo 2016). Cold ocean planets might be in abundance in the galaxy (Zeng et al. 2019), and they are significant targets for habitability and astro-biology studies. Several candidate worlds of this nature have been proposed. For example, Quick et al. (2023) identified 17 planet candidates that might be icy, including TRAPPIST-1f, g and h, according to their measured densities and equilibrium temperatures. Another planet whose properties are compatible with a water-rich planet covered in an icy surface is LHS 1140 b (Quick et al. 2020, 2023; Cadieux et al. 2024).

Several studies have characterized the possible water-bearing interior of TRAPPIST-1f, g and h planets (see e.g., Barr et al. 2018; Dobos et al. 2019). Acuña et al. (2021) studied their putative hydrospheres by modeling their pressure-temperature profiles, and they conclude that water on the planets would be in solid phases. The planets’ nonzero orbital eccentricity and short orbital periods suggest that tidal heating might be significant in their thermal budget (Luger et al. 2017). Barr et al. (2018) studied tidal heating in the -1 system, and they conclude that the interior of the planets can have similar geodynamics as solar system icy moons.

Cryovolcanism is the eruption or extrusion of volatile materials, such as water, aqueous solutions, or other ices, in liquid or vapor form that would otherwise be frozen at the surface of icy planetary bodies where temperatures fall below the desorption temperature of most volatiles (Geissler 2015). If the outer TRAPPIST-1 planets have similar geological features to the solar system icy moons, they might also exhibit cryovolcanism. Cryovolcanism has been detected on Enceladus (Porco et al. 2006; Dougherty et al. 2006), Triton (Croft et al. 1995), and other icy bodies in the Solar System (Fagents et al. 2022) and it is a source of water outgassing that contributes to the formation of exospheres, among other processes such as surface sputtering. In addition, cryovolcanic activity has long been hypothesized for Europa (Fagents et al. 2000; Fagents 2003). HST Ly-α observations (Roth et al. 2014) and Galileo magnetic and plasma wave signatures (Jia et al. 2018) were previously interpreted as evidence for plume activity. A more recent analysis, however, suggests that the HST Ly-α observations are more consistent with an atomic hydrogen exosphere than active plumes (Roth et al. 2026). In the Solar System, several icy moons possess tenuous exospheres or outgassing activity, such as Europa, which has an oxygen-dominated exosphere (Smyth & Marconi 2006; Oza et al. 2019b), and Enceladus, which possesses plumes of water vapor and ice grains, containing traces of salts, silicates, and organic compounds (Postberg et al. 2011; Porco et al. 2006). These plumes feed the E ring, a diffuse torus around Saturn which extends several satellite radii from Enceladus (Baum et al. 1981).

Our knowledge of the atmospheres of the TRAPPIST-1 planets is limited, including whether or not they possess atmospheres in the first place. The extended pre-main-sequence (PMS) phase of M-dwarf stars, like TRAPPIST-1, makes the planets subject to atmospheric loss due to their star’s prolonged period of high luminosity and XUV emission (Luger & Barnes 2015). Nevertheless, efforts have been made to observationally characterize the TRAPPIST-1 atmospheres; de Wit et al. (2018) ruled out cloudless, hydrogen-dominated atmospheres for planets d, e and f. In addition, Zieba et al. (2023) found evidence of either tenuous or no atmosphere for planet c, and the photometric secondary eclipse observations of Greene et al. (2023) hint that planet b is most likely a bare rock. Motivated by these results, in this work we adopt a simplified approach in which the outer TRAPPIST-1 planets are treated as essentially bare icy bodies, and we model only a tenuous, cryovolcanically supplied water exosphere, neglecting any additional background atmosphere.

Tidal models of varying complexity have been applied to study bodies and moons within the Solar System (Rovira-Navarro et al. 2022; Bierson & Steinbrügge 2021) and exoplanet systems (Henning & Hurford 2014; Shoji & Kurita 2014; van Woerkom & Kleisioti 2024). They range from models that assume tidal parameters that are independent of the orbital period and interior structure (constant-Q approach), and more complicated multilayered approaches, which take into account the viscoelastic behavior of each layer material properties (Kleisioti et al. 2023; Rovira-Navarro et al. 2021). Since the discovery of the TRAPPIST planets, several authors have attempted to constrain their interior structure and potential dynamics assuming that the outer planets TRAPPIST-1f, g and h, are cold ocean worlds. Quick et al. (2023) obtained ice shell thicknesses, tidal heating rates and cryovolcanic column densities. They used the constant-Q approach to assess tidal heating, which did not take the dependency of the tidal parameters on each of the planets’ specific interior structure and orbit into account. Barr et al. (2018) also modeled tidal heating in the system. However, they did not use the updated mass and radius constraints of Agol et al. (2021), resulting in significantly higher predicted water contents.

In this work, we examine the interior structure compatible with the updated density estimates using a multilayer, viscoelastic approach and the potential for cryovolcanism and its detectability. By constraining the heat generated within each layer, we determined the range of internal structures compatible with a given ice shell thickness and assessed the total internal heat budget. We focused on the possibility that the outer planets TRAPPIST-1f, g and h contain water and do not possess thick atmospheres, with the goal of assessing possible thermal equilibrium states, and whether exospheres, similar to the solar system icy moons, can be detectable with JWST. We modeled fully differentiated planets with an icy shell, sub-surface ocean, high pressure ice layer (HPI), mantle, and a core, and obtained interior structures in thermal equilibrium. We studied the possible thermal states of an icy shell and a subsurface ocean on TRAPPIST-1f, g and h, using thermal models that have been applied to solar system icy satellites. To do so, we used the planet masses and radii of Agol et al. (2021), together with a tidal model that captures the behavior of material properties as a function of the orbital period. In Section 2, we describe our interior, thermal, and tidal models. In Section 3, we present the partitioning of tidal dissipation within the planet interior, and the resulting interior constraints for TRAPPIST-1f. Sections 4 and 5 focus on the resulting cryovolcanic activity and its predicted transmission spectroscopy signatures. Finally, Section 6 compares our results with previous studies and discusses whether interior structures in thermal equilibrium can produce cryovolcanic signatures detectable with JWST.

2 Interior structure and tidal model

We assumed a fully differentiated body, with radially homogeneous layer properties that consist of an iron core, a silicate mantle, a high pressure ice (HPI) layer that accounts for the possible presence of the water ice polymorphs II through VII (e.g., Jaccard 1976), a liquid, subsurface ocean and an ice I shell with densities ρc, ρm, ρHPI, ρocean, and ρice,I (Dobos et al. 2019; Barr et al. 2018). For our work, we require a model to determine feasible interior structures of the planets that are compatible with the mass and radius values (see Table 2), under the assumption of thermal equilibrium. Here we describe the thermal model (Section 2.1) used to calculate the equilibrium states of the ice I shell, and the tidal model (Section 2.2) to determine the heating in the dissipative layers of the planets’ interior. We describe the planets’ thermal equilibrium model in Section 2.3.

2.1 Heat transfer model

For the heat transfer model of the ice I shell, we followed the approach of Sohl et al. (2003), who modeled the thermal structure of Titan’s ice shell. The ice I shell was divided into a conductive outer layer and a convective sub-layer, as is typically modeled for solar system icy moons (see e.g., Hussmann et al. 2002; Green et al. 2021). The heat that is transported through the conductive stagnant lid layer is

Qcond4πRP2=k(Ttop−TS)Dcond,Mathematical equation: ${{{Q_{cond}}} \over {4\pi R_P^2}} = k{{\left( {{T_{top}} - {T_S}} \right)} \over {{D_{cond}}}},$(1)

where RP is the planet’s radius, Dcond is the thickness of the conductive layer, TS and TTop are the surface and top temperature of the convective layer, and k is the thermal conductivity coefficient, assumed to be equal to 3.3 W m−1 K−1 (Spohn & Schubert 2003).

The equation for heat transport in the convective, sub-layer is

Qconv4π(RP−Dcond)2=k(Tm−Ttop)Dconv(RaRac)β.Mathematical equation: ${{{Q_{conv}}} \over {4\pi {{\left( {{R_P} - {D_{cond}}} \right)}^2}}} = k{{\left( {{T_m} - {T_{top}}} \right)} \over {{D_{conv}}}}{\left( {{{Ra} \over {R{a_c}}}} \right)^\beta }.$(2)

Where Dconv is the thickness of the convective layer, Tm is the melting temperature of the ice I, and Ra and Rac are the Rayleigh and the critical Rayleigh number (in the order of 103 (Schubert et al. 1979)). The Rayleigh number is given by Equation (3). Convection becomes a significant mechanism for heat transport in the ice shell when Ra ≥ Rac. If Ra ≤ Rac, we assume that the ice shell is purely conductive (Equation (1)).

Ra=αρice,Igp(Tm−Ttop)Dconv3κη(Tint),Mathematical equation: $Ra = {{\alpha {\rho _{ice,I}}{g_p}\left( {{T_m} - {T_{top}}} \right)D_{conv}^3} \over {\kappa \eta \left( {{T_{int}}} \right)}},$(3)

where α and κ are the ice I thermal expansion coefficient and thermal diffusivity equal to 1.6 × 10−4K−1 and 1.47 × 10−6m2s−1, respectively (Sohl et al. 2003). gP is the surface gravity of each planet and η(Tint) is the viscosity of the convective layer’s mean temperature (Tint).

In thermal equilibrium, the heat flux (Qeq) that is transferred through the conductive layer is equal to the flux convected through the sublayer:

Qcond=Qconv=Qeq.Mathematical equation: ${Q_{cond}} = {Q_{conv}} = {Q_{eq}}.$(4)

To obtain the thermal equilibrium state of the ice I shell, we solved Equations (1), (2), (3) and (4) iteratively for the layer thicknesses (Dstag, Dconv) and the equilibrium heat flux Qeq for different values of total ice I layer thickness (Dice). To do so, we firstly assumed a value for Dice, which is the free parameter in our interior models. This set the pressure, and thus, the melting temperature at the bottom of the convective layer (Tm), at the depth Rice,I as seen in Figure 1. The melting temperature was calculated as a function of pressure, as shown in Figure A.1 for an ammonia containing ocean of 5 wt.% (Grasset et al. 2000). These concentrations are consistent with those inferred for icy Solar System bodies, like Titan, where ammonia acts as an antifreeze facilitating the existence of subsurface oceans under cold conditions (Sohl et al. 2003). Furthermore, formation beyond the water ice line of TRAPPIST-1f, g and h (see Section 6.2) would promote the accretion of volatile compounds, including ammonia (Öberg et al. 2011). With the assumption that the ocean is isothermal, we calculated the pressure at the bottom of the ocean, and thus the thickness of the ocean, and the radius RHPI . This process is visualized in the red segment of Figure 2. For more details on the iterative method, we refer the reader to Sohl et al. (2003).

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Schematic overview of the interior model of TRAPPIST-1f, g and h. The ice I layer is divided into a conductive and a convective sublayer. The thermal model determines the thickness of the conductive and convective layers, Rice,I, RHPI, and Qeq. The tidal model computes the tidal heating rates in the HPI, and the mantle. The layer boundaries are not to scale. Note that the convective layer is only present for a sufficiently high (≳103) Rayleigh number.
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Schematic of the parameter calculation process. The red segment visualizes the heat transfer model in the outer ice shell, which is described in Section 2.1. The gray segment describes the process of calculating the thermal equilibrium in the entire planet (Section 2.3). The latter thermal equilibrium dictates that Qeq = Qtidal,HPI + Qrad.

2.2 Tidal response model

The amount of tidal heat dissipated in the interior of a planet is given by Equation (5) (e.g. Segatz et al. 1988).

E˙=−212Im(k2)(nR)5Ge2,Mathematical equation: $\dot E = - {{21} \over 2}{\rm{Im}}\left( {{k_2}} \right){{{{\left( {nR} \right)}^5}} \over G}{e^2},$(5)

where n is the planet’s mean motion, e the orbital eccentricity of the planet, and G the universal gravitational constant. We assumed that the planets are rotating synchronously (Lobo et al. 2023; Turbet et al. 2018). In addition, we assumed that the planets have zero obliquity, and therefore have zero obliquity tides.

Im(k2) is the imaginary part of the k2 Love number, which describes the deformation of the planet’s gravity field due to the tidal interactions with the host star. Im(k2) depends on the interior structure and properties (and hence the mass of the dissipative layers), rheology of the planet and its orbital period (Segatz et al. 1988; Renaud & Henning 2018; Shoji & Kurita 2014). Unlike studies that assume a constant k2/Q (where Q is the tidal quality factor) value for all planets (e.g., Quick et al. 2023), in our model (Im(k2) = k2/Q) is computed self-consistently and varies with the planet’s interior structure, rheology, and orbital period. We computed the tidal dissipation using the viscolestic theory for self gravitating bodies (see e.g., Peltier 1974; Sabadini 2016) with the Maxwell rheological law, which connects stress and strain, and has already been applied to study solar system moons (Hussmann & Spohn 2004; Fischer & Spohn 1990; Tobie et al. 2019; Roberts & Nimmo 2008; Rovira-Navarro et al. 2022) and exoplanet systems (Rovira-Navarro et al. 2021; Dobos & Turner 2015; Hay & Matsuyama 2019). We note that more complex rheological models, such as the Andrade rheological law, have also been used to investigate tidal dissipation in planetary interiors (see e.g., Kleisioti et al. 2023; Henning & Hurford 2014). By solving the equations of motion for the deformation of each layer in the Fourier domain via the correspondence principle (Peltier 1974), we obtained the gravitational potential at the surface and thus, the Love number. We use the method described by Rovira-Navarro et al. (2024) for spherically symmetric bodies, and used the the LOV3D software package (Rovira-Navarro et al. 2024) to compute the tidal dissipation in each layer for a given planetary structure.

The rheological and physical properties that were used for the thermal and tidal models are presented in Table 1, while the orbital properties of each planet are shown in Table 2. We computed tidal dissipation in the mantle, HPI layer, and ice I shell for each interior structure. For the equilibrium heat-balance closure in Eq. (6), we initially neglected the mantle and ice I contributions and subsequently verified that this approximation is valid over the explored parameter space (Section 3.1).

For the tidal response model, we assumed an Earth-like mantle and core composition (iron-rich). The values for mantle viscosity, can range from 1015 to 1022 Pa s, with 1016 Pa s representing partially molten cases (Barr et al. 2018). A shear modulus of 50 GPa is commonly adopted for the mantle at temperatures below the solidus (Hussmann et al. 2002; Barr et al. 2018). Viscosities of the high-pressure ice (HPI) layer, range from 1014 to 1020 Pa s (Barr et al. 2018), with 1014 being close the melting point. The HPI shear modulus was taken to be μΗΡΡ ≈ 3.5× 1010 Pa, assumed independent of temperature, as in Barr et al. (2018). Viscosities of Ice I range between 1013 and 1016 Pa s, and shear modulus estimates for Ice I typically range from 3.3 to 4.8 GPa (Hussmann et al. 2002; Barr et al. 2018).

2.3 Interior heat balance and structure calculation

The heat transfer model described in Section 2.1 provides the heat flux for a given value of the ice I layer thickness (Dice), and the associated thermal equilibrium heat output (Qeq). In this manner, Qeq is a function of Dice. For the entire planet to exist in equilibrium, Qeq should be equal to the amount of heat that is produced in the interior. Our aim was to determine the mass and thickness of each interior layer – the core (mc), mantle (mm), and high-pressure ice (HPI) layer (mHPI) – for a range of ice shell thicknesses (Dice), such that the total internal heat production is consistent with Qeq and measured planetary mass and radius.

We focused on two dominant sources of interior heat: tidal heating in the HPI layer and radiogenic heating in the silicate mantle, and assumed that the heat flux through the ice shell (Qeq) must equal the sum of the radioactive decay (Qrad) that is produced by the accumulation of radioactive elements in the silicate mantle, and the tidal heat that is produced in the high pressure ice layer (Qtidal,HPI). Equation (6) was used as a closure condition to identify interior structures in thermal equilibrium. We then verified its validity by recomputing tidal dissipation in all layers (mantle, HPI, and ice I) for the derived structures and by exploring a wide range of rheological parameters via a Monte Carlo sensitivity analysis (Section 2.4).

Qeq=Qtidal,HPI+Qrad.Mathematical equation: ${Q_{{\rm{eq}}}} = {Q_{{\rm{tidal,HPI}}}} + {Q_{{\rm{rad}}}}.$(6)

For the calculation of Qrad, we followed the approach of Schubert et al. (1986), and used Equation (1) from Hussmann & Spohn (2004), which expresses the radiogenic heat production as a function of mantle mass (mm). We assumed a mantle with chondritic composition and a system age of 7.6 Gyr (Burgasser & Mamajek 2017). We calculated Qtid,HPI using Equation (5), which provides a direct relationship between Qtid,HPI and the mass of the high-pressure ice layer (mHPI). We note that additional mechanisms, such as hemispheric tectonic driven by strong day-night temperature contrasts on exoplanets, may also influence heat transport (Meier et al. 2021).

To constrain the full interior structure, we applied the following constraints, which describe the mass conservation (Equation (7)), and the geometry of the system (Equation (8)), while assuming homogeneous layers with constant density.

Mp=mc+mm+mocean+mice,I+mHPI.Mathematical equation: ${M_p} = {m_c} + {m_m} + {m_{ocean}} + {m_{ice,I}} + {m_{HPI}}.$(7)

We solved for RHPI, defined as the radius of the boundary between the high-pressure ice layer and the ocean, using the volumes of the layers interior to this interface:

RHPI=(34π)13(mHPIρHPI+mmρm+mcρcore)13.Mathematical equation: ${R_{HPI}} = {\left( {{3 \over {4\pi }}} \right)^{{\textstyle{1 \over 3}}}}{\left( {{{{m_{HPI}}} \over {{\rho _{HPI}}}} + {{{m_m}} \over {{\rho _m}}} + {{{m_c}} \over {{\rho _{core}}}}} \right)^{{\textstyle{1 \over 3}}}}.$(8)

Where mice,I and mocean is the mass of the ice I shell, and the subsurface ocean. The rest of parameters in Equations (7), (8) are described in Table 1.

Since Qtid,HPI is a function of mHPI (see Section 2.2) and Qrad is a function of mm, we solved the system of Equations (6), (7), (8) for mm, mHPI, and mc for a range of Dice. Thus, we calculated mc, mm, mHPI, Qrad, and Qtidal,HPI as a function of Dice. This allowed us to derive the Core Mass Fractions (CMF) and Water Mass Fractions (WMF) values as a function of Dice. The gray boxed region of Figure 2 visualizes the calculation process described in this Section.

We finally assessed whether the assumption made in Equation (6) was valid. We recomputed the expected tidal heating rates in the different layers, as seen in Figure 3 for the computed interior structures with varying HPI, mantle, and core masses. In all cases, we concluded that the tidal heating rate in the mantle (Qtidal,m) is at least one order of magnitude lower than both the tidal heating in the HPI and the expected radiogenic heat (see Figure 3), while the tidal heat in the ice I shell (Qtidal,iceI) is typically four orders of magnitude smaller than the total heat generated in the interior (see Figure B.1), validating our initial assumption of neglecting Qtidal,m and Qtidal,iceI in Equation (6).

Table 1

Average physical and rheological properties of the interior model layers.

Table 2

Planet parameters used in the model.

Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Comparison of the heat contribution from several sources for TRAPPIST-1f as a function of ice shell thickness. Tidal heat in the mantle is negligible in the entire parameter space, while radiogenic heat and tidal heat in the ice shell both contribute significantly to the ice sheet equilibrium heat. TRAPPIST-1f due to its size and eccentricity has the highest tidal heat of the three planets.

2.4 Sensitivity analysis

To evaluate the robustness of our tidal dissipation results and address uncertainties in interior rheology and internal heating, we performed a Monte Carlo sensitivity analysis over a broad range of model parameters. Layer viscosities were varied independently using log-uniform distributions spanning several orders of magnitude. Mantle viscosities were sampled over the range log10(ηm/Pas) = 18–21, HPI viscosities over log10(ηHPI/Pas) = 14–18, and ice I shell viscosities over log10(ηice I/Pas) = 13–16. These ranges represent values commonly adopted in previous studies of icy bodies (Barr et al. 2018; Hussmann et al. 2002). In addition, the radiogenic heat flux was varied between half and twice the nominal value used in the reference model, corresponding to H = (0.25–1.0) × 10−12 W m−2. For each realization, we computed the tidal dissipation rates in the mantle, HPI layer, and ice I shell, as well as the total tidal dissipation Qtidal and the equilibrium heat flux Qeq.

3 Tidal dissipation and interior constraints for TRAPPIST-1f, g and h

We first demonstrate that neglecting tidal dissipation in the silicate mantle and the ice I shell does not significantly affect the equilibrium heat balance over the parameter space considered here, by quantifying how tidal dissipation is partitioned among the planet’s interior layers. We then assess which layers dominate the tidal energy budget and present the interior structures that can exist in thermal equilibrium.

3.1 Sensitivity analysis results

Figure 4 shows the fractional contributions of tidal dissipation in the mantle, high-pressure ice layer, and outer ice I shell, normalized by the total tidal dissipation, as a function of ice shell thickness Dice. Across the entire parameter space explored, tidal dissipation within the ice I shell accounts for only a small fraction of the global energy budget. In all realizations, the maximum fraction of tidal dissipation occurring in the ice I shell is fice = Qtidal,ice I/Qtidal ≲ 0.1. Because radiogenic heating is a significant fraction of the total energy budget, the relative contribution of ice I tidal dissipation to the total heat budget is even smaller, yielding Qtidal,ice I /Qtotal ≲ 0.04. The largest fractional contributions from the ice I shell occur primarily for models with relatively stiff HPI layers (see Figure 5 which illustrates the dependence of tidal dissipation partitioning on layer viscosities). In these cases, reduced deformation in the HPI layer allows a greater fraction of the tidal strain to be accommodated within the overlying ice I shell. Nevertheless, even under these conditions, tidal dissipation in the ice I shell is a small portion of the total heat budget.

Enhanced tidal dissipation in the mantle occurs only in a subset of models with small core mass fractions (CMF ≲ 0.1) and relatively low mantle viscosities. In these cases, the mantle contribution can reach Qtidal,m/Qeq ~ 0.2, since radiogenic heating in the mantle also contributes significantly to the total heat budget. Such conditions are mostly achieved for stiff HPI layers combined with smaller mantle viscosities (Figure 5). If partial melting were to occur in the mantle, the fraction of tidal dissipation localized there would likely increase further.

Overall, these results demonstrate that, for the interior structures considered here, tidal dissipation is not dominated by the ice I shell, in contrast to inferences drawn from Solar System satellites (see e.g., Sohl et al. 2003). This difference arises primarily from the specific interior configurations explored in this study, in particular the presence of very thick high-pressure ice (HPI) layers, in some cases exceeding ~2000 km in thickness. Such a thick viscoelastic HPI layer represents a large volume of deformable solid and can therefore contribute substantially to the total tidal dissipation budget, reducing the relative importance of the overlying ice I shell.

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Fractions of QtidalQeqMathematical equation: ${{{Q_{{\rm{tidal}}}}} \over {{Q_{{\rm{eq}}}}}}$ of the mantle, high-pressure ice layer, and outer ice I shell, as a function of Dice. The lighter and darker shaded regions indicate the 75% and 95% uncertainty intervals derived from the Monte Carlo distribution, respectively. Note that the y-axis scales differ between panels.

3.2 Possible interior structures of the planets TRAPPIST-1f, g and h

In Figure 6a, we present the interior structure of TRAPPIST-1f (Rcore, Rmantle, RHPI) as a function of Dice (see Section 2.3). Figures C.1, C.2 present the corresponding Figures for TRAPPIST-1g and h. A given Dice thickness can transport an amount of heat equal to the radiogenic heat generated in the mantle and the tidal heat produced in the HPI layer. Thicker ice I layers can transport smaller amounts of heat. The latter translates to a structure with higher MHPIMMantleMathematical equation: ${{{M_{HPI}}} \over {{M_{Mantle}}}}$ ratios, as seen in Figure 3. Thus, as Dice increases, the ratio MHPIMMantleMathematical equation: ${{{M_{HPI}}} \over {{M_{Mantle}}}}$ also increases (Figure 6). The various interior structures are consistent with the observed densities of the planets. As a result, the calculated range for ice shell thickness is limited; thicker ice shells reach thermal equilibrium at lower Qeq values. Achieving such equilibrium requires a larger core, since a smaller mantle reduces radiogenic heating. This, in turn, prevents the total internal heat, radiogenic from the mantle and tidal from the high-pressure ice (HPI) layer, from exceeding the heat transfer capacity of the ice shell. Thus, the high-pressure ice (HPI) layer expands rapidly with Dice, compensating for the increased core size.

The right horizontal-axis limit of Figure 6 corresponds to the case of a planet without a mantle layer, meaning that the iron core would be in direct contact with the HPI (CMF = 0.75). The left horizontal-axis limit is consistent with an interior structure that lacks a HPI layer. Such an interior structure corresponds to the minimum WMF. Smaller WMF’s would not be compatible with TRAPPIST-1f, g and h having an icy shell with a subsurface liquid ocean under the assumptions of our model. In particular, changing the assumptions for the physical and rheological parameters in Fig. 1 will change the specific numerical values we obtain, but the trends and overall magnitude of the quantities we present will be representative for a body with the interior layer structure as defined in our model. Without an HPI layer, the subsurface ocean would directly interact with the mantle. This interaction is crucial for astrobiology and habitability studies, as it allows essential elements for the development of life to potentially enrich the subsurface ocean, similar to what is believed to occur on Europa (Hand et al. 2022). Recently (Hernandez et al. 2022) suggested that even when high-pressure ice layers are present, these do not fully block chemical communication. This suggests that interior-ocean exchange of electrolytes and essential elements may remain possible across a wide range of planetary conditions.

Figure 6b shows the thickness of the stagnant lid, convective layer, and ocean layer, as a function of Dice. If TRAPPIST-1f has a CMF of 0.5 and an ice shell in thermal equilibrium, the ocean, approximately 100 km thick, would lie about 6 km beneath an almost entirely conductive ice shell. The extremely thin convective shell thickness supports the assumption of a conducting ice I layer made in this work (Section 2.2), in contrast to the convective, thicker ice shells observed on icy moons within the solar system (see e.g., Hussmann et al. 2002).

We present the thermal equilibrium interior structure parameters for TRAPPIST-1f, g and h in Tables 3, 4, 5 respectively for CMF=0.18, 0.32, and 0.5 assuming that the planets have not been exposed to events that could have stripped away their mantle, such as mantle evaporation (Cameron 1985). For the planets TRAPPIST-1f and g, even in the case of a large CMF equal to 0.5, a subsurface ocean would exist below an ice shell of thickness ≲ 10 km. Such a thin ice I shell hints that the ocean is more likely to be in direct contact with the surface through cryovolcanism. As a result, a plume would likely be composed of material originating directly from the subsurface ocean, similarly to the cryovolcanic mechanism on Enceladus (Quick et al. 2023; Porco et al. 2014). TRAPPIST-1h, the least massive of the three planets, exhibits equilibrium interior structures with an ice layer thickness ≳10 km. This raises the possibility that any cryovolcanism on TRAPPIST-1h could resemble that of Europa, coming from water-pockets in the ice (Quick et al. 2023). The ice shell I thicknesses that we obtain are mostly in agreement with the results of (Quick et al. 2023). For some interior configurations, particularly those with thicker Ice I shells, our modeled tidal heating rates for TRAPPIST-1f (≈10 TW) fall within the same order of magnitude as previous estimates by Bolmont et al. (2026), who constrain tidal heating in TRAPPIST-1f between ≈10–30 TW. In the following section we discuss whether cryovolcanism could be observed with spectroscopic observations on TRAPPIST-1f, g and h.

Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Top: fice as a function of different layers’ viscosity. Bottom: fmantle as a function of different layers’ viscosity. Enhanced tidal dissipation occurs in the ice I shell for stiffer HPI layers. Tidal dissipation in the mantle contributes more to the heat budget for stiffer HPI layers and lower mantle viscosity.
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Interior structure as a function of Dice for the planet TRAPPIST-1f. (a) Radii of different layers as a function of Dice. (b) Thickness of the stagnant lid, convective and ocean layers as a function of Dice.

4 Modeling observations of cryovolcanic activity

Cryovolcanic activity has been observed on Saturn’s moon, Enceladus (Porco et al. 2006; Spencer et al. 2018), and is also suspected on Jupiter’s moon, Europa, based on multiple detections (Roth et al. 2014; Sparks et al. 2017; Paganini et al. 2020). In the case of Enceladus, a subsurface ocean is located ≲ 10 km beneath the south pole (Iess et al. 2014). The thin ice shell of Enceladus allows cracks to propagate fully downward to the ocean due to the long-term evolution of the ocean pressure and ice shell thickness (Rudolph et al. 2022). Cryovolcanic activity on Enceladus is thought to be continuous, while its intensity appears to be correlated with the orbital period (Ingersoll et al. 2020) and tidal stresses (Hedman et al. 2013; Soucek et al. 2024). On the other hand, cryovolcanic eruptions on Europa are sporadic (Sparks et al. 2017). In particular, (Paganini et al. 2020; Phillips et al. 2000) show that large cryovolcanic eruptions are likely rare, and Fagents et al. (2000) suggest that smaller eruptions are likely to occur more frequently. Given an ice shell of ≈25–30 km (Quick & Marsh 2015), cryovolcanic activity on Europa is thought to originate from water pockets that are located within the ice shell (Fagents 2003; Manga & Wang 2007).

Similar geological features, including volcanism and cryovolcanism, may be present on exoplanet systems (Kleisioti et al. 2024; Quick et al. 2023, 2020; Oza et al. 2019a; Meyer zu Westram et al. 2024; Guenther & Kislyakova 2019). However, predicting cryovolcanic outgassing on cold ocean exoplanets is challenging because it can be influenced by poorly constrained parameters and interior processes, including the total internal heat production, possibly the fraction of that heat generated by tidal dissipation, and the thickness of the outer ice shell. While surface fractures, like Enceladus’ tiger stripes, may form solely due to tidal deformation of the ice shell, radiogenic heating can also contribute to maintaining subsurface liquid reservoirs (Leone et al. 2021). If tidal dissipation is sufficient to create cracks, radiogenic heat may further sustain or enhance internal melting, and together these heat sources can maintain subsurface oceans that enable cryovolcanic activity (Ingersoll et al. 2020; Neveu & Rhoden 2019).

The total volcanic outgassing Mvolc that is driven by ĖTotal can be expressed in terms of the gravitational binding energy of a water molecule to the planet (U):

M˙Volc≈ηmH20E˙TotalU.Mathematical equation: ${{\dot M}_{Volc}} \approx \eta {m_{{H_2}0}}{{{{\dot E}_{Total}}} \over U}.$(9)

Where η is a factor describing the efficiency that the internal energy is used to power volcanic outgassing (Oza et al. 2019a), mH20Mathematical equation: ${m_{{H_2}0}}$ is the mass of one water molecule. The gravitational binding energy U is equal to GMpmH20RpMathematical equation: ${{G{M_p}{m_{{H_2}0}}} \over {{R_p}}}$ (Johnson et al. 2015). An efficiency factor of η = 1 would represent the theoretical upper limit, where all the available energy is used for the outgassing of water. In reality, this value is much lower. Using a mass loss rate of 300 kg s−1 (Villanueva et al. 2023) for Enceladus and an internal power of 35 GW (Hemingway & Mittal 2019), one obtains an efficiency factor for Enceladus on the order of 10−3.

Given their relatively close-in orbits and the active nature of TRAPPIST-1, exospheres may form not only through outgassing from cryovolcanic activity, but also via the mass loss through thermal escape caused by XUV activity (Wordsworth et al. 2018; Luger & Barnes 2015) (see also Section 6.2). Quick et al. (2023) modeled this latter phenomenon by calculating the water mass loss rate under the XUV flux of TRAPPIST-1, and translating it into column densities produced by the XUV-driven escape.

The total H2O column density was computed assuming a steady-state balance between continuous surface and interior supply and XUV-driven loss processes. Specifically, the observable column density is given by

Ntotal=M˙volcτXUV4πRp2mH2O+M˙subτXUV2πRp2mH2O,Mathematical equation: ${N_{{\rm{total}}}} = {{{{\dot M}_{{\rm{volc}}}}{\tau _{{\rm{XUV}}}}} \over {4\pi R_p^2{m_{{{\rm{H}}_{\rm{2}}}{\rm{O}}}}}} + {{{{\dot M}_{{\rm{sub}}}}{\tau _{{\rm{XUV}}}}} \over {2\pi R_p^2{m_{{{\rm{H}}_{\rm{2}}}{\rm{O}}}}}},$(10)

where Ṁvolc and Ṁsub are the mass rates of water vapor due to cryovolcanic outgassing and thermal sublimation, respectively.

In Eq. (10), volcanic outgassing is globally distributed while surface sublimation occurs only over the illuminated hemisphere. The quantity τXUV represents the effective molecular lifetime against XUV-driven photochemical loss.

We parameterized the sublimation flux, assuming it occurs primarily on the illuminated hemisphere and using the analytical form

M˙sub=2πasubTeq1/2exp(−bTeq),Mathematical equation: ${{\dot M}_{sub}} = 2\pi {a_{sub}}T_{eq}^{1/2}{\rm{exp}}\left( { - {b \over {{T_{eq}}}}} \right),$(11)

where asub = 1.9 × 1032cm−2 s−1 and b = 8500 K are empirically chosen coefficients to represent a low-sublimation case as was simulated for H2O at Ganymede’s leading hemisphere (Leblanc et al. 2017; Oza 2017). A second more realistic case, for pure water ice sublimation following Johnson et al. (1981) considers a canonical b = 6146K for the same asub. This substantially increases the expected column density of water vapor, allowing for a steady-state atmosphere to accumulate. Deviations from this column density in future observations are therefore constraints on the purity of ice on the surfaces of icy exoplanets in principle. We tested both sublimation cases and consider our low-sublimation results as representing a conservative, lower bound on the expected atmospheric column density from pure water ice sublimation.

Following Quick et al. (2023), we approximated this lifetime using the H2O photoionization timescale scaled from Europa, as described in Equation (12).

τXUV=4.54×107s(aaEuropa)2,Mathematical equation: ${\tau _{{\rm{XUV}}}} = 4.54 \times {10^7}{\rm{s}}{\left( {{a \over {{a_{{\rm{Europa}}}}}}} \right)^2},$(12)

where a is the planet’s semi-major axis and aEuropa is Europa’s orbital distance. This scaling assumes that the relevant ionizing flux varies as a−2 scaled to Europa’s reference value (Shematovich et al. 2005). Future works considering the bare rock or bare ice scenarios of airless exoplanets and exomoons should consider the detailed electron-impact dissociation, and plasma destruction rates of water as is estimated in 3-D exo-sphere models for Europa (Oza et al. 2019b).

We considered three limiting atmospheric configurations for a fixed volatile inventory ntot: (1) hydrostatic, representing the maximum possible collisional redistribution, (2) globally uniform sputtered exospheres, (3) and spatially localized sputtered exospheres, to explore how different distributions of cryovol-canically produced water affect observable transmission spectra. These scenarios span the range of plausible density structures that may arise from cryovolcanic activity and allowed us to assess how different distributions of cryovolcanic outgassed material influence detectability.

Table 3

Interior and thermal parameters of TRAPPIST-1f in thermal equilibrium for different values of CMF.

Table 4

Interior and thermal parameters of TRAPPIST-1g in thermal equilibrium for different values of CMF.

Table 5

Interior and thermal parameters of TRAPPIST-1h in thermal equilibrium for different values of CMF.

4.1 Hydrostatic atmospheres

We first examined the hydrostatic limit of an atmosphere composed exclusively of water vapor. Under this assumption, the total surface pressure reduces to the partial pressure of water,

Ptotal=NtotalmH2OgP,Mathematical equation: ${P_{{\rm{total}}}} = {N_{{\rm{total}}}}{m_{{{\rm{H}}_2}{\rm{O}}}}{g_P},$(13)

where gP is the planetary surface gravity.

Assuming an isothermal atmosphere in hydrostatic equilibrium, and an exponential decay model for the density of the planet’s atmosphere, we derived a density profile of the outgasses water. We set the atmospheric temperature equal to the planetary equilibrium temperature Teq (Table 2). This hydrostatic density structure was used to compute transmission spectra with the radiative transfer code petitRADTRANS (Mollière et al. 2019, 2020; Alei et al. 2022), employing the POKAZATEL linelist (Polyansky et al. 2018) for water opacities. Spectra were modeled over the wavelength range 0.2–2.5 μm, covering the dominant H2O absorption features, and were subsequently forwarded to the PandExo exposure time calculator (Batalha et al. 2017) to generate mock JWST data points using the JWST/NIRISS Substrip 256 instrument.

4.2 Globally uniform sputtered exospheres

We next considered the case in which the water molecules are uniformly distributed across the planetary surface, producing a globally homogeneous collisionless exosphere. The quantity Nvolc represents the vertical column density of H2O, i.e., the number of water molecules per unit surface area. To obtain the total number of outgassed water molecules, ntot, we assumed that the volatile material is uniformly distributed over the planetary surface. In this case, the relevant area is the full spherical surface, A=4πRp2Mathematical equation: $A = 4\pi R_p^2$, and the total number of molecules is

ntot=4πRp2Ntotal.Mathematical equation: ${n_{{\rm{tot}}}} = 4\pi R_p^2{N_{{\rm{total}}}}.$(14)

This ntot is then used as input to the radiative transfer tool PROMETHEUS (Gebek & Oza 2020) to generate synthetic transmission spectra of a water exosphere as follows.

PROMETHEUS computes transmission spectra for gaseous media in arbitrary geometries, assuming line-of-sight optical depth integration over a user-specified 3D density distribution. Because a cryovolcanic plume expands collisionlessly above the exobase and rapidly becomes non-hydrostatic, we treated the outgassed water as an extended exosphere rather than as an isothermal hydrostatic atmosphere. The density distribution was therefore modeled using a power-law scenario, which assumes a spherically symmetric exosphere of the form

n(r)=n0(Rpr)q,Mathematical equation: $n\left( r \right) = {n_0}{\left( {{{{R_p}} \over r}} \right)^q},$(15)

where q is the power-law index, r the radial distance from the planet’s center, and n0 is determined by the total number of absorbing molecules in the exosphere.

We adopted a power-law index q = 6 (Gebek & Oza 2020), which controls how steeply the exosphere density falls with radius, and a pseudo-temperature of 200 K for the H2O absorber, which determines the thermal Doppler width of individual lines. Steep density profiles with indices ≈6 are indicative of strong ion–neutral scattering, where neutral molecules interact efficiently with the surrounding plasma and are removed on short timescales (Johnson 1990). We computed the transmission spectra for the globally uniform sputtered exopsheres using the same opacity treatment, wavelength coverage, and JWST instrument simulation setup as the one adopted for the hydrostatic atmospheres in Section 4.1.

4.3 Effect of spatial localization

In the preceding analysis, we assumed that material outgassed via cryovolcanic activity is uniformly distributed over the planetary surface, corresponding to a globally averaged column density. This assumption provides a conservative lower limit on detectability. However, cryovolcanic activity are expected to be spatially structured, producing localized denser areas analogous to plume-like configurations. To quantify the effect of such localization, we introduced a surface filling factor f, defined as the fraction of the planetary surface area 4πRp2Mathematical equation: $4\pi R_p^2$ over which the sputtered material is distributed. For a fixed number of cryovolcanically produced molecules (nvolc), the effective column density associated with a localized cloud is then given by

Nvolc,loc=nvolc4πRp2f,Mathematical equation: ${N_{{\rm{volc,loc}}}} = {{{n_{{\rm{volc}}}}} \over {4\pi R_p^2f}},$(16)

such that f = 1 corresponds to a globally uniform cloud, while f ≪ 1 represents increasingly localized, plume-like distributions. The total effective column density becomes

Ntotal,loc=M˙volcτXUV14πRp2mH2O1f+M˙subτXUB2πRp2mH2O.Mathematical equation: ${N_{{\rm{total,loc}}}} = {{{{\dot M}_{ & {\rm{volc}}}}{\tau _{{\rm{XUV}}}}1} \over {4\pi R_p^2{m_{{{\rm{H}}_2}{\rm{O}}}}}}{1 \over f} + {{{{\dot M}_{{\rm{sub}}}}{\tau _{{\rm{XUB}}}}} \over {2\pi R_p^2{m_{{{\rm{H}}_2}{\rm{O}}}}}}.$(17)

Concentrating a fixed inventory into a smaller surface area increases the line-of-sight column density and therefore enhances the observability of spectral features. When f = 1 , the expression becomes equivalent to Equation (10), corresponding to the globally uniform redistribution case. Throughout this work, this limiting case was adopted as a conservative baseline. Physically plausible values of f are uncertain, since they depend on the number, geometry, and temporal variability of active vents. We adopted a value of f = 0.004 as a limiting case of a localized configuration corresponding approximately to the fractional surface area associated with an individual tiger-stripe source region (Porco et al. 2006). More generally, values of f ~ 10−2−10−1 may represent broad regional localization, whereas f ≲ 10−3 would correspond to highly localized vent-like sources.

5 Transmission spectroscopy signatures of cryovolcanic outgassing

Figure 7 presents model transmission spectra of cryovolcanic outgassing on TRAPPIST-1f, generated using the radiative transfer codes and exposure time calculator described in Section 4. The spectra are shown for three representative column densities and are compared to synthetic JWST/NIRISS data (black error bars). The top panel illustrates a tenuous hydrostatic atmosphere and serves as a detectability benchmark, while the middle and bottom panels correspond to a non-hydrostatic, globally uniform sputtered cloud for different column densities and number of transits. If JWST were to observe every transit of TRAPPIST-1f (P = 9.2 days) over a 12-year mission lifetime, corresponding to approximately 450 transits, a hydrostatic water column density of Nvolc ≃ 1.9 × 1019 molecules m−2 would become marginally detectable, as the 1.4 μm water feature exceeds the synthetic JWST uncertainties (Figure 7a).

To evaluate whether the column densities inferred from the outgassing models presented in Section 4 are physically plausible, Figure 8 shows the total water column density, Ntotal, as a function of the volcanic mass loss rate, Ṁvolc, and the corresponding efficiency factor η for TRAPPIST-1f. Vertical dashed lines indicate representative volcanic mass loss rates inferred for Enceladus, Europa, Earth, and Io. The light gray shaded region indicates the column densities required for detectability based on hydrostatic atmosphere spectra (Figure 7a), while the darker gray region corresponds to the detectability benchmarks derived from non-hydrostatic, sputtered-cloud spectra. For a globally uniform sputtered cloud, detectability requires a substantially larger total volatile inventory than the hydrostatic benchmark. The darkest shaded region represents the low-sublimation benchmark case for Ntot = 1.1 × 1024, molecules m−2, and for which detection could be achieved with approximately 20 transits (Figure 7b).

The top panel of Figure 8 corresponds to a globally uniform Ntotal, which represents a conservative detectability scenario for non-hydrostatic exospheres (f = 1 ). Even when assuming an observing campaign of 450 transits with JWST, the column densities required for detection fall within a regime that demands efficiencies η > 1, indicating that a spatially uniform sputtered cloud is energetically unphysical under the assumptions of our low-sublimation model (dashed blue line).The shaded red region highlights this unphysical parameter space. In contrast, the high-sublimation case (continuous blue line) lies well above the JWST detectability threshold corresponding to ~20 transits. This suggests that the signal could be detectable even under the conservative assumption of a globally uniform exosphere (f=1) and for physically realistic efficiency factors. A detection in this regime could thus provide quantitative constraints on the sublimation rate.

The bottom panel shows the same models under the assumption that the sputtered material is spatially confined to a surface fraction f = 0.004, yielding plume-like configurations that enhance the effective line-of-sight column density (Equation (17)). Such localization allows the low-sublimation modeled columns to approach values higher than the hydrostatic detectability benchmark. Because these configurations are still modeled as non-hydrostatic, collisionless exospheres, their detectability is more consistently assessed relative to the sputtered-cloud benchmark. When evaluated against this more appropriate benchmark, localized sputtering represents a substantially more favorable scenario than a globally uniform cloud. For Ṁvolc between Earth-like and Io-like values, a region of observability for low-sublimation emerges that requires 20 transits, and depends on the efficiency factor η, which quantifies how effectively internal energy is converted into volcanic outgassing and remains poorly constrained. High sublimation column densities, on the other hand, could be observed for both uniformly distributed and localized sputtered clouds with ≈20 transits. Any detection at column densities that would otherwise require η > 1 therefore, would imply that the absorbing material is spatially confined rather than globally distributed, consistent with localized sputtering or plume-like activity or high water sublimation.

Taken together, these results show that (i) hydrostatic water atmospheres provide an upper bound on the detectability of cryovolcanic water vapor, (ii) globally uniform sputtered exospheres with low sublimation rates are energetically inconsistent with detectable transmission signals, and (iii) spatial localization of outgassed material reduces the energetic requirements of a detection (iv) a detection would provide quantitative constraints on the sublimation rate of an icy exoplanet. These findings place constraints on the detectability of cryovolcanic water vapor on TRAPPIST-1f with JWST and motivate future observations with next-generation facilities.

Cryovolcanic and sputtered neutral environments in the Solar System provide useful context for interpreting the column densities and geometries considered here. In a simple way, Europa hosts either a sputter-dominated water exosphere or a sublimation-dominated exosphere in which neutrals follow largely ballistic trajectories, producing extended clouds of neutral water dissociation products with column densities of order 1012–1014molecules cm−2, depending on the orbital phase (Oza et al. 2019b). For instance, it is well known the darker trailing hemisphere of Europa is far hotter and also coincident with plasma ram (Roth 2021). Ganymede, by contrast, possesses an intrinsic magnetosphere that enhances ion-neutral interactions resulting in more compact distributions following Ganymede’s field lines (Leblanc et al. 2023). Shallow profiles (q ~ 3) imply extended, Europa-like distributions, whereas steeper profiles (q ~ 6), as used in this work, confine the same inventory closer to the planet, more analogous to plasma-modified environments such Ganymede’s magnetic exosphere.

For comparison with Solar System literature, the water geysers at Enceladus were probed by JWST to have a total of 𝒩tot ~ 1034.9 water molecules (Villanueva et al. 2023), roughly a third which begin orbiting Saturn (Cassidy & Johnson 2010). In the scenario of a high fraction of atmospheric escape then, a circumstellar water torus may also be amenable for detection about the ultra-cool dwarf star. A detection of a tenuous circumstellar torus would be indicative of a cryovolcanic system venting directly to space (Villanueva et al. 2023). A larger, total inventory of ntot ~ 1038 water molecules corresponds to an average stellar-disk line-of-sight (LOS) column density of 〈N〉 ≈ 1 × 1015.6molecules cm−2 when normalized by the stellar disk area, following the definition of 〈N〉 used in Oza et al. (2019a). This value is comparable to columns discussed for neutral column densities at Europa (Oza et al. 2019b; Roth et al. 2014) and Ganymede’s water exosphere (Leblanc et al. 2023).

Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Synthetic transmission spectra of water atmospheres and exo-spheres on TRAPPIST-1f, taking into account cryovolcanic outgassing, sublimation, and volatile loss from XUV irradiation. The black data points are generated with PandExo for JWST/NIRISS. (a) Spectrum of a tenuous hydrostatic water atmosphere for Ntotal = 1.9 × 1019 (hydrostatic benchmark) produced using petitRADTRANS (purple line). Mock JWST data points assume 450 transits and a binned resolution of 12. (b) Spectrum of a uniformly distributed sputtered cloud produced using PROMETHEUS (blue line). The darker blue data point are binned observations, demonstrating how the signal at ≈1.9 microns gets diluted. This sub-figure corresponds to a Ntotal = 1.1 × 1024 molecules m−2 (low sublimation case). Mock JWST data points assume 20 transits. (c) Same as (b) but for Ntotal = 1.6 × 1026 molecules m−2 which corresponds to the benchmark of the high sublimation case. Mock JWST data points assume 15 transits.
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Ntotal as a function of Ṁvolc for TRAPPIST-1f. The shaded red region indicates parameter space requiring efficiencies η > 1, which are unphysical. The hydrostatic and sputtered cloud benchmarks represent the column densities observable with JWST/NIRISS with 450 transits for a hydrostatic atmosphere and a collisionless exosphere, respectively. The darkest gray region indicates an exosphere’s column density observable with 20 transits. Top: detectability of water column densities for a globally uniform sputtered cloud. Bottom: detectability of water column densities assuming the sputtered material is confined to a surface fraction f = 0.004, yielding enhanced line-of-sight columns, which represent plume-like configurations.

6 Discussion

In this section, we place our modeling results in the context of previous studies, discuss the detectability of cryovolcanic activity on TRAPPIST-1f, g and h, and outline directions for future observational and theoretical work.

6.1 Water mass fraction comparison with previous work

In addition to refining the mass and radius measurements for the TRAPPIST-1 planets, Agol et al. (2021) theorized over their possible interior structures and water content. They estimated their WMFs assuming a condensed water layer of 300 K at 1 bar. Acuña et al. (2021) also modeled the interior and hydrospheres of the Trappist-1 planets, taking into account Equations of State (EOS) for different phases of water, including water vapor and high pressure ice phases. In Figure 9, we show the calculated WMFs for CMF = 0.25 from our model, taking into account the mass and radius uncertainties (Table 2), and we compare them to those derived by Agol et al. (2021); Acuña et al. (2021). We compare our results solely with the latter two studies, as previous work (Barr et al. 2018; Dorn et al. 2018) on the interior structures of the system have not used the updated values for mass and radius. In particular, they have modeled the interior of the TRAPPIST-1 planets using mass and radius measurements, constrained by Delrez et al. (2018); Grimm et al. (2018). Our results (Figure 9) are mostly within the uncertainties of Acuña et al. (2021), and there is a larger discrepancy with the results of Agol et al. (2021), whose results are significantly lower (for planets f and h). This discrepancy can be explained by our different modeling strategies. In order to reach the given density constraints, dictated by the mass and radius measurements, Agol et al. (2021) used a condensed water layer, without taking into account any equation of state for different phases, thus assuming lower density for the water layer. The agreement of our results with Acuña et al. (2021) is to be expected, as we included similar phases for water (high pressure ice polymorphs). Our results can be interpreted as one out of the multiple interior structures that Acuña et al. (2021) predict, since their error-bars are calculated for CMF values that span between 0.2 and 0.3, and we only show the interior structure for CMF = 0.25.

Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Estimated water mass fraction, WMF, together with their uncertainties for TRAPPIST-1f, g and h with our model (green), Agol et al. (2021, purple), and Acuña et al. (2021, red).

6.2 Water content in the TRAPPIST-1 system and other possible interior structures

The extended Pre-Main-Sequence (PMS) phase of M-dwarf stars, like TRAPPIST-1, can cause atmospheric mass loss via thermal escape, and water loss due to the increased luminosity and XUV activity (see e.g., Wordsworth et al. 2018; Luger & Barnes 2015). This water mass loss depends on the planet’s orbital distance (Bourrier et al. 2017), making close-in planets more prone to volatile escape (Bolmont et al. 2017). Gialluca et al. (2024) calculated the water loss rates for the TRAPPIST-1 planets, and studied their implications on recent observations (Zieba et al. 2023; Greene et al. 2023), which imply a tenuous oxygen atmosphere or no atmosphere for planets b and c, respectively. According to Gialluca et al. (2024), however, the loss of atmospheres on the inner planets does not preclude the retention of substantial water inventories on the outer planets, where lower stellar irradiation could allow water to persist. Recent studies have assessed the complete loss of a water reservoir at the Galilean satellites, namely Io (Bennacer et al. 2026) and Europa, both which experience atmospheric sputtering (Johnson 1990), or plasma-driven escape of sodium and potassium (Gebek & Oza 2020). The sputter escape of water has yet to be modeled in detail at Europa and Ganymede (Oza et al. 2019b; Leblanc et al. 2023), in contrast to Io’s observed sodium escape by atmospheric sputtering (Burger et al. 2001), however is to first order, photoionization limited.

As seen from Figure 8, tidal heating rates and, thus Nvolc is highly dependent on the CMF that is assumed for the planets. We explore the range of probable CMFs for the planets that are compatible with their formation mechanism. The most favored scenario for the TRAPPIST-1 planets is formation in a protoplanetary disk, possibly beyond the water ice line, with a subsequent inward migration, which is supported by the resonant chain and compactness of the system (Teyssandier et al. 2022). Planet formation within a protoplanetary disk results in planets with similar ratios of relative refractory elements, meaning similar CMFs, for fully differentiated bodies (Bond et al. 2010; Elser et al. 2012), but with different volatile element content (Öberg & Bergin 2016). Thus, the WMF of the planets is also a function of their formation mechanism. If the TRAPPIST-1 planets are best described by the CMF of the most dense planet (planet c with 22–31% CMF), their densities can be fitted with different light-element content or water, according to Dorn et al. (2018). With our model, the latter CMF ratios would mean WMF values of ≈4 – 8 × 10−2, 3 – 6 × 10−2, and 1 – 11 × 10−2 for TRAPPIST-1f, g and h, respectively.

Varying the WMF in order to obtain the observed densities is compatible with the results of Agol et al. (2021, Figure 9), who found that all TRAPPIST-1 planets fall within one compositional line that is less dense than solar system terrestrial planets. Other possible compositions for the TRAPPIST-1 planets include either an interior with less iron compared to the terrestrial content, or oxidized iron in the mantle (Agol et al. 2021).

6.3 Future work

A first question to address in the future would be how to confirm that an observation of water comes from cryovolcanism. Exogenic sources of water are expected to remain constant over time (Paganini et al. 2020), whereas endogenic cryovolcanism may vary depending on the balance of interior processes such as tidal heating and radiogenic heat generation. Repeated observations could potentially contribute to the confirmation of the origin of water. The relationship between cryovolcanism and the different interior heat sources remains an open question, with significant implications for interpreting observations of cryovolcanic activity.

A second point of discussion would be a possible connection of observations of cryovolcanic activity to interior constraints. If tidal heating is the dominant driver of cryovolcanism, planets with low CMF values are more likely to exhibit detectable cryovolcanic signatures due to enhanced tidal dissipation. Conversely, if the total internal heat budget is dominated by radiogenic sources, planets with higher CMF values might exhibit stronger cryovolcanic activity.

Further investigation is needed to determine feasible interior models for a particular water column density observation. A more detailed analysis would need to consider uncertainties in interior parameters and ice rheology, the temperature dependence of rheological properties and the identification of secondary equilibrium states in mantle dynamics, as in Kleisioti et al. (2023). While our study outlines the potential sources of cryovolcanism and their implications, the detailed uncertainties inherent to interior modeling and observational constraints fall outside the scope of this work. Future research is needed to refine these models and better constrain the interplay between planetary structure, heat generation, and cryovolcanic activity.

7 Conclusions

In this work, we combine thermal and tidal heating models, previously applied to icy satellites in the Solar System on TRAPPIST-1f, g and h. We identify interior configurations compatible with the presence of subsurface oceans and show that, for all three planets, tidal heating within the high-pressure ice layer and radiogenic heating dominate over tidal dissipation in the silicate mantle and ice I shell under the assumptions of our model. Our dissipation partitioning differs from Solar System icy satellites, where dissipation is often concentrated in the ice I shell, because our TRAPPIST-1 solutions can contain very thick viscoelastic HPI layers which represent a large volume of deformable solid and therefore contribute substantially to the total tidal dissipation budget.

In addition, we constrain the ice I and subsurface ocean layer thicknesses for all three planets. The interior structures of the inner two planets, TRAPPIST-1f and g, would be in thermal equilibrium only in cases where the outer ice I shell is relatively thin, between ≈4 and 10 km, under the assumptions of our model. We also calculate the corresponding WMFs and CMFs that are compatible with such interior structures. We find that for an Earth-like CMF (0.32) TRAPPIST-1f, g and h, would possess a subsurface ocean in thermal equilibrium for a WMF of 8.5 × 10−2, 5.9 × 10−2 and 11.8 × 10−1.

Moreover, we investigate the detectability of cryovolcanic water vapor via transmission spectroscopy. We show that hydrostatic water atmospheres provide an upper bound on detectability, whereas non-hydrostatic, globally uniform sputtered exospheres require energetically unphysical efficiencies to produce detectable transmission signatures with JWST when pure water ice sublimation rates are low. Spatial localization of outgassed material in plume-like configurations substantially enhances effective line-of-sight column densities and reduces the energetic requirements for a detection. Such configurations are detectable on TRAPPIST-1f for volcanic mass fluxes between Earth-like and Io-like values with ≈20 transits, although such detections depend sensitively on the assumed energy conversion efficiencies and spatial localization of outgassed material. When high water ice sublimation rates are assumed for TRAPPIST-1f non-hydrostatic, globally uniform and localized sputtered clouds produce signals above the JWST/NIRISS detectability threshold corresponding to ≈20 transits, under the assumptions of our models. Since the high sublimation rate is standard for pure water ice (Johnson et al. 1981), our simulations suggest a potential detection (or non-detection) of water vapor is able to inform the pure water mass fraction of ice on TRAPPIST-1f’s surface.

Taken together, our results place constraints on the observability of cryovolcanic water vapor on TRAPPIST-1f, g and h with current facilities. While JWST observations need favorable conditions to detect cryovolcanic activity, future observations with next-generation space- (e.g., Habitable World Observatory) and ground-based telescopes (e.g., ELT/ANDES), combined with improved modeling of interior-surface interactions, will be essential for assessing the prevalence and observational signatures of cryovolcanism on icy exoplanets (Quick et al. 2023).

Acknowledgements

The authors would like to thank the anonymous reviewers for their constructive reviews which greatly strengthened the manuscript. The Python programming language (Rossum 1995) and the open-source Python packages numpy (van der Walt et al. 2011), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2013), petitRADTRANS (Mollière et al. 2019, 2020; Alei et al. 2022) and PandExo (Batalha et al. 2017) were used. For the calculation of tidal heating the LOV3D software package (Rovira-Navarro et al. 2024) was used. This research has been supported by the PEPSci Program (Planetary and ExoPlanetary Science Program), NWO, the Netherlands.

References

  1. Acuña, L., Deleuil, M., Mousis, O., et al. 2021, A&A, 647, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Agol, E., Dorn, C., Grimm, S. L., et al. 2021, PSJ, 2, 1 [NASA ADS] [Google Scholar]
  3. Alei, E., Konrad, B. S., Angerhausen, D., et al. 2022, A&A, 665, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Barr, A. C., Dobos, V., & Kiss, L. L. 2018, A&A, 613, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Batalha, N. E., Mandell, A., Pontoppidan, K., et al. 2017, PASP, 129, 064501 [Google Scholar]
  7. Baum, W., Kreidl, T., Westphal, J., et al. 1981, Icarus, 47, 84 [Google Scholar]
  8. Bennacer, Y., Mousis, O., & Hue, V. 2026, ApJ, 997, 70 [Google Scholar]
  9. Bierson, C. J., & Steinbrügge, G. 2021, PSJ, 2, 89 [Google Scholar]
  10. Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bolmont, E., Sastre, M., Revol, A., Kervazo, M., & Tobie, G. 2026, A&A, 707, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bond, J. C., O’Brien, D. P., & Lauretta, D. S. 2010, ApJ, 715, 1050 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121 [NASA ADS] [CrossRef] [Google Scholar]
  14. Burgasser, A. J., & Mamajek, E. E. 2017, ApJ, 845, 110 [Google Scholar]
  15. Burger, M. H., Schneider, N. M., de Pater, I., et al. 2001, ApJ, 563, 1063 [Google Scholar]
  16. Cadieux, C., Doyon, R., MacDonald, R. J., et al. 2024, ApJ, 970, L2 [NASA ADS] [CrossRef] [Google Scholar]
  17. Cameron, A. 1985, Icarus, 64, 285 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cassidy, T. A., & Johnson, R. E. 2010, Icarus, 209, 696 [Google Scholar]
  19. Croft, S. K., Kargel, J. S., Kirk, R. L., et al. 1995, in Neptune and Triton, eds. D. P. Cruikshank, M. S. Matthews, & A. M. Schumann, 879 [Google Scholar]
  20. de Wit, J., Wakeford, H. R., Lewis, N. K., et al. 2018, Nat. Astron., 2, 214 [NASA ADS] [CrossRef] [Google Scholar]
  21. Delrez, L., Gillon, M., Triaud, A. H. M. J., et al. 2018, MNRAS, 475, 3577 [NASA ADS] [CrossRef] [Google Scholar]
  22. Dobos, V., & Turner, E. L. 2015, ApJ, 804, 41 [NASA ADS] [CrossRef] [Google Scholar]
  23. Dobos, V., Barr, A. C., & Kiss, L. L. 2019, A&A, 624, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Dorn, C., Mosegaard, K., Grimm, S. L., & Alibert, Y. 2018, ApJ, 865, 20 [NASA ADS] [CrossRef] [Google Scholar]
  25. Dougherty, M. K., Khurana, K. K., Neubauer, F. M., et al. 2006, Science, 311, 1406 [NASA ADS] [CrossRef] [Google Scholar]
  26. Ehrenreich, D., des Etangs, A. L., Beaulieu, J.-P., & Grasset, O. 2006, ApJ, 651, 535 [Google Scholar]
  27. Elser, S., Meyer, M., & Moore, B. 2012, Icarus, 221, 859 [Google Scholar]
  28. Fagents, S. A. 2003, J. Geophys. Res. Planets, 108, 5139 [Google Scholar]
  29. Fagents, S. A., Greeley, R., Sullivan, R. J., et al. 2000, Icarus, 144, 54 [Google Scholar]
  30. Fagents, S. A., Lopes, R. M., Quick, L. C., & Gregg, T. K. 2022, in Comparative Planetology, 1, Planetary Volcanism across the Solar System, eds. T. K. Gregg, R. M. Lopes, & S. A. Fagents (Elsevier), 161 [Google Scholar]
  31. Fischer, H.-J., & Spohn, T. 1990, Icarus, 83, 39 [Google Scholar]
  32. Gebek, A., & Oza, A. V. 2020, MNRAS, 497, 5271 [Google Scholar]
  33. Geissler, P. E. 2015, in The Encyclopedia of Volcanoes, 2nd edn., eds. H. Sigurdsson, B. Houghton, S. R. McNutt, H. Rymer, & J. Stix (Amsterdam: Academic Press), 763 [Google Scholar]
  34. Gialluca, M. T., Barnes, R., Meadows, V. S., et al. 2024, PSJ, 5, 137 [Google Scholar]
  35. Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [Google Scholar]
  36. Grasset, O., Sotin, C., & Deschamps, F. 2000, Planet. Space Sci., 48, 617 [NASA ADS] [CrossRef] [Google Scholar]
  37. Green, A. P., Montesi, L. G. J., & Cooper, C. M. 2021, J. Geophys. Res. Planets, 126, e2020JE006677 [Google Scholar]
  38. Greene, T. P., Bell, T. J., Ducrot, E., et al. 2023, Nature, 618, 39 [NASA ADS] [CrossRef] [Google Scholar]
  39. Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Guenther, E. W., & Kislyakova, K. G. 2019, MNRAS, 491, 3974 [Google Scholar]
  41. Hand, K. P., Phillips, C. B., Murray, A., et al. 2022, PSJ, 3, 22 [Google Scholar]
  42. Hay, H. C. F. C., & Matsuyama, I. 2019, ApJ, 875, 22 [Google Scholar]
  43. Hedman, M., Gosmeyer, C., Nicholson, P., et al. 2013, Nature, 500, 182 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hemingway, D. J., & Mittal, T. 2019, Icarus, 332, 111 [Google Scholar]
  45. Henning, W. G., & Hurford, T. 2014, ApJ, 789, 30 [Google Scholar]
  46. Hernandez, J.-A., Caracas, R., & Labrosse, S. 2022, Nat. Commun., 13, 3303 [Google Scholar]
  47. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hussmann, H., & Spohn, T. 2004, Icarus, 171, 391 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hussmann, H., Spohn, T., & Wieczerkowski, K. 2002, Icarus, 156, 143 [NASA ADS] [CrossRef] [Google Scholar]
  50. Iess, L., Stevenson, D. J., Parisi, M., et al. 2014, Science, 344, 78 [NASA ADS] [CrossRef] [Google Scholar]
  51. Ingersoll, A. P., Ewald, S. P., & Trumbo, S. K. 2020, Icarus, 344, 113345 [NASA ADS] [CrossRef] [Google Scholar]
  52. Jaccard, C. 1976, J. Glaciol., 17, 155 [Google Scholar]
  53. Jia, X., Kivelson, M. G., Khurana, K. K., & Provan, G. T. 2018, Nat. Astron., 2, 459 [NASA ADS] [CrossRef] [Google Scholar]
  54. Johnson, R. E. 1990, Energetic Charged-Particle Interactions with Atmospheres and Surfaces (Berlin, Heidelberg: Springer-Verlag) [Google Scholar]
  55. Johnson, R. E., Lanzerotti, L. J., Brown, W. L., & Armstrong, T. P. 1981, Science, 212, 1027 [Google Scholar]
  56. Johnson, R. E., Oza, A., Young, L. A., Volkov, A. N., & Schmidt, C. 2015, ApJ, 809, 43 [NASA ADS] [CrossRef] [Google Scholar]
  57. Kleisioti, Dirkx, D., Rovira-Navarro, M., & Kenworthy, M. A. 2023, A&A, 675, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Kleisioti, E., Dirkx, D., Tan, X., & Kenworthy, M. A. 2024, A&A, 687, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Krissansen-Totton, J., & Fortney, J. J. 2022, ApJ, 933, 115 [NASA ADS] [CrossRef] [Google Scholar]
  60. Lau, H. C. P., Mitrovica, J. X., Austermann, J., et al. 2016, J. Geophys. Res. Solid Earth, 121, 6991 [NASA ADS] [CrossRef] [Google Scholar]
  61. Leblanc, F., Oza, A. V., Leclercq, L., et al. 2017, Icarus, 293, 185 [NASA ADS] [CrossRef] [Google Scholar]
  62. Leblanc, F., Roth, L., Chaufray, J., et al. 2023, Icarus, 399, 115557 [NASA ADS] [CrossRef] [Google Scholar]
  63. Leone, G., Bieger, K., & Soto, M. 2021, PSJ, 2, 29 [Google Scholar]
  64. Lobo, A. H., Shields, A. L., Palubski, I. Z., & Wolf, E. 2023, ApJ, 945, 161 [NASA ADS] [CrossRef] [Google Scholar]
  65. Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119 [Google Scholar]
  66. Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nat. Astron., 1, 0129 [Google Scholar]
  67. Manga, M., & Wang, C.-Y. 2007, Geophys. Res. Lett., 34, L07202 [NASA ADS] [CrossRef] [Google Scholar]
  68. Meier, T. G., Bower, D. J., Lichtenberg, T., Tackley, P. J., & Demory, B.-O. 2021, ApJ, 908, L48 [NASA ADS] [CrossRef] [Google Scholar]
  69. Meyer zu Westram, M., Oza, A. V., & Galli, A. 2024, J. Geophys. Res. Planets, 129, e2023JE007935 [Google Scholar]
  70. Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67 [Google Scholar]
  71. Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131 [Google Scholar]
  72. Neveu, M., & Rhoden, A. R. 2019, Nat. Astron., 3, 543 [Google Scholar]
  73. Nimmo, F., & Pappalardo, R. T. 2016, J. Geophys. Res. Planets, 121, 1378 [Google Scholar]
  74. Öberg, K. I., & Bergin, E. A. 2016, ApJ, 831, L19 [Google Scholar]
  75. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 [Google Scholar]
  76. Oza, A. 2017, PhD thesis, Sorbonne Université [Google Scholar]
  77. Oza, A. V., Johnson, R. E., Lellouch, E., et al. 2019a, ApJ, 885, 168 [Google Scholar]
  78. Oza, A. V., Leblanc, F., Johnson, R. E., et al. 2019b, Planet. Space Sci., 167, 23 [Google Scholar]
  79. Paganini, L., Villanueva, G. L., Roth, L., et al. 2020, Nat. Astron., 4, 266 [Google Scholar]
  80. Peltier, W. R. 1974, Rev. Geophys., 12, 649 [NASA ADS] [CrossRef] [Google Scholar]
  81. Phillips, C. B., McEwen, A. S., Hoppa, G. V., et al. 2000, Geophys. Res. Planets, 105, 22579 [Google Scholar]
  82. Pichierri, G., Morbidelli, A., Batygin, K., et al. 2024, Nat. Astron., 8, 1408 [NASA ADS] [CrossRef] [Google Scholar]
  83. Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597 [NASA ADS] [CrossRef] [Google Scholar]
  84. Porco, C., DiNino, D., & Nimmo, F. 2014, AJ, 148, 45 [NASA ADS] [CrossRef] [Google Scholar]
  85. Porco, C. C., Helfenstein, P., Thomas, P. C., et al. 2006, Science, 311, 1393 [Google Scholar]
  86. Postberg, F., Schmidt, J., Hillier, J., Kempf, S., & Srama, R. 2011, Nature, 474, 620 [NASA ADS] [CrossRef] [Google Scholar]
  87. Quick, L. C., & Marsh, B. D. 2015, Icarus, 253, 16 [Google Scholar]
  88. Quick, L. C., Roberge, A., Mlinar, A. B., & Hedman, M. M. 2020, PASP, 132, 084402 [NASA ADS] [CrossRef] [Google Scholar]
  89. Quick, L. C., Roberge, A., Mendoza, G. T., Quintana, E. V., & Youngblood, A. A. 2023, ApJ, 956, 29 [Google Scholar]
  90. Renaud, J. P., & Henning, W. G. 2018, ApJ, 857, 98 [Google Scholar]
  91. Roberts, J. H., & Nimmo, F. 2008, Icarus, 194, 675 [Google Scholar]
  92. Rossum, G. 1995, Python Reference Manual, Tech. rep., Centrum voor Wiskunde en Informatica (CWI) [Google Scholar]
  93. Roth, L. 2021, Geophys. Res. Lett., 48, e2021GL094289 [Google Scholar]
  94. Roth, L., Saur, J., Retherford, K. D., et al. 2014, Science, 343, 171 [NASA ADS] [CrossRef] [Google Scholar]
  95. Roth, L., Retherford, K. D., Saur, J., et al. 2026, A&A, 709, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Rovira-Navarro, M., van der Wal, W., Steinke, T., & Dirkx, D. 2021, PSJ, 2, 119 [Google Scholar]
  97. Rovira-Navarro, M., Katz, R. F., Liao, Y., van der Wal, W., & Nimmo, F. 2022, J. Geophys. Res. Planets, 127, e2021JE007117 [Google Scholar]
  98. Rovira-Navarro, M., Matsuyama, I., & Berne, A. 2024, PSJ, 5, 129 [Google Scholar]
  99. Rudolph, M. L., Manga, M., Walker, M., & Rhoden, A. R. 2022, in LPI Contributions, 2678, 53rd Lunar and Planetary Science Conference, 2562 [Google Scholar]
  100. Sabadini R., Vermeersen B., & Cambiotti G. 2016, Global Dynamics of the Earth: Applications of Viscoelastic Relaxation Theory to Solid-Earth and Planetary Geophysics, 2nd edn. (Berlin: Springer) [Google Scholar]
  101. Schoonenberg, D., Liu, B., Ormel, C. W., & Dorn, C. 2019, A&A, 627, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Schubert, G., Cassen, P., & Young, R. 1979, Icarus, 38, 192 [NASA ADS] [CrossRef] [Google Scholar]
  103. Schubert, G., Spohn, T., & Reynolds, R. T. 1986, in IAU Colloq. 77: Some Back-ground about Satellites, eds. J. A. Burns & M. S. Matthews (University of Arizona Press), 224 [Google Scholar]
  104. Segatz, M., Spohn, T., Ross, M. N., & Schubert, G. 1988, Icarus, 75, 187 [NASA ADS] [CrossRef] [Google Scholar]
  105. Shematovich, V., Johnson, R., Cooper, J., & Wong, M. 2005, Icarus, 173, 480 [Google Scholar]
  106. Shields, A. L., Meadows, V. S., Bitz, C. M., et al. 2013, Astrobiology, 13, 715 [NASA ADS] [CrossRef] [Google Scholar]
  107. Shoji, D., & Kurita, K. 2014, ApJ, 789, 3 [NASA ADS] [CrossRef] [Google Scholar]
  108. Smyth, W. H., & Marconi, M. L. 2006, Icarus, 181, 510 [Google Scholar]
  109. Sohl, F., Hussmann, H., Schwentker, B., Spohn, T., & Lorenz, R. D. 2003, J. Geophys. Res. Planets, 108 [Google Scholar]
  110. Soucek, O., Běhounková, M., Lanzendörfer, M., et al. 2024, Nat. Commun., 15, 7405 [Google Scholar]
  111. Sparks, W. B., Schmidt, B. E., McGrath, M. A., et al. 2017, ApJ, 839, L18 [Google Scholar]
  112. Spencer, J. R., Nimmo, F., Ingersoll, A. P., et al. 2018, in Enceladus and the Icy Moons of Saturn, eds. P. M. Schenk, R. N. Clark, C. J. A. Howett, A. J. Verbiscer, & J. H. Waite (University of Arizona Press), 163 [Google Scholar]
  113. Spohn, T., & Schubert, G. 2003, Icarus, 161, 456 [NASA ADS] [CrossRef] [Google Scholar]
  114. Teyssandier, J., Libert, A.-S., & Agol, E. 2022, A&A, 658, A170 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Tobie, G., Grasset, O., Dumoulin, C., & Mocquet, A. 2019, A&A, 630, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  117. van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
  118. van Woerkom, Q. B., & Kleisioti, E. 2024, A&A, 684, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Villanueva, G. L., Hammel, H. B., Milam, S. N., et al. 2023, Nat. Astron., 7, 1056 [NASA ADS] [CrossRef] [Google Scholar]
  120. Wordsworth, R. D., Schaefer, L. K., & Fischer, R. A. 2018, AJ, 155, 195 [NASA ADS] [CrossRef] [Google Scholar]
  121. Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proc. Natl. Acad. Sci. U.S.A., 116, 9723 [Google Scholar]
  122. Zieba, S., Kreidberg, L., Ducrot, E., et al. 2023, Nature, 620, 746 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A Melting temperature as a function of pressure

Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Melting temperature as a function of pressure for an ammonia containing ocean of 5 wt.% (Grasset et al. 2000; Sohl et al. 2003)

Appendix B Tidal heat generated in the Ice I outer

shell

Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Heat source distribution in the interior of TRAPPIST-1f. Qtidal,HPI Qtidal,m, and Qtidal,iceI is the tidal heat generated in the HPI layer, mantle and ice I shell, respectively. Qrad is the radiogenic heat produced in the mantle. The tidal heat that is generated in the ice I layer at least 4 orders of magnitude smaller than the total heat generated in the interior (Qeq).

Appendix C Interior structures as a function of Dice for the planets TRAPPIST-1g, h

Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

a. Interior structure as a function of Dice for the planet TRAPPIST-1g.
Thumbnail: Fig. C.2 Refer to the following caption and surrounding text. Fig. C.2

(a) Interior structure as a function of Dice for the planet TRAPPIST-1h.

All Tables

Table 1

Average physical and rheological properties of the interior model layers.

In the text
Table 2

Planet parameters used in the model.

In the text
Table 3

Interior and thermal parameters of TRAPPIST-1f in thermal equilibrium for different values of CMF.

In the text
Table 4

Interior and thermal parameters of TRAPPIST-1g in thermal equilibrium for different values of CMF.

In the text
Table 5

Interior and thermal parameters of TRAPPIST-1h in thermal equilibrium for different values of CMF.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Schematic overview of the interior model of TRAPPIST-1f, g and h. The ice I layer is divided into a conductive and a convective sublayer. The thermal model determines the thickness of the conductive and convective layers, Rice,I, RHPI, and Qeq. The tidal model computes the tidal heating rates in the HPI, and the mantle. The layer boundaries are not to scale. Note that the convective layer is only present for a sufficiently high (≳103) Rayleigh number.
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Schematic of the parameter calculation process. The red segment visualizes the heat transfer model in the outer ice shell, which is described in Section 2.1. The gray segment describes the process of calculating the thermal equilibrium in the entire planet (Section 2.3). The latter thermal equilibrium dictates that Qeq = Qtidal,HPI + Qrad.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Comparison of the heat contribution from several sources for TRAPPIST-1f as a function of ice shell thickness. Tidal heat in the mantle is negligible in the entire parameter space, while radiogenic heat and tidal heat in the ice shell both contribute significantly to the ice sheet equilibrium heat. TRAPPIST-1f due to its size and eccentricity has the highest tidal heat of the three planets.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Fractions of QtidalQeqMathematical equation: ${{{Q_{{\rm{tidal}}}}} \over {{Q_{{\rm{eq}}}}}}$ of the mantle, high-pressure ice layer, and outer ice I shell, as a function of Dice. The lighter and darker shaded regions indicate the 75% and 95% uncertainty intervals derived from the Monte Carlo distribution, respectively. Note that the y-axis scales differ between panels.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Top: fice as a function of different layers’ viscosity. Bottom: fmantle as a function of different layers’ viscosity. Enhanced tidal dissipation occurs in the ice I shell for stiffer HPI layers. Tidal dissipation in the mantle contributes more to the heat budget for stiffer HPI layers and lower mantle viscosity.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Interior structure as a function of Dice for the planet TRAPPIST-1f. (a) Radii of different layers as a function of Dice. (b) Thickness of the stagnant lid, convective and ocean layers as a function of Dice.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Synthetic transmission spectra of water atmospheres and exo-spheres on TRAPPIST-1f, taking into account cryovolcanic outgassing, sublimation, and volatile loss from XUV irradiation. The black data points are generated with PandExo for JWST/NIRISS. (a) Spectrum of a tenuous hydrostatic water atmosphere for Ntotal = 1.9 × 1019 (hydrostatic benchmark) produced using petitRADTRANS (purple line). Mock JWST data points assume 450 transits and a binned resolution of 12. (b) Spectrum of a uniformly distributed sputtered cloud produced using PROMETHEUS (blue line). The darker blue data point are binned observations, demonstrating how the signal at ≈1.9 microns gets diluted. This sub-figure corresponds to a Ntotal = 1.1 × 1024 molecules m−2 (low sublimation case). Mock JWST data points assume 20 transits. (c) Same as (b) but for Ntotal = 1.6 × 1026 molecules m−2 which corresponds to the benchmark of the high sublimation case. Mock JWST data points assume 15 transits.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Ntotal as a function of Ṁvolc for TRAPPIST-1f. The shaded red region indicates parameter space requiring efficiencies η > 1, which are unphysical. The hydrostatic and sputtered cloud benchmarks represent the column densities observable with JWST/NIRISS with 450 transits for a hydrostatic atmosphere and a collisionless exosphere, respectively. The darkest gray region indicates an exosphere’s column density observable with 20 transits. Top: detectability of water column densities for a globally uniform sputtered cloud. Bottom: detectability of water column densities assuming the sputtered material is confined to a surface fraction f = 0.004, yielding enhanced line-of-sight columns, which represent plume-like configurations.
In the text
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Estimated water mass fraction, WMF, together with their uncertainties for TRAPPIST-1f, g and h with our model (green), Agol et al. (2021, purple), and Acuña et al. (2021, red).
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

Melting temperature as a function of pressure for an ammonia containing ocean of 5 wt.% (Grasset et al. 2000; Sohl et al. 2003)
In the text
Thumbnail: Fig. B.1 Refer to the following caption and surrounding text. Fig. B.1

Heat source distribution in the interior of TRAPPIST-1f. Qtidal,HPI Qtidal,m, and Qtidal,iceI is the tidal heat generated in the HPI layer, mantle and ice I shell, respectively. Qrad is the radiogenic heat produced in the mantle. The tidal heat that is generated in the ice I layer at least 4 orders of magnitude smaller than the total heat generated in the interior (Qeq).
In the text
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

a. Interior structure as a function of Dice for the planet TRAPPIST-1g.
In the text
Thumbnail: Fig. C.2 Refer to the following caption and surrounding text. Fig. C.2

(a) Interior structure as a function of Dice for the planet TRAPPIST-1h.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.