ZEPHYRUS model overview

ZEPHYRUS models two channels of atmospheric mass loss for rocky exoplanets coupled to the PROTEUS interior-atmosphere framework: the continuous, bulk hydrodynamic escape driven by stellar XUV irradiation, and the impulsive erosion caused by giant impacts during accretion. The continuous channel implements an energy-limited (EL) formalism following Watson et al. (1981) ¹ and Lopez & Fortney (2013) ²; it is called at each PROTEUS time step with the current planetary radius and mass, the stellar XUV flux supplied by MORS, and the escape radius computed from the atmospheric structure produced by AGNI or JANUS. The mass-loss rate it returns is distributed across atmospheric species according to their elemental mass mixing ratios, so the atmosphere is depleted in bulk without elemental fractionation. The impulsive channel implements the giant-impact erosion scaling laws of Kegerreis et al. (2020) ¹¹ and Roche et al. (2026) ¹², which return the fraction of the target's atmosphere removed by a single collision.

A model parameter reference can be found here.

Energy-limited escape

The mass-loss rate is computed by escape.EL_escape as

where \(\epsilon\) is the escape efficiency factor (epsilon), \(R_\mathrm{XUV}\) is the planetary radius at which the atmosphere becomes optically thick to stellar XUV photons (Rxuv), \(F_\mathrm{XUV}\) is the XUV flux received at the planet (Fxuv) supplied by MORS, \(M_p\) is the planetary mass (Mp), \(G\) is the gravitational constant, and \(K_\mathrm{tide}\) is the tidal correction factor described below. The efficiency \(\epsilon\) quantifies the fraction of incident XUV energy that is converted into work against gravity to drive the outflow; canonical values for rocky planets lie in the range \(0.1 \leq \epsilon \leq 0.3\), although ZEPHYRUS accepts any \(\epsilon \in (0, 1]\).

Radius scaling

The cubic radius factor in the numerator of Eq. (1) is selected at runtime by the scaling argument of escape.EL_escape:

Both forms reduce to \(R_p^3\) when \(R_\mathrm{XUV} = R_p\), which is the conservative lower bound on the mass-loss rate adopted by Luger & Barnes (2015) ⁴ and Moore et al. (2023) ⁵. Allowing \(R_\mathrm{XUV} > R_p\) increases the effective XUV-absorbing area and therefore the escape rate. In PROTEUS, \(R_\mathrm{XUV}\) is recomputed at each time step from the atmospheric pressure–temperature profile at a user-specified reference pressure \(P_\mathrm{XUV}\).

Tidal correction \(K_\mathrm{tide}\)

When the tidal_contribution flag is True, the effective gravitational potential is reduced by the host star's tidal field following the tidal reduction factor of Erkaev et al. (2007), eq. 17 ³:

with the Hill radius

where \(a\) is the planetary semi-major axis, \(e\) is the orbital eccentricity, and \(M_\star\) is the stellar mass. Factoring the numerator gives \(K_\mathrm{tide} = (\xi - 1)^2\,(2\xi + 1) / (2\xi^3)\), which is non-negative for every \(\xi > 0\) with a double root at \(\xi = 1\). In the physical regime \(\xi > 1\) it lies in \((0, 1)\), rising toward 1 for \(\xi \gg 1\) (the XUV radius well inside the Hill sphere) and falling toward 0 as the atmosphere expands toward the Roche lobe at \(\xi = 1\); because the escape rate divides by \(K_\mathrm{tide}\), the rate is enhanced by the tidal correction and diverges as \(\xi \to 1\). The tidally corrected rate is therefore defined only for \(\xi > 1\): ZEPHYRUS raises a ValueError for \(\xi \le 1\), where the atmosphere reaches the Roche lobe and the energy-limited approximation no longer applies. When tidal_contribution is False, \(K_\mathrm{tide} = 1\) is enforced.

Giant-impact atmospheric erosion

A giant impact removes part of the target planet's atmosphere in a single event. ZEPHYRUS computes the eroded fraction with collision.mass_loss, which implements the scaling law of Kegerreis et al. (2020), their Eq. 1 ¹¹:

