The Implications of Thermal Hydrodynamic Atmospheric Escape on the TRAPPIST-1 Planets

JWST observations of the seven-planet TRAPPIST-1 system will provide an excellent opportunity to test outcomes of stellar-driven evolution of terrestrial planetary atmospheres, including atmospheric escape, ocean loss, and abiotic oxygen production. While most previous studies use a single luminosity evolution for the host star, we incorporate observational uncertainties in stellar mass, luminosity evolution, system age, and planetary parameters to statistically explore the plausible range of planetary atmospheric escape outcomes. We present probabilistic distributions of total water loss and oxygen production as a function of initial water content, for planets with initially pure water atmospheres and no interior–atmosphere exchange. We find that the interior planets are desiccated for initial water contents below 50 Earth oceans. For TRAPPIST-1e, f, g, and h, we report maximum water-loss ranges of 8.0−0.9+1.3 , 4.8−0.4+0.6 , 3.4−0.3+0.3 , and 0.8−0.1+0.2 Earth oceans, respectively, with corresponding maximum oxygen retention of 1290−75+75 , 800−40+40 , 560−25+30 , and 90−10+10 bars. We explore statistical constraints on initial water content imposed by current water content, which could inform evolutionary history and planet formation. If TRAPPIST-1b is airless while TRAPPIST-1c possesses a tenuous oxygen atmosphere, as initial JWST observations suggest, then our models predict an initial surface water content of 8.2 −1.0+1.5 Earth oceans for these worlds, leading to the outer planets retaining >1.5 Earth oceans after entering the habitable zone. Even if TRAPPIST-1c is airless, surface water on the outer planets would not be precluded.


INTRODUCTION
As M dwarfs are the most common stars in our local stellar neighborhood (Cantrell et al. 2013), whether their planetary systems can harbor life is a key question in astrobiology that may be amenable to observational tests in the near term.Due to their small size and low temperature, it is considerably easier to conduct observations of terrestrial planets in M dwarf systems with current and near-term technology than for any other main sequence host star (Tarter et al. 2007).Terrestrial planetary targets of interest for atmospheric characterization with M dwarf hosts may be accessible with the JWST (Gardner et al. 2006;Seager et al. 2009;Deming et al. 2009;Morley et al. 2017;Lustig-Yaeger et al. 2019;Wunderlich et al. 2019;Gialluca et al. 2021), and the ground-based extremely large telescopes (ELTs; e.g., E-ELT, TMT, GMT) (Gilmozzi & Spyromilio 2007;Snellen et al. 2013;Marconi et al. 2016;Ezaki et al. 2016;Fanson et al. 2018;López-Morales et al. 2019;Currie et al. 2023;Hardegree-Ullman et al. 2023).In the more distant future, the flagship Habitable Worlds Observatory (HWO), as defined by the recent Astronomy 2020 Decadal Survey, may include a small subset of M dwarf planets in its observing sample, that will be accessible via transmission and/or direct imaging (LUVOIR Team et al. 2019;Gaudi et al. 2020;National Academies of Sciences, Engineering, and Medicine et al. 2021).
Despite the observational advantages of M dwarf systems, they face a number of habitability challenges due to the proximity of the habitable zone (HZ) to the parent star.This leads to strong gravitational interactions, such as tidal locking, that can result in synchronous rotation and extreme day-night temperature gradients (Showman et al. 2013;Yang et al. 2014;Barnes 2017;Lobo et al. 2023).The M dwarf star's relative proximity to its HZ planets also enhances radiative interactions, and these could pose the most significant habitability and characterization challenges.The early and pronounced luminosity evolution of the star may produce a runaway greenhouse on the planets (Ingersoll 1969;Pollack 1971;Nakajima et al. 1992;Goldblatt et al. 2013; Barnes et al. 2013;Turbet et al. 2019), potentially promoting an extended magma ocean phase (Hamano et al. 2013;Schaefer et al. 2016;Barth et al. 2021) and severe atmospheric mass loss through thermal atmospheric escape (e.g., Hunten 1973).In addition to luminosity evolution, heightened stellar activity also increases the stellar XUV of M dwarf stars, which enhances atmospheric loss, and may complicate characterization and interpretation of spectra from these planets by driving planetary photochemistry that can enhance, destroy or even mimic proposed biosignatures (Segura et al. 2005;Meadows et al. 2018).
During the extended pre-main sequence (PMS) phase, an M dwarf is at its brightest as it contracts to its main sequence size (Baraffe et al. 2015).The protostar reaches radiative equilibrium but not nuclear fusion, leading to extended Kelvin-Helmholtz contractions to compensate for the loss of energy through radiation (Hayashi 1961;Lamers & Levesque 2017).These contractions lead to a decrease in radiative opacity, allowing energy to escape more rapidly while the star continues to accrete mass from the surrounding disk -both effects contribute to the star being orders of magnitude brighter than on the main sequence (e.g., Baraffe et al. 1998;Dotter et al. 2008;Reid & Hawley 2013;Lamers & Levesque 2017).Additionally, the PMS phase is considerably longer for less massive stars; while a Sun-like star may spend tens of millions of years in the PMS, an M dwarf may spend on the order of a billion years in the PMS (Baraffe et al. 2015).
This increased luminosity and XUV activity during the PMS phase of an M dwarf results in prolonged incident radiation that is well in excess of the runaway greenhouse limit on planets that will be in the HZ during the main sequence.This early PMS phase may drive thermal hydrodynamic atmospheric escape, which can strongly abrade a planetary atmosphere, and modify the redox state (Jeans 1925;Hunten 1973;Lammer et al. 2003;Volkov et al. 2011;Luger & Barnes 2015).In the hydrodynamic escape regime, incident radiation leads to water photolysis in the atmosphere and the hydrogen is subsequently lost, leading to permanent water loss and the potential for abiotic oxygen build-up (Wordsworth & Pierrehumbert 2014;Luger & Barnes 2015;Lincowski et al. 2018).Since oxygen production from ocean-loss could produce thousands of bars of O 2 (Luger & Barnes 2015;Lincowski et al. 2018), with at least a few bars retained after atmospheric, crustal, and magma ocean losses (Schaefer et al. 2016;Barth et al. 2021), hydrodynamic escape on planets orbiting M dwarfs could produce significant abiotic oxygen that is potentially detectable and can be confused for biologically-generated O 2 (e.g., Segura et al. 2007;Schwieterman et al. 2016a,b;Meadows 2017;Meadows et al. 2018).Consequently, it is highly important to understand the mechanisms that generate abiotic O 2 , and to quantify the likely range of oxygen production that may occur during this hydrodynamic escape process.
The 7-planet TRAPPIST-1 system (Gillon et al. 2016(Gillon et al. , 2017) ) provides an excellent, nearby (12 pc), and relatively observationally accessible example of a transiting M dwarf system with which we can explore and test the outcomes of enhanced atmospheric escape for terrestrial-sized planets.As a late-type M8V star, TRAPPIST-1 likely provided an early heightened XUV environment (Baraffe et al. 2015), and its 7 transiting terrestrial planets, including 3 -4 within the current habitable zone, enable tests of atmospheric loss processes as a function of distance from the parent star.As several of these planets may follow similar evolutionary pathways as Venus or Earth, there is also a prime opportunity for comparative planetology, particularly in questioning what makes a potentially habitable planet inhospitable (Kane et al. 2019).
TRAPPIST-1 is a high priority target for JWST General and Guaranteed Time Observations (Gillon et al. 2020;Greene et al. 2023;Zieba et al. 2023), and quantifying likely water loss ranges and abiotic oxygenation histories will be needed to interpret current and upcoming observations.While possible post-ocean and atmospheric loss atmospheres, and habitable environments in this system have already been modeled in several previous studies (Lincowski et al. 2018;Krissansen-Totton et al. 2018;Lustig-Yaeger et al. 2019;Barth et al. 2021), until recently (Birky et al. 2021; Krissansen-Totton & Fortney 2022) a lack of constraints on the XUV evolution have prohibited comprehensive statistical studies of PMS evolution outcomes.
Whereas previous studies (e.g., Luger & Barnes 2015;Tian 2015;Schaefer et al. 2016;Bolmont et al. 2017;Bourrier et al. 2017;Lincowski et al. 2018;Johnstone 2020) have explored implications for planetary evolution by analyzing planetary escape rates driven by single or tens of stellar luminosity evolution tracks, more comprehensive statistical studies are needed to constrain the broader family of evolutionary solutions that are consistent with preexisting data.For example, while a single stellar evolutionary track may provide a plausible history for water loss on a planet, it cannot be used to explore the phase space of water loss and oxygen production outcomes, quantify confidence in the single interpretation, or place upper limits on maximum water loss and abiotic O 2 production.By leveraging the computationally efficient planetary evolution model, VPLanet (Barnes et al. 2020), and developing support scripts (Vconverge) to automatically run simulation suites to convergence, here we explore the full observational uncertainty of the stellar evolution (Birky et al. 2021) and planetary parameters (Agol et al. 2021) of the TRAPPIST-1 system.To define the plausible phase space of water loss and oxygen production outcomes, we conduct hundreds of thousands of simulations of the water evolution history of the TRAPPIST-1 planets across a grid of initial water contents from 1 terrestrial ocean (TO) to 500 TO -spanning the full range that has been suggested as plausible in the literature (Schoonenberg et al. 2019;Agol et al. 2021;Acuña et al. 2021;Raymond et al. 2022).A terrestrial ocean, or Earth ocean, is our standard unit of measurement for water content and is equivalent to 1.39 × 10 21 kg.Though Krissansen-Totton & Fortney (2022) adopted a similar statistical approach considering the system constraints from Birky et al. (2021) and Agol et al. (2021), our study conducts a greater number of simulations to more comprehensively explore the phase space of initial water content, although at the expense of a coupled atmosphere and interior.We present predictions with 1σ uncertainty on water loss as a function of initial water abundance for each TRAPPIST-1 planet.This study additionally identifies the maximum water loss possible for the HZ planets regardless of increasing initial water mass, assuming initially pure water atmospheres and no exchange with the interior.Furthermore, given this maximum water loss and our understanding of oxygen drag during escape, we calculate the maximum amount of abiotic oxygen that could be produced through this process, with perfectly ineffective oxygen sinks at the surface, to provide upper limits on total atmospheric oxygen build-up.
Through this statistical study of all 7 known planets in the system, escape evolution as a function of orbital distance is explored, methods of inferring system evolution from observation are developed and demonstrated, model intercomparisons are explicitly quantified (Lincowski et al. 2018;Krissansen-Totton & Fortney 2022), and we begin to interpret the most recent JWST observations (Greene et al. 2023;Zieba et al. 2023) in the context of our findings.To infer system evolution, we show how our knowledge of water loss over time can be used in conjunction with present-day water mass fraction (WMF) measurements of the TRAPPIST-1 planets to retrieve plausible evolutionary histories and probable initial water contents.The retrieval of evolutionary history will aid in developing a relationship between initial and final water content that can connect formation models to observations.In regards to JWST, we consider recent observations (Greene et al. 2023;Zieba et al. 2023) in the context of our results by analyzing the planets individually, and then by showing the significant additional information that can be discerned when considering multiple bodies in the system and using a comparative planetology approach.In particular, we provide insights into initial water content of the system and current plausible conditions on planets that have yet to be observed.
Section 2 describes our modeling methods, including a description of the atmospheric escape mechanisms we consider ( §2.1), construction of model input ( §2.2), how convergence is achieved through multiple simulations ( §2.3), and a comparison to the model before these most recent modifications (as presented in Lincowski et al. 2018) ( §A.1).Section 3 presents the key results of our study, including the evolutionary pathways observed in our simulations ( §3.1), the full statistical spread of all models ( §3.2), identification of the most probable initial water content to reproduce final water mass fractions presented in Agol et al. (2021) ( §3.3.1),quantification of oxygen retention efficiency ( §3.4), and a comparison of our results to a similar study published in Krissansen-Totton & Fortney (2022) ( §3.5).The interpretation, significance, and limitations of these results and comparison to past work is discussed in Section 4, notably including an analysis of recent JWST observations in the context of our study ( §4.3.2).We summarize our conclusions in Section 5. Please note, throughout the paper, in referring to individual planets, we adopt the shorthand of 'T1-i' for 'TRAPPIST-1i', where i would refer to the planet in question (b, c, d, e, f, g, or h).

