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 law of Kegerreis et al. (2020) ¹¹, which returns 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.

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 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. 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; substantially thicker atmospheres, which can cushion the impactor, fall outside the law's regime.