capped at 1 for total erosion, where subscript \(i\) denotes the impactor, \(t\) the target, \(M_\mathrm{tot} = M_i + M_t\), and \(b \equiv \sin\beta\) is the dimensionless impact parameter for impact angle \(\beta\) (0 head-on, 1 fully grazing). The prefactor and exponent are least-squares fits to the paper's suite of 259 SPH simulations, each with an uncertainty of 0.01. The mutual escape speed of the pair at contact is

and \(f_M(b)\) is the fractional interacting mass of the pair (their Eq. B1), built from density-weighted spherical caps of common height \(d = (R_t + R_i)(1 - b)\):

where \(V_{t,i}\) are the full body volumes. At equal bulk densities \(f_M\) reduces exactly to the fractional interacting volume of their Eq. B2. The common-height caps are a linearised bookkeeping: outside the fitted geometry, for a much denser and much smaller impactor near head-on, the raw \(f_M\) can leave \([0, 1]\) and vary non-monotonically with \(b\), so ZEPHYRUS clamps \(f_M\) to \([0, 1]\). Within the fitted domain the clamp never engages.

Three input conventions follow the paper and must be honoured by the caller: \(v_c\) is the speed at first contact, not the relative speed at infinity; the masses and radii exclude any atmosphere, with radii taken at the base of the atmosphere; and the densities are bulk values of the atmosphere-free bodies.

The returned fraction applies to the target's atmosphere as a whole. Consistent with the bulk-removal treatment of the continuous channel, the caller partitions the lost mass across atmospheric species without elemental fractionation.

The Roche et al. (2026) law

Roche et al. (2026) ¹² generalize atmospheric erosion scaling to account explicitly for envelope mass fraction \(f_\mathrm{atm} \equiv M_\mathrm{atm} / M_\mathrm{t}^\mathrm{tot}\). The scaling law splits atmospheric loss into near-field erosion (\(X_\mathrm{NF}\), loss near the impact site) and far-field erosion (\(X_\mathrm{FF}\), loss through ground motion):

where \(X_\mathrm{atm}\) is clamped to \([0, 1]\). The fraction of the atmosphere located in the near-field region, \(f_\mathrm{NF}\), is parameterized by a generalized logistic function:

clamped to \([0, 1]\), and the far-field envelope fraction is \(f_\mathrm{FF} = 1 - f_\mathrm{NF}\). Here \(R_\mathrm{t}^\mathrm{r}\) and \(R_\mathrm{i}^\mathrm{r}\) are the refractory (atmosphere-free) core-plus-mantle radii.

Near-field loss \(X_\mathrm{NF}\) is described by:

with a velocity floor at \(v_\mathrm{c} / v_\mathrm{esc} = 1.0\), clamped to \([\max(0, X_\mathrm{NF}(v_\mathrm{c} = v_\mathrm{esc})), f_\mathrm{NF}]\).

Each coefficient vector (\(\boldsymbol{\zeta}\), \(\boldsymbol{\xi}\), \(\boldsymbol{\psi}\)) is evaluated as an empirical function combining power-law, polynomial, and logarithmic dependencies on impact parameter \(b\), envelope mass fraction \(f_\mathrm{atm}\), refractory target mass \(M_\mathrm{t}^\mathrm{r} / M_\oplus\), and refractory impactor mass fraction \(\gamma \equiv M_\mathrm{i}^\mathrm{r} / (M_\mathrm{i}^\mathrm{r} + M_\mathrm{t}^\mathrm{r})\):