METHODS
We conduct our simulations using VPLanet, an open source code for modeling star, planet, and planetary system interactions and their impact on evolution (Barnes et al. 2020).VPLanet is broken up into modules for different physical processes; in this work, the STELLAR (stellar evolution, §2.1.1)and AtmEsc (atmospheric escape, §2.1.2) modules were used, and we specifically modify the treatment of diffusion-limited escape, the critical drag XUV flux, and the allowed behavior of the ratio of oxygen to hydrogen loss.A mathematically detailed description of our modifications and how they operate in AtmEsc can be found in Appendix A.
Input planetary and stellar parameters for our models come from Agol et al. (2021) and Birky et al. (2021), respectively.Using these updated constraints, we build and preserve self-consistent relationships between posteriors (see §2.2) allowing us to randomly sample a string of self-consistent planetary and stellar inputs for each simulation, as opposed to each input individually.Typically, 3500 -5500 simulations are conducted per initial water content per planet to build distributions of final water and oxygen content.The exact number of simulations in a particular suite is determined by model convergence (found via the Kolmogorov-Smirnov, or K-S, test), which is described in Section 2.3.In Appendix A.1, we illustrate the effects of our modifications to AtmEsc by comparing it to a simulation run with an older, validated version of the model (Luger & Barnes 2015;Lincowski et al. 2018Lincowski et al. , 2022)).
We note that a small number of outliers stemming from ∼1.24% of the stellar mass samples were identified and removed from our distributions.These outlying stellar mass samples were greater than 4 standard deviations from the average and were likely a product of the machine learning method implemented in Birky et al. (2021).A more detailed description of outlier removal is given in Appendix B.

Stellar Evolution (STELLAR)
The STELLAR module evolves the host star characteristics over time by performing a bicubic spline integration across mass and time of the evolutionary tracks in Baraffe et al. (2015) to track the stellar radius, radius of gyration, effective temperature, and luminosity from the PMS phase to present-day.It then uses the empirical broken power law formula of Ribas et al. (2005) to track the XUV output as a function of bolometric luminosity, which is particularly relevant for M dwarfs such as TRAPPIST-1: where L XU V is the XUV luminosity, L Bol is the bolometric luminosity, t is the time, t sat is the XUV saturation time, β XU V is the exponential decay rate of the XUV luminosity, and f sat is the XUV saturation fraction, which is the ratio of XUV to bolometric luminosity during the saturation time; the latter 3 of these parameters are user defined inputs (see §2.2).The STELLAR module was used to fit the TRAPPIST-1 stellar input parameters that we use in this study (Fleming et al. 2020;Birky et al. 2021), and a more thorough and mathematical description of STELLAR can be found in Barnes et al. (2020); STELLAR was not modified in this work.The output stellar parameters from the STELLAR module become the input for the atmospheric escape model, AtmEsc, at each timestep; the stellar XUV characteristics, evolved by STELLAR, are the main driver of hydrodynamic thermal atmospheric escape.

Atmospheric Escape (AtmEsc)
Atmospheric escape can typically be split up into the thermal and non-thermal regimes; due to the heightened XUV environment, the TRAPPIST-1 planets are most likely subjected to primarily thermal escape processes during the PMS phase.For the innermost planets that never enter the habitable zone, thermal escape dominates throughout the lifetime of the atmosphere, but for the remaining planets, once the star reduces in luminosity sufficiently for them to enter the habitable zone, non-thermal escape processes have the potential to be more significant (e.g., Dong et al. 2018).Within the thermal regime we can consider Jeans escape or hydrodynamic blow-off -the latter occurring when the internal energy of the escaping species approaches the kinetic energy required for escape (Hunten 1973); hydrodynamic loss is likely to provide a larger impact on the TRAPPIST-1 planetary atmospheres than Jeans escape, due to its extremely high mass loss rates.
Hydrodynamic escape can further be categorized as energy-or diffusion-limited.Energy-limited escape occurs when the escape flux is limited by the amount of energy incident on the planet, whereas diffusion-limited escape occurs when the escape flux is limited by the speed at which the escaping gas (typically hydrogen) can diffuse through a static background gas.During energy-limited escape, the escaping species may also drag a heavier species with it -the efficiency of which is highly dependent on mixing ratios in the atmosphere, the incoming XUV flux, and the mass of the species being dragged; thus, it can lead to observable effects such as mass fractionation of the heavier dragged species (e.g., Hunten et al. 1987;Zahnle et al. 1990), but this process produces only very weak fractionation of the hydrodynamically escaping species (i.e., D/H ratio, 18 O/ 16 O) (Mandt et al. 2009;Lammer et al. 2020).This drag can also be a considerable source of atmospheric mass loss, particularly during a star's most active period.In the diffusion-limited regime, the escape flux must be less than the maximum upward flux for which the background atmosphere is static (Hunten 1973); thus, by definition, oxygen drag cannot occur under diffusion-limited escape.In the diffusion-limited regime mass fractionation may occur (Hunten et al. 1989;Zahnle et al. 1990), enhancing the atmosphere's D/H ratio, but is likely only considerable on large timescales (Mandt et al. 2009).
Here we use VPLanet's AtmEsc module to model both energy-and diffusion-limited escape for the TRAPPIST-1 planets, and a full description of this module can be found in Barnes et al. (2020), with a more succinct mathematical description of its implementation in the context of this work provided in Appendix A. In comparison to previous work with AtmEsc (Luger & Barnes 2015;Lincowski et al. 2018Lincowski et al. , 2022;;Barnes et al. 2020), we modify the treatment of the switch to diffusion-limited escape, implement a calculation of the critical drag XUV flux (Eq.A5), and change the allowed behavior of η (ratio of oxygen escape to production, Eq.A8) to be less restrictive; these changes follow Schaefer et al. (2016).In this work, the switch from energy-to diffusion-limited escape occurs when the mixing ratio of free oxygen (i.e., oxygen not contained in water molecules) is greater than or equal to the mixing ratio of water.The critical drag XUV flux (F Crit XU V , Eq. A5) is defined as the minimum incoming XUV flux required for oxygen drag to occur; if the incoming XUV flux is less than this, oxygen drag will be halted regardless of whether or not the escape is energy-limited.Finally, η is defined as the ratio of oxygen escape to production, or the ratio of oxygen to hydrogen loss.In past AtmEsc implementations, this was restricted to be ≤1, meaning, at most, 1 oxygen atom could be lost for every 2 hydrogen atoms.Here, we allow η to vary up to ≤2 as we note that an accumulated oxygen reservoir may be depleted under highly favorable drag conditions, effectively meaning oxygen is lost faster than it is produced for certain periods of time.Appendix A.1 shows a comparison of the modified model used in this work to the previous version of the model (Luger & Barnes 2015;Lincowski et al. 2018Lincowski et al. , 2022)); in general, our modifications lead to slightly higher water loss (∼1 -8%) and thus slight higher oxygen build-up (∼10%).As simplifications are made within this implementation, Section 4.4.1 discusses future model improvements that could be made for a more realistic escape treatment.
Our simulations make several assumptions, including initially pure water atmospheres, no surface oxygen sinks or geological activity, and that all water is on the planetary surface and available for escape by 1 Myr.As in Wordsworth et al. (2018), we assume no primordial H 2 -He envelope to be conservative in finding upper limit values of water loss and oxygen production.We explore initial water contents from 1 terrestrial ocean (TO) to 500 TO on a modified log space grid.We choose the upper limit for initial water content simulated here to be 500 TO, based on several recent studies which constrain the initial and present-day water mass fractions and suggest these are around 10%wt (e.g., Schoonenberg et al. 2019;Agol et al. 2021;Raymond et al. 2022), which would be equivalent to ∼430 TO on a 1 Mearth planet; 500 TO is approximately bracketed between 10%wt of the larger planets, T1-b, c, and g, and the smaller planets, T1-d, e, f, and h.As we assume all oxygen not lost through escape remains in the atmosphere, our standard unit of measurement for oxygen production is bars of atmospheric pressure, which considers the mass and radius of a particular planet when converting from kg of oxygen.We additionally assume uniform mixing ratios for gases in the atmosphere (no vertical structure) with an isothermal temperature profile and do not consider photochemical processes besides the division of water molecules through photolysis.The effects our assumptions may have on the results is explored in Section 4.4.

Input Description
In determining simulation inputs, we consider not only observational uncertainty in the stellar luminosity evolution and physical planetary parameters (Birky et al. 2021;Agol et al. 2021), but also that planetary parameters such as mass, radius, orbit and insolation can be self-consistently dependent on the assumed stellar environment, as planetary mass and radius are derived from planet-to-star ratios.To account for this, we construct self-consistent chains of input parameters for the entire system, such that one simulation randomly samples an input chain as opposed to each individual parameter separately.
Here we detail the stellar and planetary parameters we consider, and how the input chain construction for a single simulation is accomplished, which is then repeated for simulation suites at a particular initial water content.The stellar inputs used in our work (from Birky et al. 2021) are the mass of the star (M * ), the XUV saturation fraction (f sat ), the XUV saturation time (t sat ), the exponential decay rate of the XUV luminosity (β XU V ), and the age (τ ).For each of the 7 planets, the input parameters used (from Agol et al. 2021) are the planetary mass ratio (M p /M * ), the planetary radius ratio (R p /R * ), the eccentricity (e), and the orbital period (P ).More specifically, Agol et al. (2021) presents a transit timing and photodynamical analysis of the system; the former, transit timing variations (TTV), are highly sensitive to the ratios of planetary to stellar mass, along with the eccentricity and orbital period.The photodynamical analysis provides sensitive constraints on the stellar density (ρ * ) and planet-to-star radius ratios by combining the mass ratios and orbital parameters derived from the TTV analysis to improve the derivation of planetary and stellar densities from Spitzer photometry (Agol et al. 2021).Thus, the actual planetary mass and radius of a single input chain will be dependent on the sampled stellar parameters.
To construct one input chain, one suite of stellar, TTV, and photodynamical inputs are randomly sampled.The selected stellar mass is then combined with the stellar-to-planetary mass ratios from the TTV sample to produce the planetary masses: where i refers to the planet in question (b, c, d, e, f, g, or h).Using Eq. ( 2) to find the planetary mass creates a relationship between the stellar inputs and the TTV inputs.The same stellar mass is then used to find the stellar radius using the stellar density from the photodynamical inputs (Agol et al. 2021): . (3) Finally we can take this stellar radius and combine it with our stellar-to-planetary radii ratios to get the planetary radii: Thus, by using one sampled stellar mass and Equations ( 2) through (4), a self-consistent relationship between all priors is built.For our constructed set of input chains, no single stellar, TTV, or photodynamical sample from Birky et al. (2021) and Agol et al. (2021) is used more than once; and likewise, no one input chain is used more than once in a given simulation suite (i.e., at a given initial water content).The process of sampling input chains is facilitated by VSPACE (Barnes et al. 2020), an open source support routine for VPLanet that creates a suite of input files for the model across a parameter space defined by the user; sampling from user-defined input distributions is new functionality that has been added to VSPACE in this work.Finally, Tables 1 and 2 show the 50 th percentile values and 1-σ uncertainties of the system's physical parameters, found from 20,900 input chains using the probability distributions described in Agol et al. (2021) and Birky et al. (2021).Since we fold in uncertainty in the stellar mass (as reported by Birky et al. 2021), our planetary mass and radii uncertainties reported in Table 1 may differ slightly than those reported in Agol et al. (2021), despite using identical mass and radius ratios, because they adopted a fixed 0.09 M .

Convergence
To create probability distributions of the final water loss and oxygen production on each planet per initial water content, we systematically add simulation outcomes to a distribution until convergence is achieved, i.e., when the addition of simulations to the sample no longer changes the 3σ confidence intervals of the resulting distribution.This process begins with an initial set of 500 simulations, whose final water and oxygen mass comprise the initial distribution of escape outcomes.A step is then taken, defined as the addition of 100 simulation outcomes to the distribution.The difference between the probability distributions from before and after the step is then quantified by the two-sample K-S test, which essentially determines how likely it is that the two distributions originated from the same (unknown) parent distribution.This process of taking steps is continued until the K-S statistic is 0.0004 or less, which we empirically determined to be converged, meaning the addition of more simulations no longer appreciably changes the results.As an added precaution, the model is required to report convergence for 3 consecutive steps before it declares global convergence has been achieved and the result finalized.Typically this process requires approximately 3500 -5500 simulations per planet per initial water content.This process is facilitated by Vconverge, an open source support routine for VPLanet (Barnes et al. 2020) that has been created for this work.Vconverge uses the updates made to VSPACE (described in the previous Section 2.2) to create the input files needed for the initial simulation set, and each subsequent step of 100 simulations.This support routine then automates the process described in this section to find convergence.
We run simulation suites over a modified log grid of initial water contents spanning 1 -500 TO.From 1 -20 TO, we use a step size of 1 TO; from 20 -100 TO, we use a step size of 10 TO; then from 100 -500 TO, we use a step size of 100 TO.In total, 32 initial water contents are considered.In the case of T1-e, an additional 9 simulation suites were run between 100 -200 TO, with a step size of 10 TO, to achieve the results shown in Section 3.3.1.These additional simulation suites are only considered with respect to that result, as their inclusion would not affect any other paper result.

