Characterisation of the hydrospheres of TRAPPIST-1 planets

Context. Planetary mass and radius data are showing a wide variety in densities of low-mass exoplanets. This includes sub-Neptunes, whose low densities can be explained with the presence of a volatile-rich layer. Water is one of the most abundant volatiles, which can be in the form of different phases depending on the planetary surface conditions. To constrain their composition and interior structure, it is required to develop models that calculate accurately the properties of water at its different phases. Aims. We present an interior structure model that includes a multiphase water layer with steam, supercritical and condensed phases. We derive the constraints for planetary compositional parameters and their uncertainties, focusing on the multiplanetary system TRAPPIST-1, which presents both warm and temperate planets. Methods. We use a 1D steam atmosphere in radiative-convective equilibrium with an interior whose water layer is in supercritical phase self-consistently. For temperate surface conditions, we implement liquid and ice Ih to ice VII phases in the hydrosphere. We adopt a MCMC inversion scheme to derive the probability distributions of core and water compositional parameters Results. We refine the composition of all planets and derive atmospheric parameters for planets b and c. The latter would be in a post-runaway greenhouse state and could be extended enough to be probed by space mission such as JWST. Planets d to h present condensed ice phases, with maximum water mass fractions below 20%. Conclusions. The derived amounts of water for TRAPPIST-1 planets show a general increase with semi-major axis, with the exception of planet d. This deviation from the trend could be due to formation mechanisms, such as migration and an enrichment of water in the region where planet d formed, or an extended CO$_{2}$-rich atmosphere.