The far-field loss \(X_\mathrm{FF}\) accounts for ground motion and scales with the modified specific impact energy \(Q'_\mathrm{R}\):

with \(Q'_\mathrm{R}\) expressed in \(\mathrm{MJ\,kg^{-1}}\). Part of the far-field loss does not depend on impact energy whenever \(\psi_1 + \psi_3 > 0\), driven by the \(s_{16} \log_{10}(f_\mathrm{atm})\) term of \(\psi_1\). At zero impact energy the empirical fit yields \(X_\mathrm{FF} = f_\mathrm{FF}\operatorname{clip}(\psi_1 + \psi_3, 0, 1)\), which is strictly positive when \(\psi_1 + \psi_3 > 0\). Within the fitted parameter box, \(\psi_1 + \psi_3 > 0\) occurs above threshold values of \(\gamma\) that depend on target mass:

Below \(f_\mathrm{atm} = 0.01\) this positive zero-energy region expands. Consequently, the 'X_FF_zero_energy' diagnostic flag can be raised for scenarios within the fitted range, as in two of the authors' reference impacts: one with \(\gamma = 0.5\) at \(f_\mathrm{atm} = 0.01\) on a \(1\,M_\oplus\) target and one with \(\gamma = 0.4\) at \(f_\mathrm{atm} = 0.05\) on a \(0.35\,M_\oplus\) target. ZEPHYRUS evaluates the authors' fit and reports this zero-energy far-field contribution in diagnostics['X_FF_zero_energy'] and the 'X_FF_zero_energy' flag. For the authors' own Moon-forming impact scenarios (Table D1 of Roche et al. 2026, evaluated with radii scaling as \(M^{1/4}\,R_\oplus\)), the law gives the following loss fractions at \(f_\mathrm{atm} = 10^{-4}\): CA01 0.208, R12 0.310, CS12 0.467, C12 0.466, LS18a 0.613, and LS18b 0.363. For \(f_\mathrm{atm}\) from \(10^{-4}\) to \(10^{-6}\) the law gives CA01 0.208 to 0.296 and LS18a 0.613 to 0.764; Roche et al. (2026, Sect. 4.3) estimate about 20% to 30% for the canonical impact (CA01) and about 70% to 80% for the most energetic synestia-forming scenario (Lock et al. 2018) for an atmosphere of 100 bar or less. The law reaches 0.70 for LS18a between \(f_\mathrm{atm} = 10^{-5}\) and \(10^{-6}\) (0.688, 0.764), yielding 0.613 at \(10^{-4}\). For CA01 at \(f_\mathrm{atm} = 10^{-4}\), the zero-energy far-field part is 0.057 of 0.208 (27.5% of the loss, 75% of the far-field part). Over the 790 fitting rows, the mean and maximum absolute misfits against simulation data are 0.0396 and 0.2224.

The forms on this page follow the authors' published code, which reproduces their Fig. 3. The printed Roche et al. (2026), Eq. 6 has \(\zeta_3 R - \zeta_2\) in the exponent, the printed Roche et al. (2026), Eq. 10 has \(M_\mathrm{t}^\mathrm{tot}\) in the mass ratio, and the \(X_\mathrm{NF}\) floor at \(v_\mathrm{c} = v_\mathrm{esc}\) is in the code and not in the paper; Table C prints some values with only 2 significant digits (\(s_{15}, k_{15}, s_{25}, s_{33}\)), and the constant term \(q_{41}\) of \(\zeta_4\) (161.057) and \(\zeta_5\) (-161.022) almost cancel, so values from the published dataset are needed.

Impact energy and escape speed definitions

The mutual escape speed at contact uses the total target mass and refractory impactor mass with refractory radii:

The modified specific impact energy \(Q'_\mathrm{R}\) follows the interacting-mass formulation of Leinhardt & Stewart (2012) ¹⁴ and Roche et al. (2025) ¹³, Eqns. 13 to 18:

where \(\mu \equiv M_\mathrm{i}^\mathrm{r} M_\mathrm{t}^\mathrm{tot} / (M_\mathrm{i}^\mathrm{r} + M_\mathrm{t}^\mathrm{tot})\) and \(\mu_\alpha \equiv \alpha M_\mathrm{i}^\mathrm{r} M_\mathrm{t}^\mathrm{tot} / (\alpha M_\mathrm{i}^\mathrm{r} + M_\mathrm{t}^\mathrm{tot})\). Here \(\alpha\) is the interacting mass fraction of the impactor:

with \(l \equiv R_\mathrm{t}^\mathrm{r} + R_\mathrm{i}^\mathrm{r} - (R_\mathrm{t}^\mathrm{r} + R_\mathrm{i}^\mathrm{r}) b\) if \((R_\mathrm{t}^\mathrm{r} + R_\mathrm{i}^\mathrm{r}) b + R_\mathrm{i}^\mathrm{r} > R_\mathrm{t}^\mathrm{r}\), and \(\alpha = 1\) when the impactor is completely intercepted.

Symbols and units

Coefficient origin

The scaling coefficients comprise 61 stored values (Roche et al. 2026, Appendix C):

• The 23 near-field mass coefficients (\(q_{ij}\), \(\zeta_5\), \(\zeta_6\)) are fitted to a separate suite of initialised SPH planets without impacts (radius ratios 0.1 to 1.0 in steps of 0.1 plus 0.001, \(M_\mathrm{t}^\mathrm{r} \in [0.01, 5.0]\,M_\oplus\), \(f_\mathrm{atm} \in [0.01, 0.2]\)); \(q_{14} = q_{34} = 2\) are fixed exponents. The impact dataset reaches a radius ratio of 1.015.
• The 38 loss coefficients (\(k_{ij}\), \(s_{ij}\)) are fitted to 296 of the 300 new impact simulations (96 at \(f_\mathrm{atm} = 0.01\), 100 at 0.1, and 100 at 0.2, all on a single target mass of about \(1\,M_\oplus\)) and the 494 simulations at \(f_\mathrm{atm} = 0.05\) of Roche et al. (2025) ¹³ (spanning multiple target masses), giving 790 rows total. The loss coefficients (\(k, s\)) were fitted on \(R_\mathrm{i}^\mathrm{r} / R_\mathrm{t}^\mathrm{r} \ge 0.498\) only; the 0.001 lower bound of the radius-ratio range comes from the separate \(f_\mathrm{NF}\) planet-initialisation suite. The velocity exponent is shared across near-field and far-field regimes (\(k_{15} = k_{25} = k_{35}\)). Away from \(1\,M_\oplus\), the impact data cover only \(f_\mathrm{atm} = 0.05\).

All 61 values are transcribed at full precision from the published dataset ¹². In Table C of Roche et al. (2026), some coefficients are printed with only 2 significant digits (\(s_{15}, k_{15}, s_{25}, s_{33}\)), making the dataset values necessary for reproduction. Terms fixed in Tables C1 to C3 (\(q_{14} = q_{34} = 2\), \(q_{37} = q_{48} = 1\), and unlisted parameters fixed to 0) are held in internal constants.

Calibrated domain and stability policy

The simulation suite constrains the law over the following parameter space:

• Target refractory mass: \(M_\mathrm{t}^\mathrm{r} \in [0.35, 5.0]\,M_\oplus\)
• Atmosphere mass fraction: \(f_\mathrm{atm} \in [0.01, 0.2]\)
• Impactor mass fraction: \(\gamma \in [0.1, 0.5]\)
• Impact parameter: \(b \in [0.0, 0.9]\)
• Radius ratio: \(R_\mathrm{i}^\mathrm{r} / R_\mathrm{t}^\mathrm{r} \in [0.001, 1.015]\) (from the 0.001 lower bound of the initialised suite to the 1.015 maximum of the impact data)
• Contact velocity: \(v_\mathrm{c} \in [1.0, 3.0]\,v_\mathrm{esc}\)

To support planetary evolution and accretion calculations where conditions cross these empirical boundaries, zephyrus.collision.impact_loss applies a structured evaluation policy:

1. Diagnostic range flags. When inputs fall outside the calibrated range with a 1% relative tolerance (\(f_\mathrm{atm} \notin [0.01, 0.2]\), \(M_\mathrm{t}^\mathrm{r} \notin [0.35, 5.0]\,M_\oplus\), \(\gamma \notin [0.1, 0.5]\), \(b > 0.9\), \(R_\mathrm{i}^\mathrm{r} / R_\mathrm{t}^\mathrm{r} \notin [0.001, 1.015]\), or \(v_\mathrm{c} / v_\mathrm{esc} > 3.0\)), the loss fraction is computed and the out-of-range parameter names are recorded in ImpactLossResult.flags ('f_atm', 'M_t_earth', 'gamma', 'b', 'R_ratio', 'v_ratio'). Where a stability bound equals a fitted bound (such as \(\gamma = 0.5\)), the 1% relative tolerance does not apply: any clamped value is flagged. Additional diagnostic flags record physical regimes: 'v_sub_escape' triggers when contact speed falls below 0.99 mutual escape speed (\(v_\mathrm{c} < 0.99\,v_\mathrm{esc}\), where the near-field velocity floor evaluates at \(v_\mathrm{esc}\)), and 'X_FF_zero_energy' triggers when the zero-energy far-field loss exceeds zero (\(X_\mathrm{FF,zero} > 0\), which can occur inside the fitted range). In PROTEUS the atmosphere fraction is typically \(10^{-5}\) to \(10^{-3}\), 1 to 3 decades below the calibrated 0.01, so 'f_atm' is flagged. Over a 400-impact grid in the fitted box (\(M_\mathrm{t}^\mathrm{r} \in \{0.35, 1, 2, 5\}\,M_\oplus\), \(\gamma \in \{0.1, 0.2, 0.3, 0.4, 0.5\}\), \(b \in \{0, 0.2, 0.4, 0.6, 0.8\}\), \(v_\mathrm{c} / v_\mathrm{esc} \in \{1, 1.5, 2, 3\}\), radii from equal bulk density), \(X_\mathrm{NF}\) reaches \(f_\mathrm{NF}\) in 172/400 impacts at \(f_\mathrm{atm} = 10^{-2}\), 300/400 at \(10^{-3}\), and 400/400 at \(10^{-4}\), while \(X\) keeps rising through the far-field term (for CA01, \(X\) evaluates to 0.056 at \(10^{-2}\), 0.131 at \(10^{-3}\), 0.208 at \(10^{-4}\), 0.252 at \(10^{-5}\), and 0.296 at \(10^{-6}\)). An airless target (\(f_\mathrm{atm} = 0\)) returns no flags.
2. Stability clamping. The bounds match the parameter range where Roche et al. (2026, Sect. 4.1) verified numerical stability (\(f_\mathrm{atm} \in [10^{-6}, 0.4]\), \(M_\mathrm{t}^\mathrm{r} \in [10^{-3}, 10]\,M_\oplus\), \(\gamma \in [10^{-3}, 0.5]\)). Clamped parameter values enter only the empirical fit arguments (\(\boldsymbol{\zeta}, \boldsymbol{\xi}, \boldsymbol{\psi}\)); physical quantities (\(v_\mathrm{esc}\), \(v_\mathrm{c} / v_\mathrm{esc}\), \(Q'_\mathrm{R}\)) and the far-field mass ratio use the physical input masses and radii. diagnostics['gamma'] reports the raw impactor mass fraction. Applied bounds are recorded in diagnostics['clamped'].
3. Zero atmosphere. When \(f_\mathrm{atm} = 0.0\), the function returns \(X = 0.0\) immediately.
4. Grazing collisions. The scaling law is calibrated for impact parameters \(b \in [0.0, 0.9]\). At \(b = 1\) the geometry factor is 0, so the bracket of Eq. (13) is \(\psi_1 + \psi_3\) and \(X_\mathrm{FF}\) evaluates to the zero-energy value \(f_\mathrm{FF}\operatorname{clip}(\psi_1 + \psi_3, 0, 1)\) (0 in the benchmark case \(M_\mathrm{t}^\mathrm{r} = 1.0\,M_\oplus\), \(\gamma = 0.3\), \(f_\mathrm{atm} = 0.01\), \(v_\mathrm{c} / v_\mathrm{esc} = 1.5\), and 0.080 at \(f_\mathrm{atm} = 10^{-4}\) with \(R_\mathrm{i}^\mathrm{r} / R_\mathrm{t}^\mathrm{r} = (\gamma / (1 - \gamma))^{1/3}\)); \(X\) at \(b = 1\) is \(X_\mathrm{NF} + X_\mathrm{FF,zero}\) (0.0315 for the benchmark case), while collisions with \(b > 0.9\) trigger the 'b' diagnostic flag.

Physical caveats

The Roche et al. (2026, Sect. 4.1) scaling law does not account for pre-impact planetary rotation, surface liquid water oceans, thermal evolution, or a core mass fraction other than about 0.3 (Sect. 2.1). All simulations assume H2-He envelopes governed by the Hubbard & MacFarlane 1980 ¹⁵ equation of state.

Physical effects operate in both directions:

• Miscibility (overestimate): Magma-envelope miscibility at high pressures without a sharp boundary can cause the scaling law to overestimate the loss of massive envelopes on young planets (Roche et al. 2026, Sect. 4.1).
• Post-impact thermal loss (underestimate): The law covers the immediate shock- and vapour-plume-driven loss only; it neglects later thermally driven loss (an outflow driven by heat from the post-impact interior, Biersteker & Schlichting 2021), so for primordial H2-He envelopes the total loss can be higher (Roche et al. 2026, Sect. 4.1); that later loss becomes negligible for envelopes of higher mean molecular weight.
• Atmospheric composition (upper limit): Heavier atmospheres (\(\mathrm{CO}, \mathrm{CO}_2\)) are less susceptible to shock-driven removal, so shock-driven loss is an upper limit for a given atmosphere mass (Roche et al. 2026, Sect. 4.3). For PROTEUS atmospheres of higher mean molecular weight (\(\mathrm{H}_2\mathrm{O}, \mathrm{CO}_2, \mathrm{O}_2\)), the H2-He loss fractions likewise serve as upper limits.

Law selection with impact_loss

Callers select the desired erosion law with zephyrus.collision.impact_loss(law=...), where law is 'kegerreis2020' or 'roche2026'. The function returns an ImpactLossResult dataclass with fraction, flags, and diagnostics.

Coupling to PROTEUS

ZEPHYRUS treats atmospheric escape as a bulk process: at each PROTEUS time step the total mass-loss rate from Eq. (1) is partitioned across atmospheric species in proportion to their elemental mass mixing ratios as computed by CALLIOPE. No elemental fractionation between light and heavy species is imposed in the outflow itself; however, because only outgassed volatiles are subject to escape while dissolved species remain in the magma ocean reservoir, escape fractionates the planet's total (interior + atmosphere) volatile budget over time, preferentially retaining species that are highly soluble in silicate melts (e.g. H\(_2\)O, S\(_2\)).

More about this in its dedicated page.

Regime of validity

The EL formalism is appropriate in the high-irradiation, hydrodynamic regime that dominates atmospheric loss during the first \(\sim 10^6\)–\(10^8\) yr of evolution for close-in rocky planets ¹⁶. Outside this regime—at lower XUV fluxes or for less extended atmospheres—non-thermal escape (Jeans escape, ion pickup, charge exchange) becomes comparable to or exceeds the hydrodynamic rate, and the bulk EL prescription no longer applies. ZEPHYRUS does not currently include these processes; users should verify that the integrated XUV-driven loss exceeds non-thermal estimates (e.g. \(\sim 10^7\)–\(10^8\) g s\(^{-1}\) for an Earth-mass planet; Kislyakova et al. 2014 ⁷) before interpreting model outputs.

Similarly, the bulk-removal assumption breaks down when the hydrodynamic particle flux drops below the critical flux required to drag heavy species against gravity, at which point compositional fractionation in the outflow becomes significant ⁸⁹. Following Yoshida et al. (2022) ¹⁰, the critical flux for H\(_2\)O in an H\(_2\) background is \(\approx 1.9 \times 10^{8}\) g s\(^{-1}\).

The Kegerreis et al. (2020) giant-impact erosion law (Eq. 4) is constrained by simulations spanning target masses of roughly 0.3 to 3 \(M_\oplus\), impactor masses down to about 0.05 \(M_\oplus\), bulk densities from about half to double Earth's, contact speeds of 1 to 3 \(v_\mathrm{esc}\), all impact angles, and thin atmospheres of order 1 percent of the planet mass. The median deviation of the simulations from the law is 9 percent, rising to about 20 percent for slow, head-on impacts, whose outcomes are chaotic. For the Kegerreis et al. (2020) law, the loss depends only mildly on the atmosphere mass in this thin-atmosphere regime, with a factor of 10 less atmosphere increasing the eroded fraction by roughly 10 percent (Kegerreis et al. 2020); substantially thicker atmospheres, which can cushion the impactor, fall outside the law's regime. In contrast, under the Roche et al. (2026) law, atmospheric loss fraction exhibits a stronger dependence on envelope mass: for the canonical Moon-forming impact (CA01), \(X\) increases from 0.056 to 0.131 as \(f_\mathrm{atm}\) decreases from \(10^{-2}\) to \(10^{-3}\).