RESULTS
We first describe the possible pathways of escape observed in our simulations for each planet through time, and how they change for different initial water contents ( §3.1).We then present probabilistic distributions of final water loss and oxygen production for each planet as a function of initial water content ( §3.2).Finally, we show the results of analyses performed with the data, including methodology for predicting initial water content from a known final water content ( §3.3.1),efficiency of oxygen production through escape ( §3.4), and a comparison to existing work with more comprehensive interior modeling (Krissansen-Totton & Fortney 2022) to contextualize the atmospheric escape calculations of both studies ( §3.5).

Escape Through Time
Here, we describe the 6 distinct evolutionary pathways that the TRAPPIST-1 planets displayed in our model, due to combinations of escape physics, assumed initial water, and the interactions between the host star evolution and planetary properties (Figure 1).A given simulation's behavior as the host star luminosity and XUV evolves is primarily contingent on initial water content and orbital distance, but noticeable effects can also be attributed to the planet's mass and radius.Understanding the possible mechanisms of escape a planet could experience as a function of stellar evolution and planetary properties helps predict likely observables.Note, that "Oxygen Retained" (used on the y-axis label of Figures 1, 2 and 3) refers to the oxygen that was produced and retained through the photolysis and thermal escape of water, when considering loss through hydrodynamic drag as the only oxygen sink.
There are 3 possible states a planet may be in when the hydrodynamic escape phase ends, here termed "end states": hydrodynamic thermal escape is halted by complete desiccation (Figure 1, left column), escape is halted by sufficient decrease in incident XUV marked by entrance to the habitable zone (HZ) (Figure 1, middle column), or escape never ends and is active at present day (Figure 1, right column).Within these cases, we can also consider whether the planet experienced diffusion-limited escape, or if it remained in the energy-limited regime for the duration of the simulation.As diffusion is responsible for mass fractionation, which may be observable (Lincowski et al. 2019), knowledge of whether or not a planet may have experienced diffusion-limited escape is important for predicting the likelihood of fractionation signals in the planetary spectrum.
Planetary desiccation is due to an interplay between the planet's initial water content and its accumulated stellar insolation.Desiccation may therefore be experienced by any planet in the system regardless of orbital distance, as those planets with lower initial water inventories are more likely to become desiccated, in addition to those with larger stellar insolation.If a planet becomes desiccated, it will pass the point where the mixing ratio of water is less than that of free oxygen (i.e., the "diffusion switch"); thus, by definition planets that become desiccated will last be in the diffusion-limited escape regime.
As an M dwarf's luminosity fades with time, a planet once interior to the habitable zone may enter it (Luger & Barnes 2015), and it is assumed that this will shut off hydrodynamic escape due to the formation of a cold trap that prevents sufficient water from reaching the exobase (Kasting 1988).Only the outer planets (beyond TRAPPIST-1 d) will enter the HZ and experience this shut-off mechanism, thus it is primarily dependent on orbital position.We find the inner edge of the HZ reaches planets T1-e, f, g, and h at roughly ∼380 Myr, ∼170 Myr, ∼100 Myr, and ∼50 Myr, respectively, thus there is a short amount of time for the planets to lose water relative to the age of the system.Consequently, depending on their initial water inventories, these planets may or may not have experienced diffusion-limited escape, with lower initial water contents more likely to reach the diffusion switch before HZ entrance.
In the last end state, escape is never halted before the age of the system is reached, and the planet retains water and active escape up to present-day.Since the planets beyond TRAPPIST-1 d will have entered or passed beyond the habitable zone by the present day, simulations that end with active escape can only occur for interior planets.This end state is primarily dependent on initial water and typically only occurs for large initial water contents of >30 -60 TO, depending on orbital distance.Planets that follow this trajectory can experience diffusion-limited escape or stay in the energy-limited regime for the entire simulation, which is also dependent on initial water content and integrated insolation over time.
As uncertainty in the stellar evolution can affect what pathway a particular planet takes, and uncertainty in planetary mass and radius can sway borderline cases, this work leverages over 700,000 simulations to explore the full uncertainty space.The predicted final states of the planets and associated uncertainty from the culmination of all simulations is described in the following section.

Statistical Outcomes of Simulation Final States
Here we present the distributions from all simulations of water loss and oxygen production for each planet as a function of initial water content (Figures 2 and 3). Figure 2 shows the predicted final water loss (top row) and oxygen retained (bottom row) for the TRAPPIST-1 interior planets (T1-b, c, and d).For low initial water contents (1 -20 TO, left column), the interior planets are completely desiccated, enabling large amounts of oxygen production; T1-c retains the most oxygen against atmospheric escape, followed by d and then b.At larger water contents, T1-c continues to retain the most oxygen, but T1-b begins to retain more than T1-d due to T1-b's greater cumulative water loss at large initial reservoirs.T1-b, c, and d begin to avoid desiccation and retain fractional amounts of an ocean at initial water contents of 60, 50, and 30 TO, respectively.
For higher initial water contents (30 -500 TO, right column), the interior planets' water loss curves begin to flatten, approaching a maximum possible water loss value.Maximum possible water loss here refers to the maximum amount of water a planet can lose given its maximum possible simulation time (the system age for interior planets, the time to the habitable zone for exterior planets) and sufficiently large water contents, such that adding more initial water would negligibly affect the predicted water loss.As water loss directly correlates with oxygen production, we also see maximum values for oxygen production and retention after thermal escape corresponding to maximum possible water loss.T1-d clearly displays this flattening off behavior indicating a maximum possible water loss of ∼90 TO and oxygen production of ∼15,000 bars.The continuing upward trend of water loss and oxygen production for T1-b and c suggest that the maximum possible water loss is first achieved for initial contents larger than those considered in this work (>500 TO), but we can conclude the maximum possible water loss of T1-b and c is greater than ∼250 TO and ∼230 TO, respectively.Escape is halted by complete planetary desiccation, for an interior planet, T1-c (black), and for a HZ planet, T1-f (orange), both with 1 TO initial water inventory.In these cases the atmospheric water is removed entirely, and the planet is last in diffusion-limited escape before the simulation halts.Middle Column: Simulations of HZ or outer planets where escape is halted by entrance to the HZ and the subsequent formation of a cold trap, all for initial water inventories of 10 TO.In this case escape could end in either the diffusion-limited (for T1-e, pink) or energy-limited (for T1-g, teal) regime.Right Column: Escape never shuts off as the planetary water reservoir is sufficiently large to sustain ongoing escape through the systems lifetime.This end state can only occur for interior planets (and results for T-1 b and d with an initial water inventory of 200 TO are shown here) as planets in or beyond the HZ have escape shut off by the formation of a cold trap.Once again, escape could end in either the diffusion-limited (for T1-b, red) or energy-limited (for T1-d, blue) regime, depending on total received insolation and initial water inventory.Top Row: Shows water lost (in TO) over time for each of the column cases previously described.Bottom Row: Shows oxygen produced (in bars) over time for each of the column cases previously described.
When considering only the balance of oxygen production and loss through thermal escape of water atmospheres, we may quantitatively explore the tendency of the planets to retain the oxygen they produce, or become completely airless.In some simulations with very low initial water, oxygen drag may be so efficient that only a negligible (<1 bar) atmosphere remains post-escape, even without the consideration of additional oxygen sinks (e.g., Schaefer et al. 2016); for an initial content of 1 TO, ∼19%, ∼2.0%, and ∼2.5% of T1-b, c, and d simulations, respectively, were completely airless.This number drops to less than 1% for all planets by an initial content of 5 TO.
Figure 3 shows the water loss and oxygen production results for the HZ and outer planets (T1-e, f, g, and h) with the same layout as Figure 2. The outer planets rarely or never experience desiccation above a threshold initial water inventory; T1-e is desiccated for initial water up to 2 TO, while T1-f is only desiccated for 1 TO initial.T1-g is nearly desiccated for 1 TO, only retaining on the order of 10 −5 TO, but T1-h is never predicted to be desiccated for the initial water contents considered here (≥1 TO).
As HZ entrance halts escape on the outer planets, there is only a short time frame prior to entering the HZ (approximately 50-380 Myr) in which these planets can experience massive water loss through thermal hydrodynamic escape, which significantly limits their maximum possible water loss, and oxygen buildup and retention.T1-e, f, g, and h display maximum possible water losses of 8.0 +1.3 −0.9 , 4.8 +0.6 −0.4 , 3.4 +0.3 −0.3 , and 0.8 +0.2 −0.1 TO, respectively, with corresponding maximum possible oxygen retention amounts of 1290 +75 −75 , 800 +40 −40 , 560 +30 −25 , and 90 +10 −10 bars, respectively.For initial water inventories above 4 TO, water loss amounts tend to decrease with increasing orbital distances, with T1-e losing the most, followed by f, g, and h.However, for low initial water contents (≤2 TO), the planets lose comparable amounts of water, but T1-g retains the most oxygen due to its lower incident radiation followed by f, e, and then h.At an initial water inventory of 3 TO, T1-g moves to retain less oxygen than T1-f and e, and by 4 TO the ordering of the planets' oxygen retention tends back towards their orbital distance ordering as the closer planets lose more water overall and thus build-up more oxygen.
Unlike the interior planets, the outer planets are rarely found to lose all produced oxygen.At an initial content of 1 TO, less than 1% of simulations for all outer planets retain <1 bar of oxygen.For all 7 planets, the differences in oxygen retention ordering are influenced by the XUV environment and planetary gravity.In Section 4, we discuss these water loss ( §4.1) and oxygen build-up ( §4.2) results, along with the broader implications and importance.

Predicting Initial Water Content
A fundamental limitation to our understanding of planetary evolutionary processes is that the atmospheric histories of exoplanets are not directly observable; however, they may be inferred from other measurable quantities.Here, we show two different approaches to inferring initial water content: first by using predictions of the present-day water mass fraction of the planets (Agol et al. 2021) ( §3.3.1), and then by using constraints on atmospheric oxygen content of the interior two planets inferred from the most recent JWST observations (Greene et al. 2023;Zieba et   Within our model and its limitations, we have shown a direct relationship between the initial and final water contents of the TRAPPIST-1 planets (Figures 2 & 3); thus, it is possible to take a known or predicted final water content and determine the initial water content our model predicts that is most likely to produce that final state.As our models only consider one evolutionary process (atmospheric escape driven by evolving stellar conditions), we caution that this analysis is meant to be a proof of concept.Our results are dependent on our model assumptions and limitations, as well as the assumptions made in the final state predictions we adopt.
In this analysis, we use the present-day water mass fraction predictions for each of the TRAPPIST-1 planets from Agol et al. (2021), which are derived from observationally constrained planetary densities, assuming an Earth-like core mass fraction of 32.5%.To calculate these water mass fractions, Agol et al. (2021) assumes all water on the planets is on the surface; steam (Turbet et al. 2020) and condensed (Dorn et al. 2017) water mass-radius relationships were used for the interior (T1-b -d) and outer (T1-e -h) planets, respectively.As any surface water on the interior planets is expected to be vaporized (Turbet et al. 2019(Turbet et al. , 2020)), they conclude the interior planets should have water mass fractions drastically less than 0.01 wt%.Since a broad range of our initial water contents would result in this desiccated state for the inner planets, we cannot constrain their most probable initial water contents from this present day desiccated state alone; thus, we only attempt to constrain initial water abundance for the outer planets.The predicted present-day water constraints we adopt for the outer planets T1-e, f, g, and h are given in Table 3 (from Table 9 in Agol et al. 2021).
The method described in Appendix C is used to determine the initial water content most likely to produce these predicted final water contents in our model.Briefly, the present-day water predictions are first modeled with a gamma probability density function (PDF); we then sum the individual simulations' values on that PDF and normalize for the number of simulations per initial water content to get a likelihood of that initial content given our predicted final state.Figure 4 shows the likelihood versus initial water content per planet.The peaks of these curves indicate the initial water content on our sampling grid that best recreated the predicted final states; for T1-e, f, g, and h these are 90, 200, 300, and 80 TO, respectively.

Body Water Mass Fraction [H2O wt%]
T1-e 2.9 +1.7 Table 3. Present-day water mass fraction constraints for the outer planets (T1-e, f, g, and h) that we use to predict most probable initial water content in our model.These constraints are taken from Agol et al. (2021) when assuming an Earth-like CMF of 32.5%, and are not definitively the true water mass fractions of the planets.2021) for each of the outer planets.The peaks denote the initial water content on our sampling grid most likely to produce the final content in our predictions.An additional 9 simulation suites were ran for T1-e only, between 100 -200 TO to fill out a smooth, continuous curve, and results from these additional simulations are denoted by triangle markers.
We further derive the 50 th percentile value of probable initial water with uncertainty for T1-e, as a proof of concept.To find the most probable initial water content predicted by our model with associated confidence intervals, we need to interpolate across the curves in Figure 4, which must possess enough data be smooth and continuous.Many of the planets we've considered do not possess smooth and continuous curves, and adding more data to these requires many additional simulation suites and computational expense; because of this, and the inherent uncertainty in our predictions, we complete this additional analysis for T1-e only as proof of concept.To obtain a continuous curve for T1-e, an additional 9 simulation suites were completed for initial water contents between 100 -200 TO (step size of 10 TO); additional data points for T1-e are shown by a triangle marker in Figures 4 & 5.These additional simulation suites are only considered in this initial water content analysis, as their inclusion in the results presented in Section 3.2 would not change the conclusions.
The interpolation for T1-e is accomplished with a modified Riemann sum, illustrated by Figure 5, to determine the 50 th percentile value of predicted initial water content and associated 1σ and 2σ confidence intervals (see Appendix C for mathematical description).The 1σ confidence interval of most probable initial water content to produce our present-day water content prediction for T1-e is 90.6   3).Same layout as Figure 4 but with linear x-axis scaling; the y-axis is the likelihood that each initial water content (in TO; x-axis) produces the desired final state.The data points that were simulated are shown by the pink points, the bins used in the Riemann sum are shown in grey, the determined 50 th percentile or true predicted value is shown by the vertical dashed black line, and the 1σ and 2σ confidence intervals are shown by the shaded dark and light blue regions, respectively.The low initial water contents (1 -20 TO) are much more finely spaced (step size of 1 TO) which is why the Riemann sum bin widths are not resolvable there.Lines connecting data points are used for the tops of the Riemann sum bins.