Introduction
Ongoing space missions such as CHEOPS (Benz 2017) and TESS (Ricker et al. 2015), and their follow-up with groundbased radial velocity telescopes, are confirming the existence of low-mass exoplanets with a wide range of densities. These densities range from the values typically inferred for the Earth or Mercury to those measured in Uranus and Neptune. The exoplanets in the former class are mainly composed of a Fe-rich core and a silicate mantle, while the latter class has a layer that is rich in volatiles. Water is the most abundant and least dense volatile after H and He (Forget & Leconte 2014), which makes it a likely species to constitute the volatile reservoir in these planets. Several studies have investigated the interior structure and composition of water-rich planets (Sotin et al. 2007;Seager et al. 2007;Dorn et al. 2015;Zeng et al. 2019), but focused mainly on its condensed phases. Nonetheless, many sub-Neptunes are close to their host star and receive enough irradiation to trigger a runaway greenhouse state in which water is present as steam. In some cases, the high pressure and temperature conditions can render the hydrosphere in supercritical and plasma, or even superionic phases (Mazevet et al. 2019;French et al. 2016). Therefore, it is crucial to include the modeling of all possible phases of water to provide an accurate description of its presence on the planetary surface. Moreover, the surface conditions are determined by the greenhouse effect caused by atmospheric gases, making the modelling of radiative-convective equilibrium in atmospheres a key parameter to determine in which phase water could be present on the surface. Most of interior structure models represent the planetary atmosphere as a gas layer with a simplified isothermal temperature profile (Dorn et al. 2018(Dorn et al. , 2017b, which is very different from the temperature profile in the convective deep layers of thick atmospheres (Marcq 2012).
Multiplanetary systems are unique environments that present both planets that can hold condensed phases as well as highlyirradiated planets with steam atmospheres. In this study, we develop a planet interior model suitable for the different conditions at which water can be found in low-mass planets. Our implementation includes a supercritical water layer, introduced in Mousis et al. (2020), coupled with a 1D radiative-convective atmosphere model (Marcq 2012;Marcq et al. 2017;Pluriel et al. 2019) to calculate the total radius of the highly-irradiated planets with water self-consistently. Furthermore, for temperate planets, we have updated the interior model presented in Brugger et al. (2016Brugger et al. ( , 2017 to include ice phases Ih, II, III, V and VI. We introduce these models in a MCMC Bayesian analysis scheme adapted from Dorn et al. (2015). This allows us to derive the water mass fraction (WMF) and core mass fraction (CMF) with their associated confidence intervals that reproduce the observed radius, mass and stellar composition measurements.
We use this model to explore the possible water content of the TRAPPIST-1 system, an ultra-cool M dwarf that hosts seven low-mass planets in close-in orbits. Three of these planets are located in the habitable zone , meaning that they can hold liquid water or ice Ih on their surfaces. Although all planets in TRAPPIST-1 system have masses and radii that are characteristic of rocky planets, their differences in density indicate that each planet has a different volatile content. This makes this planetary system ideal for testing planet interior, atmospheric structure and formation scenarios.
In Sect. 2 we describe the complete interior structure model, including the new updates for the supercritical and ice phases, the coupling between the interior and the atmosphere for steam and supercritical planets, and the MCMC Bayesian algorithm. The parameters for the TRAPPIST-1 planets used in this study are summarized in Sect. 3, including mass, radius and Fe/Si molar ratio. The results of our analysis of the hydrospheres of TRAPPIST-1 planets are described in Sect. 4. We compare our results with previous works and discuss the implications of our water estimates for planet formation in Sect. 5. We finally expose our conclusions in Sect. 6.

Planetary structure model
For consistency, we recall the main principles of the interior structure model. The basis of our model is explained in Brugger et al. (2016Brugger et al. ( , 2017. The 1D interior structure model takes as input the mass and the composition of the planet, which is parameterized by the CMF and WMF. The structure of the planet is stratified in 3 layers: a core, a mantle and a hydrosphere. The pressure, temperature, gravity acceleration and density are computed at each point of the one-dimensional spatial grid along the radius of the planet. The pressure, P(r), is obtained by integrating the hydrostatic equilibrium (Eq. 1); the gravitational acceleration, g(r), by solving Gauss's theorem (Eq. 2); the temperature, T (r), with the adiabatic gradient (Eq. 3); and the density, ρ(r), with the Equation of State (EOS). m is the mass at radius r, G is the gravitational constant, and γ and φ are the Gruneisen and the seismic parameters, respectively. Their formal macroscopic definitions are shown in equation 4, where E is the internal energy and V is the volume. The Gruneisen parameter is a thermodynamic parameter that describes the dependence of the vibrational properties of a crystal with the size of its lattice. It relates the temperature in a crystalline structure to the density, which is calculated by the EOS. The seismic parameter defines how seismic waves propagate inside a material. It is related to the slope of the EOS at constant temperature (Brugger et al. 2017;Sotin et al. 2007).
The boundary conditions are the temperature and pressure at the surface, and the gravitational acceleration at the center of the planet. The value of the latter is zero. The total mass of the planet is calculated with Eq. 5, which is derived from the conservation of mass (Brugger et al. 2017;Sotin et al. 2007). Once the total input mass of the planet is reached and the boundary conditions are fulfilled, the model has converged.
Depending on the surface conditions, the hydrosphere can be present in supercritical, liquid or ice states. For each of these phases of water, we use a different EOS and Gruneisen parameter to compute their P-T profiles and density accurately. In Sect. 2.1 we describe the updates to the supercritical water layer with respect to the model depicted in Mousis et al. (2020), while in Sect. 2.2 we present the implementation of the hydrosphere in ice phases. Finally, the coupling between the atmosphere and the interior model with planets whose hydrosphere is in steam or supercritical phases is explained in Sect. 2.3, followed by the description of the MCMC algorithm in Sect. 2.4.

Supercritical water
If the planet is close enough to its host star, the upper layer of the hydrosphere corresponds to a hot steam atmosphere, whose temperature at the base is determined by the radiativeconvective balance calculated by the atmosphere model (Marcq 2012;Marcq et al. 2017). When the pressure and temperature at the surface, which is defined as the base of the hydrosphere layer, are above the critical point of water, we include a supercritical water layer extending from the base of the hydrosphere to a height corresponding to the phase change to steam (Mousis et al. 2020). We updated the EOS for this layer to the EOS introduced by Mazevet et al. (2019), which is a fit to the experimental data provided by the International Association for the Properties of Water and Steam (IAPWS) (Wagner & Pruß 2002) for the supercritical regime, and quantum molecular dynamics (QMD) simulations data for plasma and superionic water (French et al. 2009). The IAPWS experimental data span a temperature range of 251.2 to 1273 K and of 611.7 to 10 9 Pa in pressure, while their EOS can be extrapolated up to 5000 K in temperature and 10 11 Pa in pressure (Wagner & Pruß 2002). The validity range of the EOS presented in Mazevet et al. (2019) includes that of the IAPWS plus the region in which the QMD simulations are applicable, which corresponds to a temperature from 1000 K to 10 5 K and densities in the 1-10 2 g/cm 3 range. These densities are reached at high pressures, i.e., in the 10 9 -10 12 Pa range. Following Eq. 3, the adiabatic gradient of the temperature is specified by the Gruneisen and the seismic parameters. These are dependent on the derivatives of the pressure with respect to the density and the internal energy (Eq. 4). We make use of the specific internal energy and density provided by Mazevet et al. (2019) to calculate them.

Ice phases
We extended the hydrosphere in Brugger et al. (2016Brugger et al. ( , 2017 with liquid and high pressure ice VII by adding 5 more condensed phases: ice Ih, II, III, V and VI. EOS for ice Ih has been developed by Feistel & Wagner (2006) with minimization of the Gibbs potential function from the fit of experimental data. It covers all the pressure and temperature range in which water forms ice Ih. Fei et al. (1993) proposed a formalism to derive the EOS of ices II, III and V. These EOS have the form of V = V(P, T ), Table 1: EOS and reference thermal parameters for ices Ih, II, III, V and VI. This includes the reference values for the density ρ 0 , the temperature T 0 , the bulk modulus K T 0 and its derivative K T 0 , the heat capacity C p (T 0 ), and the thermal expansion coefficient α 0 . References.
(1) Gagnon et al. (1990); (2) Báez & Clancy (1995); ( Feistel & Wagner (2006) which can be found by integrating the following differential equation (Tchijov et al. 2004): where α is the thermal expansion coefficient and β the isothermal compressibility coefficient. If the relationship between the specific volume, V, and the pressure, P, at a constant temperature T = T 0 is determined, Eq. 6 can be integrated as: Fei et al. (1993) proposed the following expression for the thermal expansion coefficient α: where η is an adjustable parameter estimated from the fitting of experimental data. Its value is 1.0 for ice II and ice III (Leon et al. 2002) and 7.86 for ice V (Shaw 1986). ρ is the density, α(P 0 , T ) is the coefficient of thermal expansion at a reference pressure P 0 , K T 0 is the isothermal bulk modulus at the reference temperature T 0 , and K T 0 is the first derivative of the isothermal bulk modulus at the reference temperature. Hence, by substituting Eq. 8 in Eq. 7 and integrating, we obtain the following EOS for high-pressure ice: The final expression (Eq. 9) requires the knowledge of the variation of the specific volume, V(P, T 0 ), with pressure at the reference temperature T 0 . Moreover, the variation of the density with temperature, ρ(T ), and the bulk modulus with its derivative at the reference temperature, K T 0 and K T 0 , must also be provided.
In Table 1 we specify the data and references to obtain these parameters for each ice phase.
In the case of ice VI, we adopt the second-order Birch-Murnaghan (BM2) formulation, which is: where ρ 0 is the reference density for ice VI. We also introduce a thermal correction to the density since the pressure also depends on the temperature: where α 0 is the reference coefficient of thermal expansion. Interfaces between liquid and ice layers are established by phase transition functions from Dunaeva et al. (2010).

Interior-atmosphere coupling
We use a one-dimensional atmosphere model designed to compute radiative transfer and pressure-temperature (P, T ) profiles for water and CO 2 atmospheres (Marcq 2012;Marcq et al. 2017). The formation of water clouds is considered in the computation of the albedo. The atmosphere is in radiative equilibrium, and presents a composition of 99% water and 1% CO 2 . The density of steam is obtained using a non-ideal EOS (Haar et al. 1984).
If the surface pressure is below 300 bar, the atmosphere and the interior are coupled at the atmosphere-mantle boundary and water does not reach the supercritical regime. However, if the surface pressure is greater than 300 bar, the atmosphere and the interior are coupled at this pressure level and a layer of water in supercritical phase forms between the atmosphere and the mantle. The pressure level at 300 bar is close enough to the critical point of water at 220 bar to avoid the atmosphere model take over pressures and temperatures where the temperature profile is adiabatic.
The top-of-atmosphere pressure is set to 20 mbar, which corresponds to the observable transiting radius (Mousis et al. 2020;Grimm et al. 2018). We denote the radius and mass from the center of the planet to this pressure level the total radius and mass, R total and M total , respectively. We also define the radius and the mass that comprise the core, mantle and supercritical layers as the bulk radius and mass, R bulk and M bulk , respectively. The atmosphere model provides the Outgoing Longwave Radiation (OLR), albedo, thickness and mass of the atmosphere as a function of the bulk mass and radius, and the surface temperature. If the atmosphere of the planet is in radiative equilibrium, the OLR is equal to the radiation the planet absorbs from its host star, F abs . The OLR depends on the effective temperature since OLR Fig. 1: Structural diagram of the coupling between the interior structure model and the atmosphere model. T base is the temperature at the bottom of the steam atmosphere in radiativeconvective equilibrium. z and M atm denote the atmospheric thickness and mass, respectively. R bulk and M bulk correspond to the planet bulk radius and mass, respectively. R guess refers to the initial guess of the bulk radius, while R interior is the output bulk radius of the interior structure model in each iteration.
= σT 4 eff , where σ corresponds to the Stefan-Boltzmann constant. To calculate the absorbed radiation F abs , we first compute the equilibrium temperature, which is where A B is the planetary albedo, R and T are the radius and effective temperature of the host star, respectively. a d is the semimajor axis of the planet. The absorbed radiation is then calculated as Figure 1 shows the algorithm we implemented to couple the planetary interior and the atmosphere. The interior structure model calculates the radius from the center of the planet to the base of the steam atmosphere. For a fixed set of bulk mass and radius, the OLR depends on the surface temperature. Consequently, the surface temperature at which the OLR equals the absorbed radiation corresponds to the surface temperature that yields radiative equilibrium in the atmosphere. This is estimated with a root-finding method. Since the bulk radius is an output of the interior model (R interior ) and an input of the atmosphere model, we first need to calculate the surface temperature for a certain mass and composition with an initial guess bulk radius. Then this surface temperature is the input for the interior model, which provides the bulk radius. With this bulk radius, we can generate a new value of the surface temperature. This scheme is repeated until the bulk radius converges to a constant value, to which we add the thickness of the atmosphere, z, to get the total radius of the planet R total . The total mass M total is obtained as the sum of the bulk mass M bulk plus the atmospheric mass M atm . The tolerance used to determine if the bulk radius has achieved convergence is 2% of the bulk radius in the previous iteration. This is approximately 0.02 R ⊕ for an Earth-sized planet.

MCMC Bayesian analysis
We adapted the MCMC Bayesian analysis algorithm described in Dorn et al. (2015) to our coupled interior and atmosphere model. The input model parameters are the bulk planetary mass M bulk , the CMF and the WMF. Therefore, m = {M bulk , CMF, W MF}, following the notation in Dorn et al. (2015). Depending on the planetary system and their available data, we can have observational measurements of the total planetary mass and radius and the stellar composition, or only the total planetary mass and radius. The available data in the former case is denoted as d = {M obs , R obs , Fe/S i obs }, while the data in the latter case is represented as d = {M obs , R obs }. The uncertainties on the measurements are σ(M obs ), σ(R obs ), and σ(Fe/S i obs ).
The CMF and WMF prior distributions are uniform distributions between 0 and a maximum limit. This maximum limit is 75% for the CMF, which is derived from the maximum estimated Fe/Si ratio of the proto-Sun (Lodders et al. 2009). With this limit on the CMF, we are assuming that the exoplanets have not been exposed to events during or after their formation that could have stripped away all of their mantle, such as mantle evaporation or giant impacts. In addition, the maximum WMF is set to 80%, which is the average water proportion found in comets in the solar system (McKay et al. 2019). The prior distribution for the mass is a Gaussian distribution whose mean and standard deviations correspond to the central value and uncertainties of the observations. The MCMC scheme first starts by randomly drawing a value for each of the model parameters from its respective prior distributions. This set of values is designated as m 1 = M bulk,1 , CMF 1 , W MF 1 . The index i = 1 corresponds to the first proposed set of input values within the first chain, n = 1. The model calculates the total mass and radius and the theoretical Fe/Si mole ratio, which are the set of output parameters g(m 1 ) = {R 1 , M 1 , Fe/S i 1 }. The likelihood of a set of model parameters is then calculated via the following relationship (Dorn et al. 2015): where the normalization constant of the likelihood function C is defined as: When the Fe/Si mole ratio is not available as data, the square residual term of the Fe/Si mole ratio is removed from Eq. 14, as well as its squared uncertainty in Eq. 15.
Subsequently we draw a new set of input parameters, m 2 = M bulk,2 , CMF 2 , W MF 2 from the prior distributions within the same chain, n. We assure that the absolute difference between the values for i = 1 and i = 2 is lower than a fixed step, which is the maximum size of the perturbation. This guarantees that the new state m 2 is uniformly bounded and centered around the old state, m 1 . The maximum perturbation size is selected so that the acceptance rate of the MCMC, which is defined as the ratio between the number of models that are accepted and the number of proposed models, is above 20%. After m 2 is chosen, the forward model calculates its corresponding output parameters and obtains their likelihood L(m 2 | d), as shown in Eq. 14. The acceptance probability is estimated with the log-likelihoods l(m | d) = log(L(m | d)) as: If P accept is greater than a number drawn from a uniform distribution between 0 and 1, m 2 is accepted and the chain moves to the state characterised by m 2 , starting the next chain n + 1. Otherwise, the chain remains in the state of m 1 and a different set of model parameters is proposed as m 3 . To make sure that the posterior distributions converge and that all parameter space is explored, we run 10 4 chains. In other words, with acceptance rates between 0.2 and 0.6, the MCMC proposes between 1.6 and 5 × 10 4 sets of model inputs. Agol et al. (2020) have performed an analysis of TTVs that includes all transit data from Spitzer since the discovery of the system. We adopt these data for the mass, radius and semi-major axis in our interior structure analysis (Table 2).

System parameters of TRAPPIST-1
TRAPPIST-1 does not have available data regarding its chemical composition. However, the Fe/Si abundance ratio can be estimated assuming that TRAPPIST-1 presents a similar chemical composition to that of other stars of the same metallicity, age and stellar population. As proposed by Unterborn et al. (2018), we select a sample of stars from the Hypatia Catalog (Hinkel et al. 2014(Hinkel et al. , 2016(Hinkel et al. , 2017. We choose the set of stars by constraining the C/O mole ratio to be less than 0.8, and the stellar metallicity between -0.04 and 0.12, as this is the metallicity range calculated for TRAPPIST-1 by Gillon et al. (2017). We discard thick disk stars since TRAPPIST-1 is likely a thin disk star. Our best-fit Gaussian to the distribution of the Fe/Si mole ratio shows a mean of 0.76 and a standard deviation of 0.12. Since this Fe/Si value is an estimate based on the chemical composition of a sample of stars that belong to the same stellar population of TRAPPIST-1, we present two scenarios for each planet. In scenario 1, the only available data are the planetary mass and radius, while scenario 2 includes the estimated stellar Fe/Si mole ratio to constrain the bulk composition.
For temperate planets that cannot have a steam atmosphere, we set the surface temperature in our interior model to their equilibrium temperatures assuming an albedo zero (Table 2). Although surface temperatures for thin atmospheres are lower than that obtained with this assumption, the dependence of the bulk radius on surface temperature for planets with condensed water is low. For example, if we assume a pure water planet of 1 M ⊕ with a surface pressure of 1 bar, the increase in radius due to a change of surface temperature from 100 K to 360 K is 0.002 R ⊕ , which is less than 0.2% of the total radius, i.e., 10 times less than our convergence criterion. Additionally, the atmospheres of TRAPPIST-1 planets in the habitable zone and farther are significantly thinner than those of the highly-irradiated planets. Lincowski et al. (2018) estimate thicknesses of approximately 80 km for temperate planets in TRAPPIST-1, which is negligible compared to their total radius. Therefore, we only calculate the atmospheric parameters (OLR, surface temperature, albedo and thickness of the atmosphere) for planets that present their hydrospheres in steam phase. Tables 3 and 4 show the retrieved parameters, including the total planetary mass and radius, and the Fe/Si mole ratio. In both scenarios, we retrieve the mass and radius within the 1σ-confidence interval of the measurements for all planets. In scenario 1, where only the mass and radius data are considered, we retrieve Fe/Si mole ratios without any assumptions on the chemical composition of the host star. Although the uncertainties of these estimates are more than 50% in some cases, we can estimate a common Fe/Si mole ratio for the planetary system. This common Fe/Si range is determined by the overlap of the 1σ confidence intervals of all planets, which corresponds to Fe/Si = 0.45-0.97. This interval is compatible with the Fe/Si mole ratio of 0.76±0.12 proposed by Unterborn et al. (2018). This overlap can also be seen in Fig. 2, which presents the 1σ-confidence regions derived from the 2D marginalized posterior distributions of the CMF and WMF. The minimum value of the common CMF is determined by the lower limit of the confidence region of planet g, which is approximately 0.23, whereas the common maximum CMF value corresponds to the upper limit of planets b and c, which is 0.4. This is partially in agreement with the CMF obtained in scenario 2, when we assume the Fe/Si mole ratio proposed by Unterborn et al. (2018), which are found between 0.2 and 0.3 (Table 4). Thus, the CMF of the TRAPPIST-1 planets could be compatible with an Earth-like CMF (CMF ⊕ = 0.32).   Table 4: Output parameters retrieved by the MCMC method for all TRAPPIST-1 planets: total mass (M ret ) and radius (R ret ), CMF, WMF and Fe/Si molar ratio. In this case the Fe/Si mole ratio estimated by following Unterborn et al. (2018) is also included as data (scenario 2).

CMF and WMF posterior distributions
In scenario 1, the retrieved WMF for all planets in the system are below 20% within their uncertainties. This maximum WMF limit reduces to 10% for scenario 2. This indicates that the TRAPPIST-1 system is poor in water and other volatiles, especially the inner planets b and c. Both planets are compatible with a dry composition in both scenarios, although the presence of an atmosphere cannot be ruled out given the possible CMF range estimated in scenario 1. Figure 3 shows the OLR calculated by the atmosphere model and the absorbed radiation (Eq. 12 and 13) for planets b, c and d. For temperatures lower than ∼ T sur f = 2000 K, the OLR has little dependency on the surface temperature. This is caused by the nearly constant temperature (between 250 and 300 K) of the radiating layers in the thermal IR range (Goldblatt et al. 2013) and it is related to the runaway greenhouse effect (Ingersoll 1969). We obtain a constant OLR or an OLR limit (Nakajima et al. 1992) of 274.3, 273.7 and 254.0 W/m 2 for planets b, c and d, respectively. These are close to the OLR limit obtained by Katyal et al. (2019) of 279.6 W/m 2 for an Earth-like planet. The small difference is due to their different surface gravities. As explained in Sect. 2.3, if the atmosphere can find a surface temperature at which the OLR and the absorbed radiation are equal, their atmospheres are in global radiative balance. This is the case for planets b and c, whose surface temperatures are approximately 2450 K and 2250 K, respectively. These are above the temperatures where the blanketting effect is effective, named T ε in Marcq et al. (2017) implying that the atmospheres of planets b and c are in a post-runaway state. However, planet d is not in global radiative balance since its absorbed radiation never exceeds its OLR. This means that planet d would be cooling down, and an internal flux of approximately 33 W/m 2 would be required to supply the extra heat to balance its radiative budget. TRAPPIST-1 inner planets are likely to present an internal heat source due to tidal heating (Barr et al. 2018;Dobos et al. 2019;Turbet et al. 2018). The tidal heat flux estimated for planet d is F tidal = 0.16 W/m 2 ( Barr et al. 2018), which is one order of magnitude lower than needed for radiative-convective balance of a steam atmosphere. Due to the blanketting effect of radiation over the surface of planet d, the OLR limit is larger than the absorbed radiation and hence the planet can cool enough to present its hydrosphere in condensed phases. and 4.4 × 10 −6 , which correspond to a surface pressure of 128.9 bar and 4.85 bar, respectively. The thermal structure of their steam atmospheres are dominated by a lower, unsaturated troposphere where water condensation does not occur. Then the atmosphere consists of a middle, saturated troposphere where cloud formation would be possible, extending up to 10 mbar, and finally an isothermal mesosphere above. Since we consider a clear transit radius of 20 mbar Mousis et al. 2020) the presence of clouds above this pressure level would flatten the water features in the planetary spectrum Katyal et al. 2020). On the other hand, planets d and e could present water in liquid phase, which could be partially or completely covered in ice Ih. While the hydrosphere of planet h is not massive enough to attain the high pressures required for ice VII at its base, planets d, to g can reach pressures up to a 100 GPa.  of TRAPPIST-1 b and c for a water-dominated atmosphere in scenario 1. The total thickness of an atmosphere is related to its scale height, which is defined as H = RT/µg, where R = 8.31 J/K mol is the gas constant, T is the mean atmospheric temperature, µ the mean molecular mass and g the surface gravity acceleration. For planets b and c, the mean atmospheric temperatures are 940.4 and 499.4 K, and their surface gravities are 10.8 and 10.7 m/s 2 , respectively. The mean molecular mass for a 99% water and 1% CO 2 atmosphere is 18.3 g/mol. The mean temperature increases with surface temperature, while the mean molecular mass is determined by the composition of the atmosphere. For the same composition and surface gravity, the scale height and therefore the thickness of the atmosphere is directly correlated to the surface temperature. As shown in Fig. 5, the atmospheric thickness, z atm increases with the surface temperature T sur f . This is known as the runaway greenhouse radius inflation effect (Goldblatt 2015;Turbet et al. 2019), where a highly irradiated atmosphere is more extended than a colder one despite having similar compositions. For planet b, its atmosphere can extend up to 450 km, while planet c presents a maximum extension of 300 km. The minimum limit for the thicknesses is zero, which corresponds to the case of a dry composition. Ortenzi et al. (2020) estimated that for a planet of 1-1.5 M ⊕ the maximum atmospheric thickness due to the outgassing of an oxidised mantle is 200 km, which is compatible with the ranges we have obtained for the atmospheric thicknesses. Scenario 2 shows the same trends for the atmospheric parameters but with lower atmospheric mass and surface pressure. With their WMF posterior distributions centered in zero and low standard deviation, the surface pressure is below 1 bar and atmospheric thicknesses below 100 km in most of the accepted models, which means that in scenario 2 planets b and c are most likely dry rocky planets. Agol et al. (2020) use the interior and atmosphere models presented in Dorn et al. (2018) and Turbet et al. (2020b) to obtain the WMF estimated of the TRAPPIST-1 planets with updated and more precise radii and masses data from Spitzer TTVs (Agol et al. 2020). We thus limited the comparison to the sole results of Agol et al. (2020) with the same input values. By doing so, we can be certain that the variations in WMF estimates are due to our different modelling approach. Figure 6 shows that planets b and c are most likely dry in scenario 2, where the resulting CMF are between 0.2 and 0.3 for the whole system. We  Barr et al. (2018), where only mass and radius data were taken into account. The lower panel corresponds to scenario 2, whose CMF is constrained in a narrow range between 0.2 and 0.3, while for Agol et al. (2020) we show the WMF for a CMF of 0.25.

WMF comparison with previous works
obtain maximum estimates of 3.4 × 10 −6 and 2.7 × 10 −6 for b and c, respectively. For the same density, the estimated value of the WMF depends on the CMF that is considered. Therefore we compare WMF estimates for similar CMF between this work and Agol et al. (2020). We show our WMF in scenario 2, since the CMF of all planets spans a narrow range between 0.2 and 0.3, which are the most similar values to one of the CMF assumed by Agol et al. (2020), CMF = 0.25. Our WMF for the steam planets of the TRAPPIST-1 system are in agreement with Agol et al. (2020), who calculated a maximum WMF of 10 −5 for a constant CMF of 0.25. We are able to reduce the maximum limit of the water content of the highly irradiated planets compared to previous studies and establish the most likely WMF with our coupled atmosphere-interior model. The calculation of the total radius requires a precise determination of the atmospheric thickness. This depends strongly on the surface temperature and the surface gravity, which are obtained with radiative transfer in the atmosphere, and the calculation of the gravity profile for a bulk mass and composition in the interior self-consistently.
In the case of planet d, we estimate a WMF of 0.036 ± 0.028, while Agol et al. (2020) obtain an upper limit of 10 −5 . The latter estimate considers that water is in vapor form, which is less dense than condensed phases, while our model shows that the surface conditions allow liquid or ice phases, resulting in a higher WMF. This discrepancy in the possible water phases on the surface of planet d is due to different atmospheric compositions. We consider a water-dominated atmosphere with 1% CO 2 , while Agol et al. (2020) and Turbet et al. (2020b) assume a N 2 and H 2 O mixture. This difference in composition changes radiative balance since CO 2 is a strong absorber in the IR compared to N 2 , which is a neutral gas. Nonetheless, N 2 is subject to stellar wind-driven escape and it is unlikely to be stable for the inner planets of TRAPPIST-1, while CO 2 is more likely to survive thermal and ion escape processes (Turbet et al. 2020a).
Our WMF for planets e to h are in agreement within uncertainties with Agol et al. (2020), although their central values are significantly lower. The EOS employed to compute the density of the water layers in Agol et al. (2020) is also used in Dorn et al. (2018) and Vazan et al. (2013), which agrees well with the widely-used SESAME and ANEOS EOSs (Baraffe et al. 2008). These EOS are not consistent with experimental and theoretical data since they overestimate the density at pressures higher than 70 GPa (Mazevet et al. 2019). This yields an underestimation of the WMF for the same total planetary density and CMF.
For the specific case of scenario 1, with no assumptions on the stellar composition and the Fe/Si mole ratio, we compared our CMF and WMF with those obtained in Barr et al. (2018) (Figure 6 and Table 5). These authors use masses and radii data given by Wang et al. (2017). They obtain lower masses compared to Agol et al. (2020) while their radii are approximately similar, which would explain why Barr et al. (2018) tend to overestimate the water content of the TRAPPIST-1 planets. Moreover, most of the mass uncertainties in Wang et al. (2017) are 30-50%, while the mass uncertainties obtained by Agol et al. (2020) are 3-5%. This causes Barr et al. (2018) to calculate wider CMF and WMF 1σ confidence intervals. In addition, there are differences between our interior modelling approach and that of Barr et al. (2018). For example, according to their results, planet b can have up to 50% of its mass as water. This high WMF value is due to the assumption that the hydrosphere is in liquid and ice I phases, and high-pressure ice polymorphs (HPPs), which are more dense than the steam atmosphere we consider. In contrast, the CMF seems to be closer to our estimates, especially for planet b, d and e, where their maximum CMF is approximately 0.40, in agreement with our CMF 1σ intervals.
We can also discuss the possible habitability of the hydrospheres of the TRAPPIST-1 planets by comparing our WMF estimates with the layer structure as a function of planetary mass Planet CMF Barr et al. (2018) Table 5: Comparison between our one-dimensional 1σ confidence regions for the CMF and those of Barr et al. (2018). We show only estimates for scenario 1, since Barr et al. (2018) did not consider any constraints on the Fe/Si ratio based on stellar composition.  and water content obtained by Noack et al. (2016) . According to Noack et al. (2016), a habitable hydrosphere must be structured in a single liquid water ocean or in several ice layers that enable the formation of a lower ocean layer. This lower ocean would be formed by the heat supplied by the mantle that melts the high pressure ice in the ice-mantle boundary (Noack et al. 2016). For planet d, a surface liquid ocean would form for all its possible WMF if the atmosphere allows for the presence of condensed phases. For planets e, f and g, the hydrosphere could be stratified in a surface layer of ice Ih and a liquid or an ice II-VI layer. In the case we had low-pressure ices II-VI, their base could be melted by the heat provided by the mantle, and form a lower ocean layer as suggested by Noack et al. (2016). At WMF ≥ 0.10, less than 50% of the possible configurations enable a habitable sub-surface ocean layer, and at a WMF ≥ 0.14, the hydrosphere is uninhabitable. In scenario 1, planets e to g reach these values within uncertainties, although their minimum values extend down to 0-0.03 in WMF, which would be the habitable regime.

System formation and architecture
In the case of scenario 1, where no Fe/Si data is assumed, the WMF increases with the distance to the star with the exception of planet h, whose WMF is similar to that of planet d. In the case of scenario 2, where a common Fe/Si of 0.76 ± 0.12 is assumed for the whole system, the WMF increases with the distance to the star (Fig. 6) with the exception of planet d whose WMF is similar to that of planet f, which is more water-rich than planet e. This slight deviation from the observed trend could be explained by migration, where planet d could have formed beyond the snow line and then migrated inwards (Raymond et al. 2018). In addition, pebble ablation and water recycling back into the disk could have been less efficient in the case of planet d, compared to planet e (Coleman et al. 2019). On the other hand, the gas at the distance at which planet d formed could have been more enriched in volatiles than the outer planets, accreting more water ice than planet e in a 'cold finger' (Stevenson & Lunine 1988;Cyr et al. 1998). Pebble formation in the vicinity of the water iceline can induce important enhancements of the water ice fraction in those pebbles due to the backward diffusion of vapor through the snowline and the inward drift of ice particles. Therefore, if a planet forms from this material, it should be more water-rich than those formed further (Mousis et al. 2019). These formation scenarios could explain the high WMF of planet d when we assume that its water layer is in condensed phases. Post-formation processes could also have shaped the trend of the WMF with axis, such as atmospheric escape due to XUV and X-ray emission from their host star. Bolmont et al. (2017) estimated a maximum water loss of 15 Earth Oceans (EO) for TRAPPIST-1 b and c and 1 EO for planet d. If we assume that the current WMF are the central values of the posterior distributions we derived in scenario 1, planets b, c and d would have had an initial WMF of 2.37 × 10 −3 , 2.50 ×10 −3 and 0.085, respectively. Therefore, atmospheric escape would have decreased the individual WMF of each planet, but the increase of WMF with distance from the star would have been preserved.
In addition to the WMF-axis trend, we can differentiate the very water-poor, close-in planets, b and c, from the outer, waterrich planets, d to h. This has been reported as a consequence of pebble accretion in the formation of other systems, such as the Galilean moons. While Io is dry, Callisto and Ganymede are water-rich, with Europa showing an intermediate WMF of 8% (Ronnet et al. 2017). Pebble-driven formation can produce planets with WMF ≥ 15% if these are formed at the water ice line (Coleman et al. 2019;Schoonenberg et al. 2019). In contrast, planets formed within the ice line would present WMF less than 5% (Liu et al. 2020;Coleman et al. 2019), which is close to the mean value we calculated for planet h, 5.5%. The maximum WMF limit in the first scenario is approximately 20%. This maximum limit is significantly lower than the typical WMF generated by planetesimal accretion scenario, which is 50-40% (Miguel et al. 2020). Therefore, our results are consistent with the pebble-driven formation scenario.

Alternative atmospheric compositions
However, the atmosphere of planet d could be dominated by other atmospheric gases different from H 2 O-based mixtures, which could produce an extended atmosphere and increase the total planetary radius. Hydrogen-dominated atmospheres have been deemed unlikely as one of the possible atmospheric compositions for all planets in the TRAPPIST-1 system, both cloudfree (de Wit et al. 2016(de Wit et al. , 2018 and with high-altitude clouds and hazes Ducrot et al. 2020). Similarly, CH 4dominated atmospheres are not probable given the photometry data of the Spitzer Space Telescope (Ducrot et al. 2020). Therefore, our best candidate to explain the low density of planet d in a water-poor scenario is CO 2 . We find that a CO 2 -dominated at-mosphere with 1% water vapor in planet d would be in radiativeconvective equilibrium by computing the OLR and absorbed radiation, as we did for water-dominated atmospheres. The resulting surface temperature is approximately 950 K, which is slightly higher than the surface temperature of Venus (700 K) with a higher water vapor mixing ratio. Figure 7 introduces the mass-radius relationships for different CMF, assuming a CO 2dominated atmosphere with a surface pressure of 300 bar. It shows that planet d is compatible with a planet with a CO 2dominated atmosphere and CMF between 0.2 and 0.3, which is a very likely CMF range for TRAPPIST-1 planets based on our analysis. Surface pressures lower than 300 bar would yield lower atmospheric thicknesses, so it would be necessary to consider a lower CMF to explain the observed density of planet d. CO 2 in the case of planet d can be provided by volcanic outgassing , since its internal heat flux produced by tidal heating is in the range 0.04-2 W/m 2 , which favours plate tectonics (Papaloizou et al. 2018). Secondary CO 2 -dominated atmospheres could have traces of O 2 , N 2 and water vapor.

Conclusions
We presented an interior structure model for low-mass planets at different irradiations that is valid for a wide range of water phases, derived from the approaches of Brugger et al. (2017) and Mousis et al. (2020). For highly-irradiated planets, we couple a 1D water steam atmosphere in radiative-convective equilibrium with a high-pressure convective layer in supercritical phase. The density in this layer is computed by using an accurate EOS for high-pressure and high-temperature water phases. For temperate planets whose surface conditions allow the formation of condensed phases, we implemented a hydrosphere with liquid water and ice phases Ih, II, III, V, VI and VII. In addition, we adapted the MCMC Bayesian algorithm described in Dorn et al. (2015) to our interior model to derive the posterior distributions of the compositional parameters, WMF and CMF, given mass, radius and stellar composition data. We then applied our interior model to the particular case of TRAPPIST-1 planets using the latest mass and radius data from Spitzer (Agol et al. 2020).
The hydrospheres of TRAPPIST-1 planets have been characterised by calculating their P-T profiles and thermodynamic phases. Planets b and c are warm enough to present steam atmospheres. They could hold post-runaway greenhouse atmospheres with thicknesses up to 450 km and surface temperatures up to 2500 K, which means that they are extended enough to be suitable targets for atmospheric characterisation by future spacebased facilities such as JWST. Moreover, planets d to g present their hydrospheres in condensed phases. These hydrospheres can contain high-pressure ices that start to form at 10 9 − 10 10 Pa.
We have obtained CMF and WMF probability distributions for all planets in the system. We found that the Fe/Si mole ratio of the system is in the 0.45-0.97 range without considering any assumption on the chemical composition of the stellar host. This Fe/Si range corresponds to a CMF value in the 0.23-0.40 range, making the CMF of TRAPPIST-1 planets compatible with an Earth-like value (0.32). In addition, our WMF estimates agree within uncertainties with those derived by Agol et al. (2020), although their most likely values are considerably lower for planets with condensed phases. In the case of planets with steam hydrospheres, their densities are compatible with dry rocky planets with no atmospheres. Nevertheless, we cannot rule out the presence of an atmosphere with the Fe/Si range we have derived without any assumption on the chemical composition of the host star. When considering a possible estimate of the Fe/Si ratio of the host star (scenario 2), we obtained lower maximum limits of the WMF for planets b and c compared to previously calculated limits by Agol et al. (2020) for a similar CMF of 0.25. Our estimated WMF in steam and condensed phases are consistent with an increase of WMF with progressing distance from the host star. This trend, as well as the maximum WMF we calculated, favour pebble-driven accretion as a plausible formation mechanism for the TRAPPIST-1 system. However, planet d presents a slightly higher WMF than planet e. This could be due to processes that took place during planet formation, such as migration, a lowefficient ablation of pebbles and gas recycling or an enhancement of the water ice fraction in pebbles at the distance of the disc where planet d formed. An extended atmosphere dominated by greenhouse gases different from a water-dominated atmosphere, such as CO 2 , could also explain the low-density of planet d compared to planet e.
Future work should include more atmospheric processes and species that determine the mass-radius relations of planets with secondary atmospheres in the Super-Earth and sub-Neptune regime. These can vary the atmospheric thickness and increase the total planetary radius with varying atmospheric masses while other compositional parameters change the bulk radius. These should be integrated in one single interior-atmosphere model, combined in a MCMC Bayesian framework such as the one we used in this study. This statistical approach has been employed with interior models for planets with H/He-dominated atmospheres (Dorn et al. 2017b(Dorn et al. ,a, 2018, or dry planets (Plotnykov & Valencia 2020), but not for planets with secondary, CO 2 and steam-dominated atmospheres. The integrated model should also include a description of escape processes, such as hydrodynamic or Jeans escapes, which is particularly interesting to explore the lifetime of secondary atmospheres. Close-in, low-mass planets are likely to outgas atmospheric species, such as CO 2 , and form O 2 via photodissociation of outgased H 2 O, during their magma ocean stage or due to plate tectonics (Chao et al. 2020). Thus, a mixture of these gases should be considered to study the thermal structure of planets with secondary atmospheres. Planets b and c in the TRAPPIST-1 system could present magma oceans due to their high surface temperatures (T ≥ 1300 K) (Barr et al. 2018;Chao et al. 2020), and the maximum surface pressure we have obtained here can be used to assess the current outgassing rate in magma ocean studies (Noack et al. (2017), Baumeister et al. submitted) and better constrain the WMF for the interior magma ocean models (e.g Katyal et al. (2020)) in the future.