Using Constraints on Atmospheric Oxygen Abundance
Although the statistical inference technique ( §3.3.1) can produce an upper limit on initial water abundance, a more refined estimate can be obtained by deriving tighter constraints from existing and future observations of atmospheric oxygen abundance on T1-b and c.Recently, Greene et al. (2023) used JWST/MIRI to obtain a blackbody brightness temperature of T B = 503 +26 −27 K for T1-b from secondary eclipse observations taken at 15 µm, which they find to be an excellent match to the 508K theoretical prediction of dayside temperature for an airless world with no heat redistribution.If a significant atmosphere was present, then the subsequent redistribution of heat from the dayside to nightside, and/or absorption from atmospheric CO 2 , would result in an average dayside temperature at 15 µm that is significantly lower than that observed.Greene et al. (2023) therefore conclude that T1-b likely does not have a substantial atmosphere today.This observation is also consistent with nearly half of all evolutionary model runs from Krissansen-Totton & Fortney (2022), which predict complete atmospheric erosion on T1-b across their considered range of initial water (0.7 -300 TO).
The observational constraint that T1-b is likely airless today can be used with our grid of atmospheric evolutionary outcomes to identify limits on the planet's initial surface water.Considering only the apparent lack of water, our model predicts the planet must have begun with ≤60 TO to be in a desiccated present-day state.Although water does not remain in these simulations, an oxygen atmosphere could still be supported, especially for the higher end of the initial water range, and depending on assumed sinks.However, a dense O 2 atmosphere, with even a small amount of outgassed CO 2 , could produce an absorption feature at 15µm that is also precluded by the observational constraints (Greene et al. 2023;Ih et al. 2023;Lincowski et al. 2023).The tighter observational constraint that T1-b is currently desiccated and also likely does not have an oxygen atmosphere, combined with considerations of oxygen sinks over the lifetime of the system, will significantly lower this upper limit on initial water inventory.
Though detailed oxygen surface sink modeling is outside the scope of the current study, to give an example of how one may further constrain this upper limit we estimate a few oxygen sinks here, first considering only oxygen loss from the magma ocean and hydrodynamic escape while water loss is ongoing.Assuming magma ocean removal operates on the same timescales as hydrodynamic escape (Hamano et al. 2013), we estimate a removal of 3 TO of oxygen to the magma ocean (∼630 bars on T1-b), based on studies of the early Venusian environment and theoretical studies on magma ocean evolution for the TRAPPIST-1 planets (Gillmann et al. 2009;Schaefer et al. 2016;Wordsworth et al. 2018;Barth et al. 2021); and from our simulations we find 4 TO equivalents of oxygen would be lost to hydrodynamic drag.These combined give maximum initial surface water contents of 7 TO, to be consistent with an airless present day T1-b.However, loss of atmospheric oxygen can still continue after total loss of the water inventory, via processes such as non-thermal escape and dry crustal oxidation.Following Krissansen-Totton & Fortney (2022), we consider a non-thermal escape parameterization of ∼100 bars of total oxygen loss across the system's lifetime, and an average estimate of oxygen loss through dry crustal oxidation of ∼10 Tmol/year (∼55 bars/Gyr) after magma ocean solidification (following the perpetual runaway greenhouse scenario in Krissansen-Totton et al. 2021b).Given a system age of ∼8 Gyr (Birky et al. 2021) and initial water contents of around 15 TO, T1-b may endure hydrodynamic escape for up to 1 Gyr, leaving 7 Gyr for dry crustal oxidation; thus, we conservatively estimate 385 bars of total oxygen removal through dry crustal oxidation.With these two additional oxygen sinks (magma ocean removal and dry crustal oxidation), our upper limit constraint on initial surface water, consistent with complete atmospheric erosion on T1-b at the present day, would be lowered to 12 TO.However, note that with more rigorous modeling of oxygen removal, this upper limit would likely be higher.Overall we can conclude the upper limit on initial water content as constrained by observations of T1-b (Greene et al. 2023) is between 12 -60 TO.
Although the above upper limit constraints on initial water abundance are developed with the observational constraint from TRAPPIST-1b, habitability studies of the system may be more concerned with the lower limit of initial water, which could inform the minimum amount of water the HZ planets possess at present-day.Constraining this lower limit may be possible by considering what we know of the atmospheric state of TRAPPIST-1c in addition to that of T1-b.Recent T1-c secondary eclipse observations with JWST/MIRI show a brightness temperature of T B = 380±31 K (Zieba et al. 2023).This is 50K below, and 1.7-σ away from, the temperature expected for an airless body with zero heat redistribution (430 K).It is also 1.6-σ from an ultramafic surface with the strong weathering expected for this intensely UV-irradiated world (Zieba et al. 2023).The T1-c JWST 15µm data are best fit by either an airless world with a high albedo surface, or a thin <0.1 bar O 2 -dominated atmosphere (Zieba et al. 2023;Lincowski et al. 2023).However, other atmosphere types are still consistent with the data, including a 3 bar steam atmosphere, which is within 1.7-σ of the observed secondary eclipse value (Lincowski et al. 2023).
Given their different stellar distances, similar planetary properties and inferred geologic oxygen sinks, there is a narrow range of initial surface water contents that will simultaneously produce an airless T1-b and non-negligible oxygen on T1-c, even if both planets are found to be desiccated.As T1-b and T1-c are similarly interior to the HZ, possess comparable surface gravities (g T 1c ∼ 0.99g T 1b ), and likely formed with a similar composition (ρ T 1c ∼ 1.004ρ T 1b ) they likely experience comparable geologic and escape evolution.However, we have shown that since T1-c is slightly farther away from the host (a T 1c ∼ 1.37a T 1b ), it will retain more oxygen for equivalent amounts of water loss as T1-b due to the lower incident XUV flux (see §3.4 and Figure 7).
Considering oxygen removal through a magma ocean, dry crustal oxidation, and non-thermal escape, we run an Markov chain Monte Carlo (MCMC) to constrain the initial water content that would produce both an airless T1-b and a tenuous 0.1 bar O 2 -dominated atmosphere on T1-c, consistent with the currently available JWST data (Greene et al. 2023;Zieba et al. 2023;Ih et al. 2023;Lincowski et al. 2023).To conduct the MCMC, we use the publicly available emcee model (Foreman-Mackey et al. 2013).MCMC walkers fit for the initial water content by running VPLanet (Barnes et al. 2020) simulations of T1-b and T1-c that predict the oxygen build-up of the system.After a simulation, the estimated oxygen removal is subtracted, and the MCMC checks to see that T1-b has no oxygen remaining while T1-c is compared to a 0.1 bar oxygen atmosphere with an uncertainty of 10 bars, as a 10 bar oxygen atmosphere begins to show >3σ disagreement from the 15 µm measurement of T1-c (Zieba et al. 2023).As described earlier in the section, we estimate removal of 3 TO-worth of oxygen from a magma ocean (Gillmann et al. 2009;Barth et al. 2021), 100 bars from non-thermal escape (Krissansen-Totton & Fortney 2022), and 10 Tmol/Gyr (∼55 bars/Gyr) from dry crustal oxidation after the hydrodynamic escape period (Krissansen-Totton et al. 2021b).We vet convergence by considering the chain's autocorrelation time; in total, we run 20 walkers for 80,000 steps with 1,000 burn-in steps (for a total of >1.5 million VPLanet simulations), this results in a chain length that is several times the autocorrelation The vertical solid black lines on both panels denotes the 50 th percentile value predicted by the MCMC while the vertical dashed black lines denote the 68% (1σ) confidence interval.This analysis predicts a minimum initial water content of 8.2 +1.5 −1.0 TO to produce the desired final state, assuming both planets begin with the same initial water content.The amount of oxygen removal is meant to be an underestimate to produce a conservative lower limit, as additional oxygen removal would increase the minimum initial water required.time, suggesting convergence has been achieved (Foreman-Mackey et al. 2013).We note that our oxygen removal parameterizations are meant to be low estimates of the total oxygen removal the planets may experience, as we are searching for a minimum initial surface water mass and greater oxygen removal would serve to raise the minimum.
Figure 6 shows the results of our MCMC analysis; the left panel is an illustration of the atmospheric oxygen left after removal by the estimated oxygen sinks, and the right panel is the distribution of initial water contents reported by the MCMC that best fit an airless T1-b and a T1-c with a 0.1 bar oxygen atmosphere.With 68% confidence (1σ) we find the minimum initial water content to be 8.2 +1.5  −1.0 TO.This value with associated uncertainty is shown by vertical lines on both panels of figure 6.
As one last note, in considering the outer planets, under the assumptions of a similar magma ocean removal rate during hydrodynamic escape, followed by a removal of ∼135 bars/Gyr analogous to geological O 2 sinks on the modern Earth (including weathering, volcanism, and subduction) (Catling 2014), atmospheres lacking in O 2 are possible for all initial water contents considered.We point out, however, that these are rough estimates, and further work would be needed in tracking oxygen sinks over time to produce more accurate results.

Oxygen Retention Efficiency
Quantification of the oxygen retention efficiency allows us to evaluate the net balance between oxygen produced via water photodissociation and subsequent hydrogen escape, and oxygen loss due to drag during thermal hydrodynamic escape.Oxygen retention efficiency is given by: where M O2,produced is the total mass of oxygen that was produced by the cumulative hydrogen loss and M O2,lost is the mass of oxygen that was lost through drag during escape.Note that this only takes into account production and losses via escape mechanisms, and oxygen may be further lost through sinks not considered here (e.g., Schaefer et al. 2016;Schaefer & Fegley 2017).
Figure 7 shows the oxygen retention efficiency (Eq.5) as a function of initial water content with contours for the 50 th percentile values and 1σ confidence region for all 7 planets.The color scale of Figure 7 shows, for a particular initial water content, the percentage of simulations that had a given oxygen retention efficiency.Additionally, Table 4 gives the 50 th percentile values and 1σ uncertainties for oxygen retention efficiency for our lowest (1 TO) and highest (500 TO) initial water case for all 7 planets.This result shows that oxygen retention efficiency increases with increasing orbital distance, initial water content, and planetary gravity.At low initial water contents, there is a strong distance dependence where efficiency increases with increasing orbital separation due to the drop in XUV flux; however, planets with low gravity may break this distance trend and show lower than expected O 2 retention efficiency when compared to the distance trend, such as T1-d and h.High initial water inventories (100 TO or more) typically lead to larger retention efficiencies than those for lower water inventories, and these retention efficiencies display a weaker distance dependence; for example T1-c at its largest initial water content possesses the largest retention efficiency of any planet.This is because high initial water contents sustain escape to much later ages; as oxygen loss through drag will be higher early on with a higher XUV environment, decreasing with time as XUV decreases, the cumulative oxygen retention efficiency will be higher for planets that remain in active escape to later times.In comparison, planets where hydrodynamic escape shuts off earlier will have lower cumulative oxygen retention efficiency because all loss processes occur in a high XUV environment.This influence is further shown by the plateau of retention efficiency as planets reach initial water contents that lead to their maximum possible water loss (i.e., the maximum amount of water they may lose given the maximum possible simulation time and an effectively infinite initial reservoir), indicating that retention efficiency is directly dependent on time spent in the escape regime with evolving XUV conditions, which is in turn dependent on initial water.At low initial water contents, the interior planets (T1-b, c, and d) retain ∼10% -30% of oxygen produced, while the outer planets retain ∼40% -70%.At high initial water contents, this number increases to ∼70% -90% for all planets in the system, excluding T1-h, which has extremely low gravity and loses ≲1 TO for any initial content, thus displaying retention efficiency that remains at ∼46% even as initial water increases.
However, counter to the general trend of retention efficiency increasing with initial water content described above, we also find that there is a limited range of initial water contents for each planet where feedback from diffusion processes causes retention efficiency to decrease with increasing initial water; this decreasing behavior is particularly noticeable for T1-d near ∼100 TO and T1-e near ∼10 TO in Figure 7.For a particular planet, this decreasing behavior occurs when the initial water content approaches an amount large enough to prevent the planet from entering the diffusion-limited escape regime over its escape lifetime (either prior to entering the HZ, or up to present day).Since diffusion-limited escape prohibits oxygen loss and favors oxygen retention, avoiding a diffusion-limited escape period will decrease oxygen retention efficiency.Thus, when this regime is prevented entirely, we see a range where retention efficiency decreases with increasing initial water, as the influence of the lack of diffusion-limited escape outweighs the gain in efficiency from increasing the mixing ratio of water.However, this behavior only occurs for a narrow range of initial water abundances, and increasing efficiency with increasing initial water contents is seen on either side of that range.
When losing equal amounts of water a strong distance dependence persists of increasing retention efficiency with increasing orbital distance, even when increasing initial water content; this behavior is apparent by the results in Table 5, which show the efficiency of oxygen retention following an equal loss of 1 TO for a low (1 TO) and high (9, 10, or 11 TO) initial water case.Note, as this result compares an equal water loss, simulations that do not lose a full TO are not included.Consequently many cases for T1-g, and all cases for T1-h are excluded, as the former does not typically lose a full TO when it has a low starting inventory, and the latter does not typically lose a full TO for any starting inventory considered.Overall, decreases in orbital distance, initial water, and gravity will all decrease oxygen retention efficiency.

Comparison to Atmosphere-Interior Evolution Model
Here we compare our calculated escape rates to Krissansen-Totton & Fortney (2022), a similar study modeling the atmospheric evolution of the TRAPPIST-1 planets, but with a particular focus on the interior evolution and surface-atmosphere exchange.Krissansen-Totton & Fortney (2022) use the Planetary Atmosphere, Crust, and MANtle (PACMAN) evolutionary model to simulate terrestrial planet mantle and atmosphere evolution from an initial magma ocean phase through solidification to temperate geochemical cycling on the TRAPPIST-1 planets (validated in Krissansen-Totton et al. 2021a,b).They complete 5000 model runs total for each planet; unlike this study, their simulations sample initial water in log space from 0.7 -300 TO.For the interior planets, T1-b, c, and d, they show  and h.Coloring shows the percentage of simulations at a particular initial water content that fell in a given oxygen retention efficiency bin (e.g., each column corresponds to one initial water content and sums to 100%).Contours are given for the 50 th percentile values (solid white lines) and the 1σ confidence region (dashed white lines).half of all model runs produce >1 bar oxygen atmospheres while the other half result in airless bodies.Alike to the current study, oxygen-rich modern atmospheres become less probable as distance from the star increases, noting that while initial water has the largest effect on oxygen build-up for the inner planets, atmospheric escape physics and geological factors influence the outer planets more heavily.Overall, Krissansen-Totton & Fortney (2022) show that complete atmospheric erosion or oxygen-rich atmospheres are most likely for the inner planets, while the outer planets are more likely to possess CO 2 or CO 2 /O 2 dominated atmospheres, the latter occurring if CO 2 efficiently acts as a secondary producer of oxygen through photodissociation.2022), and our median test case with an initial water inventory of 15 TO (log space average between 0.7 and 300 TO) stays within their model uncertainties for the entire age of the system.Note, in the T1-e comparisons, our 15 and 100 TO cases display nearly identical evolution and overlap on the plot.

Oxygen Retention Efficiency
Figure 8 shows a comparison of the rates of H 2 O loss and oxygen build-up with 2σ uncertainty envelopes, as presented in Krissansen-Totton & Fortney (2022), with our work for T-1b, an interior planet, and T-1 e, a HZ planet.Both studies adopt the stellar evolution models of Baraffe et al. (2015), and XUV properties from Birky et al. (2021), meaning differences in our work are not due to dissimilar XUV environments.Since we consider initial water contents on a grid, rather than the randomly sampled input of Krissansen-Totton & Fortney (2022), we compare 3 of our initial water contents in Figure 8; we show initial contents of 1 and 100 TO to bracket the median water contents from Krissansen-Totton & Fortney (2022), and an initial content of 15 TO is shown as the average value of 0.7 and 300 TO in log space.
The agreement between our studies is favorable, with our median 15 TO case remaining within the 2σ uncertainty region of Krissansen-Totton & Fortney (2022).The average rates from Krissansen-Totton & Fortney (2022) and from all shown initial water contents in this study begin in near perfect agreement (≲10 7 years), with our 1 TO case deviating first as desiccation is reached quickly.For T1-b, our 100 TO case begins to deviate more at ∼1 Gyr, showing relatively larger fluxes of both water escape and O 2 production as more water is available to sustain escape than the majority of model runs from Krissansen-Totton & Fortney (2022).Our 15 TO T1-b case remains in close agreement for the entirety of the system's lifetime, although the shape of the curves show that escape rates in Krissansen-Totton & Fortney (2022) slow more gradually.For T1-e, entrance to the HZ strictly ends evolution for our 15 and 100 TO cases at ∼380 Myr.However, as Krissansen-Totton & Fortney (2022) did not use as stringent of a simulation halt, they still report minor water loss and oxygen production past this point, although it is significantly reduced.Before HZ entrance, our 15 and 100 TO T1-e cases behave similarly, but the 15 TO case does begin to deviate and show more quickly lowering rates suggesting it would remain in better agreement with Krissansen-Totton & Fortney (2022), as expected.In general, 2σ uncertainties on the results shown in this study are smaller because the curves each have a fixed initial water content, while random sampling of initial water content in Krissansen-Totton & Fortney (2022) contributed to larger uncertainty in simulation behavior.
When oxygen drag is efficient, a "negative gain flux" may occur, where oxygen is lost faster than it is photolytically produced, which could deplete existing atmospheric oxygen inventories.Figure 8 shows these negative O 2 gain fluxes early in the age of the system (≲200 Myr), which we were only able to reproduce with VPLanet after the model modifications implemented here.We found that closer-in planets with lower initial water contents are more likely to experience a negative O 2 gain flux.For example, with an initial water content of 1 TO for T1-b and T1-e, 62% of the T1-b and 35% of the T1-e simulations experienced at least one time step with a negative O 2 gain flux, respectively; this number dropped to 16% and 4%, respectively, when initial water was increased to 2 TO.
Negative O 2 gain fluxes occur when the mixing ratio of oxygen is as high as possible without pushing escape into the diffusion-limited regime, and must occur during periods of high incident XUV (e.g., ≲200 Myr).As the mixing ratio of oxygen increases, loss to hydrodynamic drag also increases (e.g., Eq.A8); but if the mixing ratio of O 2 exceeds that of water, the atmosphere will enter diffusion-limited escape and oxygen drag will be shut-off altogether.Negative O 2 gain flux is therefore most likely to occur during periods when the incident XUV is high (e.g., ≲200 Myr), and when the mixing ratio of oxygen approaches the maximum value needed to trigger diffusion limited escape and the shut down of oxygen drag.A negative O 2 gain flux will then continue until diffusion-limited escape starts or the XUV flux has decreased sufficiently.Considering our 1 TO cases in Figure 8, T1-b starts to lose net atmospheric O 2 almost immediately due to its close proximity to the star, whereas T1-e needs at least 10 Myr to enter a negative O 2 gain flux state, as some time is needed for the oxygen mixing ratio to increase sufficiently.

DISCUSSION
We have used the most recently available constraints on the TRAPPIST-1 system (Agol et al. 2021;Birky et al. 2021) and a hydrodynamic atmospheric escape model to statistically explore the ranges of possible water loss and oxygen buildup on the TRAPPIST-1 planets.Improving on past studies, we combine observational uncertainties on the stellar luminosity evolution and planetary parameters to construct self-consistent input chains that are randomly sampled for each simulation, fully incorporating the statistically plausible region of initial conditions for physical system parameters.Without considering outgassing, we find the interior planets (T1-b, c, and d) nearly, or completely, lose all of their initial surface water for all but the largest initial contents (when initial water is >50 -100 TO), while the outer planets (T1-e, f, g, and h) only become desiccated for the very lowest initial water contents considered (1 -2 TO).Given that our model halts hydrodynamic escape when a planet enters the habitable zone, the outer planets never lose more than 10 TO total, which results in maximum possible oxygen retention values spanning ∼90 -1290 bars.
Comparing our results on planetary water loss and oxygen production rates to those reported in Krissansen-Totton & Fortney (2022), we find that rates from an example set of 3 of our tested initial water contents remained within their reported uncertainties for most, or the entirety, of the age of the system, showing favorable agreement.
We find that atmospheric oxygen retention -what fraction of the oxygen produced by water photolysis is retained in the atmosphere -is a function of both planet-star separation and initial water inventory.Closer to the host star, we find that larger percentages of photolysis-produced oxygen are lost due to hydrodynamic drag, offsetting the substantial oxygen produced by water photolysis.For the cumulative water loss at low initial contents, our simulations show that the interior planets (T1-b, c, and d) retain ∼10% -30% of the oxygen generated, and the outer planets (T1-e, f, and g) retain ∼50% -70%.At high initial water contents, this retention efficiency approaches ∼70% -90% for all planets in the system, excluding T1-h.However, despite the enhanced loss rate of oxygen, for initial water contents of more than 4-5 TO, the overall production rate is sufficiently high that the interior planets always retain more O 2 than the outer planets.Consequently, for initial water contents ≲4 TO, the outer planets are more likely to retain atmospheric oxygen.
We retrieved the suite of plausible water evolutionary histories consistent with a known or predicted final water mass fraction (WMF) for the outer planets, demonstrating a mechanism to determine the most probable initial water content of the planet for a predicted final water content.For the maximum present-day WMFs constrained by measured planetary densities Agol et al. (2021), we find the most probable initial water contents for all outer planets to be in the range of ∼80 -300 TO.For T1-e, we interpolated from this probability data to report the 1σ confidence interval on our maximum initial water content prediction of 90.6 +70.6  −37.2 TO, based on the model we employ and the present-day water constraints we adopt.This analysis is a proof of concept for a method to retrieve probable initial water content from an assumed present-day state.However, our model assumptions will produce upper limits, and the true initial water content may be much lower, as we assume that 1) the planets' low densities are entirely due to their WMF, and 2) that water was not sequestered in the planetary interior, but was always on the surface and accessible to atmospheric loss.To elaborate on the latter point, if water was instead partially or fully trapped in the planetary interior, less initial water would be required to recreate the same WMF at present-day, as it would be protected from loss due to escape.

Behavior of Water Loss
Our work shows that water loss rates and cumulative water lost are highly dependent on the incident XUV flux and the atmospheric mixing ratio of water; implying both orbital distance and initial water content exert strong influence on the end state observed.Qualitatively, increasing orbital distance will decrease water loss by decreasing incident XUV, and decreasing the time to HZ entry and hydrodynamic escape shutoff.Conversely, increasing initial water content will increase water loss rates by creating a higher atmospheric water mixing ratio at earlier ages during greater incident XUV.Across the conditions considered in this work, the interior planets (T1-b, c and d) frequently lose their entire initial water inventory, while the outer planets (T1-e, f, g and h) are more likely to retain water to the present-day.
As the interior planets are located close enough to the star to experience escape for the entirety of the system's lifetime, their end state is most heavily influenced by initial water content (see Figure 1), which is in agreement with (Krissansen-Totton & Fortney 2022).Conversely, the end state behavior of the HZ or outer planets ( T1-e, f, g, and  h) is primarily influenced by distance from the host star, which affects total integrated XUV flux experienced by the planet, and also determines the time to habitable zone entrance (i.e., total time spent in the hydrodynamic escape regime).For the outer planets, escape is shut off at HZ entrance due to the assumed formation of a cold trap; however, we caution that the true efficacy of a cold trap to prevent hydrodynamic escape may not be this efficient, and may be influenced by other factors (see also §4.4).As HZ entrance results in a short period of escape (∼50 -400 Myr) relative to the age of the system (τ ∼8 Gyr) for the outer planets, desiccation is rarely witnessed and the maximum amount of water they may lose through hydrodynamic escape is <10 TO in all cases in this study, and as low as 0.8 TO for T1-h.Though the outer planets may better constrain initial water content when considering measurements of present-day WMF, measuring the present-day WMF is observationally challenging as water may be sequestered in the interior or in a surface ocean, and so not amenable to transmission observations, and also difficult to derive uniquely from planetary density constraints (Agol et al. 2021).However, initial water more heavily influences the present-day atmospheric oxygen abundance on the interior planets, and ocean-loss generated O 2 , or its proxy O 3 , may be accessible in JWST observations (e.g., Lustig-Yaeger et al. 2019) and next generation ground-based telescopes (e.g., Currie et al. 2023), and so may provide the first strong constraints on the initial water content of the planetary system as a whole, when considering current and near-term observational capabilities.

Comparison with Previous Work
Comparison with past studies shows agreement in our predictions that the inner planets are more likely to be desiccated, but the outer planets retain water, and that stellar XUV and activity evolution have larger impacts on total water loss than geological processes.Our results confirm numerous past studies that desiccation is possible on the interior planets of TRAPPIST-1, or a similar M dwarf system, for all but extremely large initial water contents (Luger & Barnes 2015;Schaefer et al. 2016;Bolmont et al. 2017;Bourrier et al. 2017;Johnstone 2020;Krissansen-Totton & Fortney 2022).Similarly, we also find maximum water loss on the HZ and outer planets of ∼1 -10 TO, implying these planets may still possess water at present-day (Luger & Barnes 2015;Bolmont et al. 2017;Bourrier et al. 2017;Barth et al. 2021;Krissansen-Totton & Fortney 2022).Comparison with past studies shows that predicted water loss is very sensitive to assumptions on the stellar XUV and activity evolution, implying that tighter constraints on TRAPPIST-1's stellar spectrum and activity can greatly improve predictions for planetary outcomes.For example, Bolmont et al. ( 2017) (considering T-1b and c) and Bourrier et al. (2017) (considering all 7 planets in the system) adopt XUV saturation fractions ∼1.32x and ∼2.5 -3.5x smaller than this work, respectively.Consequently, water loss predictions for the inner planets in this study are up to 3x larger than Bolmont et al. (2017) and >3x larger than Bourrier et al. (2017).Comparison of our results, which do not include geological processes, with those of Krissansen-Totton & Fortney (2022), that do, suggest that the inclusion of geologic processes (e.g., sequestration of water in the planetary interior and outgassing over time) has less of an influence on total water loss than changing XUV environment, since Krissansen-Totton & Fortney (2022) adopts the same XUV parameterizations (Birky et al. 2021) as this study and reports agreement on overall conclusions on the final planetary water contents.These comparisons highlight how marginal increases in XUV flux may produce drastic increases in water loss, and suggest more observations of the TRAPPIST-1 stellar environment may translate to significant improvements in model predictions.

Oxygen Production, Retention and Loss
The upper limits on oxygen production and retention quantified by this study provide important environmental context on atmospheric O 2 in the TRAPPIST-1 system, and while retained atmospheres less than 1 bar may constitute a false positive for an O 2 biosignature (Schaefer et al. 2016;Meadows 2017), significantly larger amounts may be a clear indication of an abiotic post-ocean loss atmosphere (Luger & Barnes 2015).On Earth-sized terrestrials, the loss of an ocean's worth of water will generate ∼250 bars of O 2 , and the inner planets, which lose more water over the system lifetime, can produce hundreds to thousands of bars atmospheric oxygen, depending on initial water content.However, not all this oxygen is retained, and significant losses can occur due to energetic hydrogen escape flows that remove the oxygen from the atmosphere.We find that cumulative oxygen retention increases with increasing distance from the star; however, for initial water contents greater than 4 -5 TO, the interior planets retain more oxygen cumulatively than the outer planets due to greater overall water loss (and the subsequent increase in oxygen production).
Similarly, we find the efficiency of oxygen retention increases with increasing distance from the star and/or increasing initial water content ( §3.4, Figure 7).The distance dependence is due to greater XUV flux at small orbital separations, resulting in more energetic escape flows and oxygen loss through drag; however, planets with low gravity, including T1-d and h, may deviate from this pattern and lose more oxygen than their closer-in neighbors.The increase in the efficiency of oxygen retention with increasing initial water content is driven by larger water inventories producing a large atmospheric water to oxygen ratio, which decreases the rate of oxygen loss through drag.
Low oxygen retention efficiencies are also more likely to occur when oxygen mixing ratios are higher, but before the oxygen fraction exceeds the water fraction in the atmosphere, triggering diffusion limited escape and oxygen retention.Extremely low initial water contents can lead to high oxygen mixing ratios early in the system's lifetime, during the greatest incident XUV, resulting in highly efficient oxygen drag and lowered retention efficiency -in some cases leading to total oxygen loss and atmospheric erosion for the interior planets.A similarly lowered retention efficiency could also occur if the initial planetary water inventory was high, but sequestered in the interior and slowly outgassed, rather than being continuously available for escape at the surface as is often assumed (e.g., Krissansen-Totton & Fortney 2022).In this case, the throttling of the water availability by outgassing would be more likely to result in lower atmospheric water mixing ratios at any given time leading to greater oxygen loss through drag.

Comparison with Previous Work
In contrast to the general agreement on water loss predictions we see across previous studies ( §4.1.1),predictions of oxygen retention for M dwarf planets vary.In this work, we witness a broad range of oxygen retention scenarios depending primarily on initial water content and orbital distance.For medium to high initial water inventories (typically above ∼5 TO) we predict comparatively high planetary oxygen retention (e.g., hundreds to thousands of bars), which is in agreement with several previous studies (Luger & Barnes 2015;Bolmont et al. 2017;Bourrier et al. 2017;Krissansen-Totton & Fortney 2022).In contrast, others show that little to no oxygen would survive hydrodynamic escape processes (e.g.zero to tens of bars) even at medium to high initial water contents (Schaefer et al. 2016;Johnstone 2020).Within these predictions, there is further disagreement as to the precise quantity of oxygen retained, and if it would remain in the atmosphere (e.g., Luger & Barnes 2015) or be removed partially or entirely by reactions with iron in the planet's mantle (e.g., Wordsworth et al. 2018;Barth et al. 2021;Krissansen-Totton & Fortney 2022).
At low water contents (≲4 TO), this work indicates that large percentages, or even the entirety of photolytically produced oxygen may be lost solely through drag during hydrodynamic escape (as demonstrated in Figures 3, 7, and  8); similar results with complete oxygen loss have been reported by several previous studies (Schaefer et al. 2016;Johnstone 2020;Krissansen-Totton & Fortney 2022).Schaefer et al. (2016) considers water contents up to 1000 TO in the GJ 1132 system (M * = 0.181 M , see Berta-Thompson et al. 2015) and show that the super-Earth GJ 1132 b would likely lose 90 -99% of all photolytically produced oxygen to space through drag.Likewise, the Schaefer et al. (2016) model was also adopted by Kreidberg et al. (2019) which reported that initial water contents <240 TO resulted in complete atmospheric loss on LHS 3488 b.Comparing GJ 1132 b and LHS 3488 b to T1-b in the TRAPPIST-1 system, they have substantially larger insolation -with GJ 1132 b receiving about 19 times more flux than Earth (Schaefer et al. 2016) and LHS 3488 b receiving 70 times more flux than Earth (Vanderspek et al. 2019), whereas T1-b receives only about 4 times more than Earth (Agol et al. 2021).This higher incident radiation helps to explain why such high oxygen loss rates are reported.Similar oxygen behavior was shown by Johnstone (2020) when considering an Earth-like planet with an M dwarf host, reporting water loss up to 120 TO during the stellar saturation time (the first 1 Gyr) with little to no planetary oxygen retention as a result of efficient drag.However, when considering medium to high initial water contents (>4 -5 TO) and also accounting for drag, our results instead predict a significant build-up of oxygen.This contradicts the work of Johnstone ( 2020), but is consistent with several other studies, which predict oxygen build-up ranging from several hundreds (Bolmont et al. 2017) to thousands of bars (Luger & Barnes 2015) on the interior planets, and several tens (Bolmont et al. 2017) to hundreds of bars (Luger & Barnes 2015;Barth et al. 2021) on the outer planets.These differences in predictions are likely because, as we have shown, the precise amount of build-up is highly dependent on duration of escape, planetary gravity, and mathematical descriptions of crossover mass and oxygen flux.Observations to constrain atmospheric oxygen on the TRAPPIST-1 planets would help constrain evolutionary scenarios and outcomes, and improve theoretical models.

Constraining the Early Water Abundance of the TRAPPIST-1 Planets
We explore the process of constraining initial water content of the TRAPPIST-1 planets, both through model inference ( §3.3.1) and by considering recent JWST measurements.By combining simulations and JWST data from multiple planets -in this case, TRAPPIST-1 b and c -we can derive much tighter constraints on initial surface water abundance ( §3.3.2).

Statistically Retrieving Planetary Initial Water Content
As we show in Section 3.3.1,given present-day water and oxygen inventories, our model suite of escape histories can provide statistical constraints on initial water content and past water loss.While using constraints on present day water inventories informed by measured planetary densities can provide a key constraint (Agol et al. 2021), future observational evidence on the presence and nature of the TRAPPIST-1 atmospheres could be used to refine initial water content and subsequent planetary evolution.For example, if an outer planet is found to be desiccated, it would place a strict upper limit on initial water content, given the much shorter escape lifetimes for outer planets.Correspondingly, desiccation on an interior planet provides weaker constraints as there is a much broader range of initial water abundances that could still result in desiccation over the system lifetime ( §4.3.2).For a desiccated interior planet, oxygen may be a better indicator of ocean loss and other evolutionary processes as it is much more tightly correlated to initial water content (see Figure 2).If a planet is observed to still possess water, the method explored in Section 3.3.1 and Appendix C demonstrates a basic framework for conducting statistical inference to achieve an initial upper limit constraint on initial water.
Statistical inference may also be leveraged as a model intercomparison tool, and to provide critical constraints on planet formation, migration and water delivery models.In the example of the current work, if multiple studies on water loss conduct an analysis to map initial to final water content, we can indirectly compare the latent space (e.g., exact equations, order of operations, timestepping, etc) predicted by those models during evolution.In a more concrete sense, our method of retrieving initial from final water content applied to the TRAPPIST-1 system can provide model constraints for formation and water delivery scenarios of the system.Even with relatively uncertain constraints on initial water content, this may inform if the planets formed in situ or past the snow line followed by inwards migration (Raymond et al. 2022).Recent observations of TRAPPIST-1b and c have demonstrated that the former may be an airless world (Greene et al. 2023;Ih et al. 2023), while the latter may be better fit by a tenuous oxygen-dominated atmosphere (Zieba et al. 2023;Lincowski et al. 2023).In Section 3.3.2,we used MCMC to fit for the initial water content that would produce this final atmospheric state when assuming basic oxygen sink estimations and that the two planets formed with comparable water inventories; the results of this analysis are shown in Figure 6.Our MCMC reported that the minimum initial water content of the planets is 8.2 +1.5 −1.0 TO.This estimated range is in agreement with Zieba et al. (2023), who suggest an initial water mass of 9.5 +7.5  −2.3 TO through an analysis of conservative CO 2 upper limit constraints, albeit with a similar atmospheric escape model (Luger & Barnes 2015;Schaefer et al. 2016).
These constraints on minimum initial water using JWST observations of T1-b and c, and our estimates of subsequent water and oxygen escape and loss, have important ramifications for the presence of water, and habitable conditions, on the TRAPPIST-1 habitable zone planets.If all planets in the system formed with comparable initial water contents, as supported by the delicate resonance chain and recent formation studies (e.g., Raymond et al. 2022), our simulations suggest that an initial water content of 8 TO would imply the HZ and outer planets would still possess water after escape is shut off upon entrance to the habitable zone.With an initial water content of 8 TO, we find T1-e, f, g, and h would possess 1.7 +0.47  −0.57, 3.6 +0.32 −0.46 , 4.8 +0.21 −0.29 , and 7.2 +0.14 −0.20 TO, respectively, after HZ entrance.Thus, a confirmation of oxygen in the atmosphere of T1-c would increase the probability that the HZ planets retained water and could be potentially habitable.
If T1-c is airless or nearly airless, this still may provide information on an upper limit on initial water of the system, for example this would narrow our previous upper limit constraint from T1-b alone; following the same removal estimates described in section 3.3.2,we find an airless T1-c would suggest an upper limit on initial surface water of 9 TO.Alternatively, in this airless T1-c scenario, Zieba et al. (2023) predict initial water content would be reduced to ∼4 TO; for this initial water, our model predicts T1-e, f, g, and h would possess 0.03 +0.05 −0.03 , 0.52 +0.15 −0.17 , 1.2 +0.14 −0.19 , and 3.22 +0.14  −0.20 TO, respectively, after HZ entrance.Although this case would imply a greater risk of water loss on T1-e and f, it would not rule out surface water on these planets and would still support considerable endowments on T1-g and h.This is in agreement with the conclusion of Krissansen-Totton (2023) that an airless T1-b (and c) does not preclude the outer planets from possessing substantial atmospheres.
If the TRAPPIST-1 HZ and outer planets are likely to possess present-day oceans, as the recent JWST observations (Greene et al. 2023;Zieba et al. 2023) interpreted with our models suggest, then these will make important targets for JWST observations.When observing these targets in transmission, however, the further work will be required to remove the effects of stellar contamination from the data, as recent observations have shown to be prevalent (Rackham et al. 2018;Lim et al. 2023).Though O 2 has been shown to be extremely difficult to detect with JWST, particularly for the outer planets, which will be restricted to transmission measurements (e.g., Lustig-Yaeger et al. 2019;Wunderlich et al. 2019;Gialluca et al. 2021Gialluca et al. , 2023;;Meadows et al. 2023), post-ocean loss O 2 atmospheres may be detectable on the interior planets by using secondary eclipse MIRI observations to search for the O 3 that is produced by O 2 photolysis.

Broader Planetary Environment & Model Limitations
Future changes to our understanding of initial conditions, model treatment of escape, and the planetary environment post-hydrodynamic loss may influence our results, and we discuss potential ramifications here.Several changes within the model itself may create a more realistic treatment of the escape process, most notable of which would be the addition of a defined atmospheric vertical structure and variable pressure-temperature profile, and the inclusion of non-thermal escape; these beneficial model changes are discussed in Section 4.4.1.In the current study we have explored a broad range of initial water contents at a fixed escape onset time of 1 Myr; however, formation and water delivery research on the TRAPPIST-1 system may provide valuable information that would allow us to narrow the allowed range of initial water contents and conserve computational expense -these considerations are discussed in Section 4.4.2.All of these considerations should be kept in mind when using this work to interpret observations of the system.

Beneficial Model Changes & Additions
Several model updates would make our treatment of the hydrodynamic escape period more realistic, and we find the overwhelming majority of which would likely serve to reduce the cumulative water loss, suggesting the results presented in this study are upper limits to water loss and oxygen accumulation.Since our simulations currently assume atmospheres with perfect mixing and constant pressure and temperature, approximations in escape physics are required as vertical transport and temperature gradients are neglected.Approximations commonly made in escape modeling (Luger & Barnes 2015;Tian 2015;Schaefer et al. 2016;Bourrier et al. 2017) may over-or under-estimate loss rates, particularly in the energy-limited formalism (Krenn et al. 2021) or at boundary conditions between escape regime changes (Volkov et al. 2011;Gronoff et al. 2020).Our description of hydrodynamic escape follows the energy-limited approximation (Hunten 1973;Watson et al. 1981) however, it has been shown that this formalism may considerably deviate from hydrodynamic mass loss models in certain scenarios.Based on Krenn et al. (2021), which assesses the scenarios in which the energy-limited approximation is applicable, all the planets in our study remain in agreement with the more rigorous hydrodynamic loss simulations of Krenn et al. (2021) within a factor of 10 for the entirety of our simulations, except T1-d, and h, and T1-e at the end of their simulations.
Other escape processes not considered in this work may also be acting upon the atmosphere, such as non-thermal escape (e.g., Dong et al. 2018), photochemical escape (e.g., Wordsworth et al. 2018), and core-powered mass-loss (e.g., Ginzburg et al. 2016;Gupta & Schlichting 2019).Though in the case of core-powered mass-loss, this is more typically explored for super-Earth and sub-Neptune populations, and will likely have a greater effect on an initial H/He envelope (Gupta & Schlichting 2019) which we did not take into consideration in this work.
In this work, no photochemistry besides photolysis, and no geological evolution or surface-atmosphere exchange was modeled.In reality, photochemical reactions play a pivotal role in atmospheric escape, particularly in UV shielding, the creation of photochemical products that may act as UV shields, and effects on atmospheric heating and cooling.Through UV shielding, particles at higher altitudes may absorb UV radiation and protect molecules lower in the atmosphere from photodissociation, including self-shielding (e.g., O 2 protecting other O 2 molecules) and shielding between molecular species (Calahan et al. 2022).Of notable importance, the absorption cross sections of water, O 2 , O 3 , and CO 2 overlap heavily (Parkinson & Yoshino 2003;Domagal-Goldman et al. 2014), thus the intense buildup of O 2 we observe in our models may serve to protect water from photodissociation, lowering mass loss rates.Furthermore, the photochemical production of ozone (O 3 ) from free O 2 would be highly active on a planet with an over abundance of free oxygen (Segura et al. 2005;Meadows et al. 2018), increasing the amount of UV shielding molecules and greenhouse gases in the atmosphere.CO 2 , prevalent on all the terrestrial planets with atmospheres in the Solar System and considered a likely constituent of any TRAPPIST-1 planetary atmospheres (Lustig-Yaeger et al. 2019;Lincowski et al. 2018Lincowski et al. , 2022;;Krissansen-Totton & Fortney 2022), is also an important constituent to consider in super heated planetary atmospheres.CO 2 may prolong escape through heating, or bottleneck escape through cooling effects in the middle and upper atmosphere (Kulikov et al. 2007;Wordsworth & Pierrehumbert 2013).
Geological evolution, particularly a magma ocean phase, will be present and active during the hydrodynamic escape regime, with the attendant processes of water trapping in the mantle, outgassing, and the removal of oxygen or other species from the atmosphere (Barth et al. 2021;Krissansen-Totton & Fortney 2022).Hydrodynamic escape of a steam atmosphere is typically thought to be accompanied by a magma ocean phase (Hamano et al. 2013;Acuña et al. 2021), as water provides a significant greenhouse effect that can prolong crustal solidification time.A magma ocean may also prolong periods of escape by slowing the supply rate of water to the atmosphere, effectively reducing cumulative water loss.However, though cumulative water loss would be reduced, the atmosphere itself may be dry for larger initial water contents when including a magma ocean due to water trapping in the mantle (Barth et al. 2021).A magma ocean would also suppress rapid build-up of atmospheric oxygen, prolonging energy-limited escape and increasing planetary oxygen retention by trapping it in the mantle.However, as initial water content increases, atmospheric oxygen buildup may still be inevitable, even with a magma ocean (Barth et al. 2021).With a defined vertical structure (e.g., 1D atmospheric model), surface-atmosphere exchange during a magma ocean phase can be modeled explicitly, which would inform upwards transport of atmospheric constituents.

Initializing & Changing Surface Water Content
As initial water content represents our broadest initial parameter space, planetary formation and late water delivery research could narrow the range of plausible initial conditions used as input, providing tighter constraints on predicted outcomes for planets in the system, and also conserving computational expense, which could enable an increase in model complexity leading to more realistic treatment.Here we conservatively assume formation and water delivery is complete by 1 Myr.However, if a longer formation and water delivery time casued escape to begin later when the incident XUV was lower, then cumulative water loss would be reduced, and the ratio of oxygen retained to water lost would be increased.However, Raymond et al. (2022) show that the growth of the TRAPPIST-1 planets was likely complete within a few Myr with the bulk of their water reservoirs likely accreted during formation, as late water delivery greater than 5% of 1 Earth mass would disrupt the orbital resonance we observe today.In this case, the timeline of our simulations is consistent with these model results, suggesting the planets would have been subjected to the earliest, and most intense, period of hydrodynamic blow-off.
Regarding initial water content, previous work shows the TRAPPIST-1 planets may have formed with or possess up to ∼10 wt% water (Schoonenberg et al. 2019;Agol et al. 2021;Acuña et al. 2021;Raymond et al. 2022), which would provide ∼430 TO on a 1 M planet, roughly consistent with the largest initial water inventories we considered.However, based on the recent T1-b and c observations Greene et al. (2023); Zieba et al. (2023), our work would suggest the planets could not have formed with several hundreds of TO.The lower limit on initial water content, however, is relatively unconstrained by formation studies and extends to 0 for most planets (Raymond et al. 2022), which would remain plausible given the findings of this study, unless T1-c is later confirmed to possess atmospheric oxygen (see §4.3.2).Overall, our timeline of water loss and initial water contents considered are within a plausible range given recent formation and water delivery scenarios, but computational expense may be conserved by exploring a narrowed range of initial water contents.

CONCLUSIONS
We have shown the plausible range of water loss and abiotic oxygen production on the TRAPPIST-1 planets due to hydrodynamic thermal atmospheric escape using an updated VPLanet model (Barnes et al. 2020) and the most recently available constraints on the TRAPPIST-1 system (Agol et al. 2021;Birky et al. 2021).We find the interior planets T1-b, c, and d, are likely desiccated for all but the largest initial water contents (>60, 50, & 30 TO, respectively), and are at the greatest risk of complete atmospheric loss due to their proximity to the host star.However, they may still retain significant oxygen at very large initial water contents depending on assumed sink processes.The outer planets likely retain some fraction of water for all but the lowest initial water contents (∼1 TO).We find T1-e, f, g, and h lose, at most, approximately 8.0, 4.8, 3.4, and 0.8 TO, respectively, corresponding to maximum oxygen production values of approximately 1290, 800, 560, and 90 bars, respectively.Additionally, we have found oxygen retention efficiency may vary widely across distance and initial water content, with rates <10% for close in planets with low initial water up to rates of >80% for distant planets, and those with high initial water.
Using currently available JWST observations for both T1-b and c (Greene et al. 2023;Zieba et al. 2023), our work constrains planetary initial water content, and suggests that the outer planets may still possess water at present-day.In the case of an airless T1-b, and a T1-c with a positive oxygen detection, an initial water content of 8.2 +1.5 −1.0 would be most likely from our study, in agreement with predictions of Zieba et al. (2023).Assuming the outer planets are endowed with comparable initial water, an initial content of 8 TO implies T1-e, f, g, and h all likely possessed ≳1.5 TO after HZ entrance, with the predicted amount increasing with orbital distance.If T1-c is found to be airless, we suggest a basic upper limit on initial surface water of 9 TO and Zieba et al. (2023) suggest an initial water content upper limit of ∼4 TO.Note, an initial water content of 4 TO does not definitively lead to desiccation on the outer planets in our work, which is in agreement with Krissansen-Totton (2023).These conclusions motivate follow-up observations to search for the presence of water vapor or oxygen on T1-c, and future observations of the outer planets in the TRAPPIST-1 system, which may possess substantial water.
Future model upgrades could include implementing a vertical atmospheric structure to enable more detailed escape rate calculations, additional atmospheric species, photochemistry, and surface sinks.As we move into a new era of observation, current work and future studies will assist in defining the plausible environments we could encounter, and statistical studies will be key to retrieving likely evolutionary histories.et al. (2018, 2022).

A.1. Model Validation
To show how the modifications to AtmEsc in this work have changed the outputs of the model from the original version (Luger & Barnes 2015;Lincowski et al. 2018Lincowski et al. , 2022;;Barnes et al. 2020), we present a comparison of one simulation from the original version of AtmEsc to the version modified in this work.Figure 9 shows an escape rate calculation for a single simulation of the TRAPPIST-1 planets with an initial water content of 20 TO from Lincowski et al. (2018Lincowski et al. ( , 2022)), compared to the same calculation done with our model upgrades using the same input parameters.Note, as Lincowski et al. (2018Lincowski et al. ( , 2022) ) used inputs for the TRAPPIST-1 system defined by Grimm et al. (2018), this comparison is meant to show the differences in the models and does not use the updated constraints on the system we adopt in the rest of our work.Additionally, we adopt the same XUV model in this calculation as Lincowski et al. (2018Lincowski et al. ( , 2022)).
There are only minor differences between the two versions of AtmEsc in our comparison.In both models T1-b and c become desiccated, and T1-d is nearly desiccated, losing ∼0.14 TO more water in Lincowski et al. (2018Lincowski et al. ( , 2022)).The three interior planets all produce ∼10 -30% more oxygen in this work; this is because Lincowski et al. (2018Lincowski et al. ( , 2022) ) predicts faster desiccation on these planets, meaning oxygen drag is more efficient during periods of greater water loss preventing build-up.The outer planets all lose more water in this work, typically by ≤0.2 TO, and as a result build up more oxygen, typically by ≤100 bars.
B. OUTLIER REMOVAL Within our full dataset across all initial water contents, a small subset of simulations were identified as outliers originating from the stellar mass samples.As discussed in Section 2.2, we sample stellar and planetary inputs from Birky et al. (2021) and Agol et al. (2021), respectively, creating 20,900 input strings a simulation could draw from.Figure 10 shows the stellar mass samples plotted on the y-axis, with an arbitrarily assigned ID number on the x-axis.The majority of samples fall within 2 standard deviations from the mean, but a small group of samples lie at stellar masses greater than 4 standard deviations from the mean, which we adopt as our outlier cut-off.Outliers account for 259/20900, or 1.24%, of all samples.Simulations using a stellar mass above this cut-off are removed from the final dataset.
These outliers originate from the process used in Birky et al. (2021) to fit the stellar mass of TRAPPIST-1 with the most recent observational data.They used a Gaussian process (GP) surrogate model to estimate the probability of an input parameter, including stellar mass, given the observed stellar properties.Sampling of the GP surrogate posterior was performed using the emcee package (Foreman-Mackey et al. 2013).We found that these outliers originated from .All 20,900 stellar mass samples available for simulations to draw from (from Birky et al. 2021).The y-axis shows the stellar mass in M of a given sample, which are arbitrarily assigned an ID on the x-axis.The red solid line denotes the mean stellar mass across all samples, while the dark blue and subsequent lighter blue regions denote 1, 2, 3, and 4 standard deviations from the mean, as labeled in the plot.The light grey shaded region denotes samples greater than 4 standard deviations from the mean, indicating that samples in this region are classified as outliers, such that simulations using stellar masses above 4 standard deviations are removed from the final data set.
the randomized initial state of the MCMC walkers, suggesting that number of burn-in samples reported in the posterior samples by Birky et al. (2021) was insufficient.The burn-in period was taken to be 2x the average autocorrelation time, which appeared to not reflect the converged state for a subsample of the MCMC walkers.Thus, these outlying stellar mass samples are likely a product of randomized initial walker states, and as they account for 1.24% of all samples, we conclude their removal is justified.

C. PREDICTING INITIAL WATER CONTENT FROM A FINAL STATE
In Section 3.3.1 we showed the most likely initial water content our model would predict to reproduce a known present-day predicted or observed water content on the planet.Here, we detail the generalized method used to achieve this result.
To start, let's assume we have constraints on a planet's present-day water content from either observation or simulation, and models that map initial water content to final water content.For each initial water content in our sampling grid, there are a sufficient amount of simulations (typically ∼3500 -5500 in our work) to encompass the uncertainty in the planet and host star's physical parameters, which in turn define the uncertainty in the final water content at a particular initial content.The question we then pose here is how to determine which initial water content on our sampling grid best produces the final state from our present-day constraints.
The first step is to model the planet's present-day water constraints with a probability density function (PDF) to assign probability weighting to final water contents based on the defined uncertainty in the constraints.In this work, we tested several common statistical functions and identified a gamma profile to best fit the present-day constraint we adopt.Gamma profiles worked well since the water mass fractions are restricted to be greater than or equal to 0 and, consequently, show asymmetric tails, but any continuous probability distribution may be used here in theory.With a gamma profile, our present-day water PDFs are described by: f i (w, a) = w a−1 e −w Γ(a) , (C12) where i denotes which planet the PDF is describing (T1-b, c, d, e, f, g, or h), w is the present-day water content, a is the shape parameter, f i (w, a) is the relative likelihood of present-day water content w on planet i, and Γ(a) is the gamma function evaluated at a particular shape parameter.
With the PDF describing the present-day water constraints, the next step is to determine the value of the PDF at the final water content predicted by each individual simulation across a particular initial water content.We then sum the PDF value across the individual simulations of a given initial water content and normalize for the total number of simulations in that suite to obtain a probability weighting for that particular initial water content.To normalize the probability weightings per initial water content, we divide by the number of simulations in the PDF value summation: , per initial water content, (C13) where f i (w z , a) is the PDF for the planet being considered (Equation C12) evaluated at the final water content of a particular simulation, w z , and N is the total number of simulations at that initial water content.
The initial water content with the highest probability weighting is the one that best produced the assumed presentday water constraints.If the initial water content grid is sampled finely enough, we can interpolate across our distribution of probability weightings to determine where the true most likely initial water content would be, and what the confidence intervals are on that prediction.This next step can only be taken if the grid of initial water contents with simulation data is fine enough that the distribution of probability weightings appears to be a smooth, continuous curve; in some cases this may require running more data once the initial peak probability weighting has been identified.With sufficient data, this interpolation can be accomplished with a Riemann sum.To perform the Riemann sum, we take the curve of probability weightings versus initial water content and create area bins based on the data points, then find at what initial water content is the area under the curve 50% of the total area -the 50 th percentile, or median value (i.e., the initial water content most likely to produce the present-day constraints).Then finding the bounds of 68% and 95% area under the curve, centered at the median value, will yield the 1σ and 2σ confidence intervals.
In our work, we created bins between adjacent data points with linear lines as the top of the bins.Thus, the area in a particular bin between two initial water contents, IW 1 and IW 2, is given by: where s is the slope of the linear fit between the two data points and b is the y-intercept.This integrand is the same equation used to solve for a true value or confidence interval bound that lies between two data points.In fitting for the PDF this work, we used SciPy's stats module (Virtanen et al. 2020), and NumPy's polyfit routine (Harris et al. 2020) for the linear fits.

Figure 1 .
Figure1.Illustration of the possible pathways the TRAPPIST-1 planets can take through escape physics in our model.Interior planets are denoted by dashed curves and habitable zone (HZ) planets are denoted with solid curves.Escape shut off times are shown by vertical dotted lines and times of switches from energy-limited to diffusion-limited escape are shown by vertical dot-dashed lines.Left Column: Escape is halted by complete planetary desiccation, for an interior planet, T1-c (black), and for a HZ planet, T1-f (orange), both with 1 TO initial water inventory.In these cases the atmospheric water is removed entirely, and the planet is last in diffusion-limited escape before the simulation halts.Middle Column: Simulations of HZ or outer planets where escape is halted by entrance to the HZ and the subsequent formation of a cold trap, all for initial water inventories of 10 TO.In this case escape could end in either the diffusion-limited (for T1-e, pink) or energy-limited (for T1-g, teal) regime.Right Column: Escape never shuts off as the planetary water reservoir is sufficiently large to sustain ongoing escape through the systems lifetime.This end state can only occur for interior planets (and results for T-1 b and d with an initial water inventory of 200 TO are shown here) as planets in or beyond the HZ have escape shut off by the formation of a cold trap.Once again, escape could end in either the diffusion-limited (for T1-b, red) or energy-limited (for T1-d, blue) regime, depending on total received insolation and initial water inventory.Top Row: Shows water lost (in TO) over time for each of the column cases previously described.Bottom Row: Shows oxygen produced (in bars) over time for each of the column cases previously described.

Figure 2 .
Figure 2. The predicted final water loss and final oxygen produced and retained as a function of initial water content for the TRAPPIST-1 interior planets T1-b (red), c (black), and d (blue).Top Row: Water lost (in TO) versus initial water (in TO); the shaded blue-grey region indicates where desiccation occurs.For water contents in the upper left panel, all three planets are desiccated and thus their curves are overlapping.Bottom Row: Oxygen produced and retained (in bars) versus initial water.Left Column: Low initial water contents (1 -20 TO).Right Column: High initial water contents (30 -500 TO).Shaded regions give the 1σ uncertainty ranges (set by the 16 th and 84 th percentile values of the final state distributions).

Figure 3 .
Figure3.The predicted final water loss and final oxygen produced and retained as a function of initial water content for the TRAPPIST-1 habitable zone and outer planets T1-e (pink), f (orange), g (teal), and h (purple).The layout of the figure is the same as Figure2.

Figure 4 .
Figure4.The likelihood of each initial water content (in TO) needed to reproduce the predicted present-day water contents fromAgol et al. (2021) for each of the outer planets.The peaks denote the initial water content on our sampling grid most likely to produce the final content in our predictions.An additional 9 simulation suites were ran for T1-e only, between 100 -200 TO to fill out a smooth, continuous curve, and results from these additional simulations are denoted by triangle markers.

Figure 5 .
Figure 5.An illustration of the method used to determine the most probable initial water content for T1-e that would reproduce our predicted present-day content from Agol et al. (2021) (Table3).Same layout as Figure4but with linear x-axis scaling; the y-axis is the likelihood that each initial water content (in TO; x-axis) produces the desired final state.The data points that were simulated are shown by the pink points, the bins used in the Riemann sum are shown in grey, the determined 50 th percentile or true predicted value is shown by the vertical dashed black line, and the 1σ and 2σ confidence intervals are shown by the shaded dark and light blue regions, respectively.The low initial water contents (1 -20 TO) are much more finely spaced (step size of 1 TO) which is why the Riemann sum bin widths are not resolvable there.Lines connecting data points are used for the tops of the Riemann sum bins.

Figure 6 .
Figure6.Left: The oxygen left in the atmospheres of T1-b (red) and T1-c (blue) after an estimated removal of 3 TO worth of oxgyen from a magma ocean, 100 bars from non-thermal escape, and 10 Tmol/yr (∼55 bars/Gyr) from dry crustal oxidation after the planet becomes desiccated (see text for estimation explanation).Right: The resulting histogram of initial water contents predicted by MCMC when fitting for the scenario that T1-b has lost all oxygen while T1-c retains 0.1 bar of atmospheric oxygen following the estimated removal.The vertical solid black lines on both panels denotes the 50 th percentile value predicted by the MCMC while the vertical dashed black lines denote the 68% (1σ) confidence interval.This analysis predicts a minimum initial water content of 8.2 +1.5 −1.0 TO to produce the desired final state, assuming both planets begin with the same initial water content.The amount of oxygen removal is meant to be an underestimate to produce a conservative lower limit, as additional oxygen removal would increase the minimum initial water required.

Figure 7 .
Figure 7. Oxygen retention efficiency as a function of initial water content in TO.Top Row: Interior planets, T1-b, c, and d.Bottom Row: Outer planets, T1-e, f, g,and h.Coloring shows the percentage of simulations at a particular initial water content that fell in a given oxygen retention efficiency bin (e.g., each column corresponds to one initial water content and sums to 100%).Contours are given for the 50 th percentile values (solid white lines) and the 1σ confidence region (dashed white lines).

Figure 8 .
Figure 8.A comparison of the water loss and oxygen production rates from Krissansen-Totton & Fortney (2022) to this study.Design of the figure is meant to recreate the style published in Krissansen-Totton & Fortney (2022); all shaded regions denote 2σ uncertainties.The red, purple, and blue lines and shaded regions show the results from this study for initial water contents of 1, 15, and 100 TO, respectively, and the grey line and shaded region shows the results presented in Krissansen-Totton & Fortney (2022).Left Column: Comparisons for the interior planet TRAPPIST-1 b (see Figure 2 in Krissansen-Totton & Fortney 2022).Right Column: Comparisons for the habitable zone planet TRAPPIST-1 e (see Figure 3 in Krissansen-Totton & Fortney 2022).Top Row: Rate of water loss due to thermal atmospheric escape [Tmol/yr] over time in years.Note, as water is being lost, the flux is negative.Bottom Row: Rates of oxygen production (i.e., gain flux) [Tmol/yr] due to water photolysis and subsequent hydrogen escape over time in years.In general, our results are in agreement with those of Krissansen-Totton & Fortney (2022), and our median test case with an initial water inventory of 15 TO (log space average between 0.7 and 300 TO) stays within their model uncertainties for the entire age of the system.Note, in the T1-e comparisons, our 15 and 100 TO cases display nearly identical evolution and overlap on the plot.

Figure 9 .
Figure9.A comparison to the escape rates of the TRAPPIST-1 planets calculated inLincowski et al. (2018Lincowski et al. ( , 2022) )  using constraints on the system fromGrimm et al. (2018).Left Plot: The water lost in TO over time in years beginning with 20 TO initially.Right Plot: The oxygen produced in bars over time in years.Solid lines denote the rates calculated using the new version of AtmEsc from this work, and dashed lines denote the rates calculated using the version of AtmEsc fromLincowski et al. (2018Lincowski  et al. ( , 2022)).
Figure10.All 20,900 stellar mass samples available for simulations to draw from (fromBirky et al. 2021).The y-axis shows the stellar mass in M of a given sample, which are arbitrarily assigned an ID on the x-axis.The red solid line denotes the mean stellar mass across all samples, while the dark blue and subsequent lighter blue regions denote 1, 2, 3, and 4 standard deviations from the mean, as labeled in the plot.The light grey shaded region denotes samples greater than 4 standard deviations from the mean, indicating that samples in this region are classified as outliers, such that simulations using stellar masses above 4 standard deviations are removed from the final data set.

Table 2 .
Birky et al. (2021)lues for all stellar parameters used as input found from 20900 samples.Uncertainties given are set by 16 th and 84 th percentile values (1σ).Values come fromBirky et al. (2021).
+70.6  −37.2TO, and the 2σ confidence value is 90.6 +155 −71.3 TO.Applications and importance of this method are further discussed in Section 4.3.1.

Table 4 .
The efficiency of oxygen retention [%] after the cumulative water loss predicted by our model at the lowest (1 TO, middle column) and the highest (500 TO, right column) initial water content.Uncertainties reported are 1σ.

Table 5 .
The efficiency of oxygen retention [%] following a 1 TO water loss.Uncertainties reported are 1σ.The middle column shows the results for planets with an initial inventory of 1 TO, and the right column shows the results for planets with an initial inventory of 9, 10, or 11 TO.