Atmospheric Loss in Giant Impacts Depends on Preimpact Surface Conditions

To calculate atmospheric and ocean loss due to ground motion, we follow a similar approach to that of Genda & Abe (2003, 2005). A hydrostatic atmosphere, and in most cases a hydrostatic ocean, is initialized at the radius of the planet's surface in a 1D spherical geometry (Figure 1(D)). The breakout of the shock wave is then simulated by giving the lower boundary a vertical velocity that generates a shock wave in the (ocean and) atmosphere that accelerates a fraction of the (ocean and) atmosphere to escape.

We adapted the 1D WONDY hydrodynamic code (Kipp & Lawrence 1982) to calculate the evolution of the atmosphere (and ocean) in response to the motion of the ground. WONDY solves the Lagrangian 1D mass, momentum, and energy equations using a finite-difference method. Artificial viscosity is used to resolve shocks by spreading the shock front over several Lagrangian cells. We have expanded the capabilities of the WONDY code by adding options for radial gravity, three additional equations of state (EOSs), and hydrostatic initialization of an atmosphere and ocean. A full description of the adapted code and sensitivity tests are given in Appendix A, but we provide a summary here.

The atmosphere and ocean were both modeled using 500 Lagrangian cells (or zones) each, and are initialized as stationary, hydrostatic, and adiabatic, isothermal or isoenergetic depending on the EOS used. By assuming that the atmosphere and ocean are stationary and hydrostatic we neglect the influence of the gravity of the other colliding body, which could deform the surface and disturb the atmospheric structure. We will address the effect of the preimpact redistribution of the atmosphere and ocean in future work. The properties of the atmosphere were set by defining a surface temperature (T₀) and pressure (p₀) and the structure of the atmosphere was determined by integrating upwards from the surface. The top of the atmosphere is treated as a stress-free boundary, implemented by using an additional mass-less, zero-pressure cell at the top of the atmosphere. The atmosphere was modeled as an ideal gas with a constant molar mass (m_a) and ratio of specific heat capacities (γ). When present, the ocean was initialized with a given depth (H_oc) and the initial structure of the ocean was calculated by integrating downwards from the ocean surface, assuming thermal equilibrium with the atmosphere at the base of the atmosphere. In most simulations, we used the water EOS of Senft & Stewart (2008) to describe the thermodynamic properties of the ocean. In order to compare our results to those of Genda & Abe (2005), we also ran simulations using the International Association for the Properties of Water and Steam (IAPWS) EOS (Wagner 2002) and the Tilloston EOS (Tillotson 1962) using the parameters from O'Keefe & Ahrens (1982). We find good agreement between our results using the Senft & Stewart tabulated EOS and the IAPWS EOS, which is not surprising as the Senft & Stewart EOS was constructed partly using the IAPWS EOS, but find significant differences when using the Tillotson EOS, which we discuss in Section 4.2.1.

As in previous work (Genda & Abe 2003, 2005), we do not directly model the shock in the planet. Instead, the propagation of the shock from the planet into the ocean/atmosphere is simulated by imposing the velocity of the lower boundary of the domain, i.e., the ground motion. The boundary is given an initial velocity (u_G) and then allowed to follow a ballistic trajectory, ignoring the influence of any forces other than gravity (see Section A.3 for more details). Imposing a ballistic boundary condition assumes that the mass of the ground is much greater than the mass of any ocean and/or atmosphere and thus the ground is not slowed significantly by transferring momentum to the ocean and/or atmosphere. For the range of oceans and atmospheres we consider in this work, this is a good approximation. Even for the most massive atmosphere and ocean combination we simulated, the mass of the ocean and atmosphere combined is equivalent to a surface layer of only ∼10 km, and is typically much lower. Using the ballistic boundary condition, if the ground velocity is below the escape velocity, the boundary eventually stops and then accelerates downwards toward its initial position. When the boundary approaches its initial position, it is gradually brought to a stop. The prescription of this later-stage evolution of the boundary rarely has an effect on the amount of loss. We discuss the effect of nonballistic motion of the boundary in Section 6. It is important to note that u_G is the velocity of the ground upon breakout of the shock into the atmosphere or ocean (i.e., the impedance-match velocity; see Section 2), and is not the particle velocity of the shock in the planet. The relationship between the strength of the shock in the planet and u_G can vary depending on the properties of the atmosphere or ocean, and we discuss this in Section 5. Furthermore, prescribing a ballistic trajectory ignores any further positive acceleration of the ground by decompression of the surface to pressures below the impedance-match solution. We discuss the implications of this simplification in Section 6.

Simulations were run for 5000 s to ensure that a plateau in atmospheric/ocean loss was achieved. A small number of runs failed before completion. Failed runs were typically for either particularly high or low ground velocities. Failure was generally due to either insufficient numerical viscosity in the first few time steps or due to complications with stopping of the boundary upon its ballistic descent late in time. If a plateau in loss had been reached prior to failure this value was taken as the final loss, otherwise the result was discounted and not considered in our results. In another subset of cases, the stopping of the boundary upon descent caused a secondary shock into the atmosphere, leading to additional loss. This secondary shock is likely unrealistic as the surface would have spalled or vaporized shortly after release. Generally, the amount of additional loss is small as the initial hydrostatic structure of the atmosphere that allowed for the acceleration of the initial shock has been disrupted. To account for this effect, we identified the plateau in loss due to the original shock and took the loss just before the passage of the second wave as the final value for loss. We tested the sensitivity of our results to the intrinsic parameters of the code and found no variation within the range of reasonable values (Appendix B.1).

We ran simulations for a wide variety of surface pressures, surface temperatures, atmospheric compositions, ocean depths, planetary masses, ground velocities, and using the three different EOS for water to explore the dependence of each of these parameters on the efficiency of loss. The details of the surface conditions used in each set of calculations are described at the relevant point in Section 4.

The impedance-match velocities and pressures between different layers were found numerically by solving for the intersection between the relevant release curve (of the ground/ocean) with the Hugoniot of the layer above. Hugoniots were calculated by iteratively finding the particle velocity, u_p, that satisfied the Rankine–Hugoniot equations, or using the analytical expressions for shocks in an ideal gas (Melosh 1989; Zel'dovich & Raizer 2002). Release curves were calculated as isentropes through the relevant EOS with the solutions found iteratively, if necessary. Jupyter notebooks, Python scripts, and a widget that can calculate the impedance match between different materials will be made available on publication of this work.

The EOSs used for water and atmospheres were the same as for the 1D hydrodynamics simulations (Section 3.1), and the ground was modeled as forsterite (Stewart et al. 2019). To allow discussion of previous work (Kegerreis et al. 2020a, 2020b), we also calculated impedance matches using the EOS for H₂–He mixtures from Hubbard & MacFarlane (1980). An early version of the Hubbard & MacFarlane (1980) EOS that was included in the SWIFT hydrodynamics code (Schaller et al. 2018; Kegerreis et al. 2019) and used in previous work (Kegerreis et al. 2018, 2020a, 2020b) contained an error in the calculation of internal energy. The Hubbard & MacFarlane (1980) EOS is defined by expressions for pressure and heat capacity as functions of density and temperature. Previous work calculated the specific internal energy as

$$\begin{eqnarray}&&\epsilon ={c}_{V}(\rho ,T)\,T\ ,\end{eqnarray} \tag{ 1 }$$

where c_V is the specific heat capacity, ρ is the density, and T is temperature, but this expression neglects the change in heat capacity as a function of temperature and density. Here, we calculate internal energy as an integration first along the T = 0 K isotherm, and then along an isochore:

$$\begin{eqnarray}&&\epsilon (\rho ,T)={\int }_{{\rho }_{0}}^{\rho }{\left.\displaystyle \frac{\partial \epsilon }{\partial \rho ^{\prime} }\right|}_{T=0}d\rho ^{\prime} +{\int }_{0}^{T}{\left.\displaystyle \frac{\partial \epsilon }{\partial T^{\prime} }\right|}_{\rho ^{\prime} =\rho }{dT}^{\prime} ,\end{eqnarray} \tag{ 2 }$$

where primes indicate integration variables. In the formulation of Hubbard & MacFarlane (1980) the first term, the integral along the isotherm, is zero and so

$$\begin{eqnarray}&&\epsilon (\rho ,T)={\int }_{0}^{T}{\left.\displaystyle \frac{\partial \epsilon }{\partial T^{\prime} }\right|}_{\rho ^{\prime} =\rho }{dT}^{\prime} ={\int }_{0}^{T}{c}_{V}(\rho ,T^{\prime} ){dT}^{\prime} ,\end{eqnarray} \tag{ 3 }$$

which can be solved analytically. This is now the method used for calculating internal energy in the Hubbard & MacFarlane (1980) EOS table included in the current and future releases of SWIFT. It is important to note that an EOS defined by expressions for heat capacity and pressure alone, such as the Hubbard & MacFarlane (1980) EOS, is nonconservative. In other words, integration of thermodynamics variable along different paths between points in phase space can give different values. There is thus no single definition of internal energy, and our choice of energy calculation only improves on previous work in that it gives a value that is consistent with the defined EOS.

Atmospheric loss in the no-ocean case has been considered in a number of previous studies (e.g., Chen & Ahrens 1997; Genda & Abe 2003; Schlichting et al. 2015). Genda & Abe (2003) explored seven cases of ideal, adiabatic atmospheres with different atmospheric compositions, surface temperatures and pressures, and atmospheric compositions, and conducted simulations using different prescriptions for γ. They found that the degree of loss is relatively minimal unless the ground velocity exceeds ∼0.5 v_esc and observed relatively little variation in the efficiency of loss between different atmospheric properties and prescriptions for γ. Schlichting et al. (2015) calculated loss for both isothermal and adiabatic atmospheres with γ = 4/3 and γ = 5/3. In agreement with Genda & Abe (2003), they found little difference in the efficiency of loss between atmospheres with different γ, but found that loss from isothermal atmospheres was somewhat less efficient than for adiabatic atmospheres at the same ground velocity (a difference in loss of up to ∼10%). Here, we revisit atmospheric loss in the absence of an ocean to ground truth our numerical model and to further explore the effect of atmospheric properties and planetary mass on the efficiency of loss.

Figure 4 shows the evolution of an atmosphere upon breakout of a shock on a planet's surface as calculated using our 1D hydrocode. The evolution follows that expected based on the physics described in Section 2. The pressure of the shock at the base of the atmosphere is set by the impedance-match solution at the particle velocity of the ground imposed in the simulation. The shock wave accelerates up the strong adiabatic density gradient of the atmosphere, heating and compressing the gas, until the shock front reaches the top of the atmosphere (at ∼3.2 s in Figure 4). When the shock reaches the low-density edge of the atmosphere the compressed gas expands rapidly, reaching speeds far in excess of the escape velocity, and the top of the atmosphere is lost from the gravitational well of the planet (Figures 4(E), (I)). Momentum transfer to the portion of the atmosphere that is lost occurs in the first few seconds to tens of seconds after the release of the shock into the atmosphere. After this point the lost portion of the atmosphere behaves almost ballistically and its velocity begins to fall as it moves further from the planet. In our simulations, the ground eventually stops (at 830 s in Figure 4) and the remaining bound atmosphere begins to fall back down to the planet. In reality, before this point the release wave from other parts of the surface could have slowed the ground motion and the rock surface could have spalled, melted, or vaporized. However, as the momentum is transferred to the lost portion of the atmosphere very early, these complications likely do not affect the efficiency of loss from the initial shock wave (see Section 6 and Kegerreis et al. 2019).

We find good agreement between our results and those of previous studies, and demonstrate more completely that the relationship between ground velocity and atmospheric loss in the absence of an ocean is relatively insensitive to atmospheric composition, surface temperature and pressure, and planetary mass. Figure 5(A) shows the efficiency of atmospheric loss for H₂ (m_a = 2 g mol⁻¹, γ = 1.4), H₂O (m_a = 18 g mol⁻¹, γ = 1.25), CO₂ (m_a = 44 g mol⁻¹, γ = 1.29) and an approximation to an Earth-like N₂- and O₂-dominated atmosphere (m_a = 29 g mol⁻¹, γ = 1.4; Genda & Abe 2003) with surface temperatures of 283 and 3000 K (H₂, CO₂, and Earth-like) or 300 (H₂O) and surface pressures of 100 bar. The black dashed line is the result from Genda & Abe (2003) for a 1 bar atmosphere of Earth-like composition and a surface temperature of 288 K. With the exception of the 3000 K, H₂ atmosphere, all of the simulations gave very similar results and were in good agreement with those of Genda & Abe (2003) and Schlichting et al. (2015). Loss of the high-temperature H₂ atmosphere is slightly less efficient for a given ground velocity (at most a few percent), which may be surprising given that the high-T H₂ atmosphere is much more extended, and so more loosely bound, than the other example atmospheres. The high-T H₂ atmosphere has a height of 2.3 R_Earth (Earth radii), compared to a maximum height of 0.07 R_Earth for the other examples. The pressure and density gradients are much lower, affecting how the shock wave accelerates through the atmosphere, and more of the mass of the atmosphere is at lower pressure. In addition, the hot H₂ atmosphere is so extended that for ground velocities less than ∼0.6 v_esc the ground has stopped and begun to fall back to its original postion even before the initial shock has reached the top of the atmosphere. The release wave from the reversal of the ground velocity propagates upwards in the atmosphere and may play a role in reducing the efficiency of loss compared to other cases where the shock wave is supported as the atmosphere is accelerated to escape.

Zoom In Zoom Out Reset image size

Figure 5(B) shows the efficiency of atmospheric loss for atmospheres of varying surface pressures on planets of between Mars and Earth mass. Points are for all combinations of atmospheric pressures of 0.1, 0.5, 1, 5, 10, 50, 100, and 500 bar and planetary masses of 0.107 (M_Mars), 0.3, 0.5, 0.7, 0.9, and 1 M_Earth, with mass indicated by color. The black dashed line is the result from Genda & Abe (2003) for an Earth-mass planet and the solid black line is a fit to our simulation results (see Section C). There is very little variation in the efficiency of loss as a function of ground velocity with atmospheric pressure and planetary mass, when ground velocity is normalized to the escape velocity. This confirms what has been assumed in other studies (Genda & Abe 2003, 2005; Schlichting et al. 2015), that the effect of planetary mass in the absence of an ocean is almost entirely accounted for by normalization to the escape velocity.

The efficiency of atmospheric loss for a given ground velocity can be significantly enhanced if the colliding bodies have part or all of their surfaces covered by water (Genda & Abe 2005). In the following sections, we compare our results to those of Genda & Abe (2005), and explore how the the efficiency of loss is dependent on the initial surface conditions (e.g., atmospheric pressure, ocean depth, etc.) and the mass of the planet.

4.2.1. Comparison to Previous Results

Figure 6(A) shows the efficiency of loss as a function of ground velocity for H₂ atmospheres of 300, 30, and 1 bar (dotted lines) above 3 km deep oceans on Earth-mass planets with surface temperatures of 300 K. The pressures at the base of ocean were approximately 600, 330, and 300 bar, respectively. For reference, loss in the no-ocean case is shown in black. The examples shown in Figure 6(A) were also explored by Genda & Abe (2005) and their results are shown as dashed lines. In their study, Genda & Abe (2005) considered H₂ atmospheres with six different atmospheric pressures, but we have chosen to only show three here to allow clear comparison with our results. At ground velocities less than ∼6 km s⁻¹ (0.54 v_esc) we find relatively good agreement between our results and those of Genda & Abe (2005), but at higher velocities our results diverge, particularly for higher-pressure atmospheres. This difference is due to the EOS for water used at high ground velocities in each study. Genda & Abe (2005) ran simulations using both the IAPWS EOS (Wagner 2002) and the Tilloston EOS (Tillotson 1962). At lower ground velocities the two EOSs gave relatively similar results, but the maximum pressure limit of the IAPWS EOS precluded its use for ground velocities above 6 km s⁻¹, where the pressure in the shocked stated exceeded the range of the EOS (marked by crosses on each loss curve in Figure 6(A)). It is therefore in the regime in which Genda & Abe (2005) only performed calculations using the Tillotson EOS in which our results significantly differ from theirs.

Zoom In Zoom Out Reset image size

To confirm that the only difference between our results and those of Genda & Abe (2005) is the water EOS used, we ran additional simulations using the Tillotson EOS for the ocean and found good agreement between our results when using the same EOS. Figure 6 shows an example set of simulations for an Earth-mass planet with a 300 bar, H₂ atmosphere over a 3 km ocean, with our simulations using the Senft & Stewart (2008) EOS shown by the dotted line, our results using the Tillotson EOS as a dashed–dashed–dotted line, and the results of Genda & Abe (2005) as a dashed line. The treatment of expanded states in the Tillotson EOS often leads to unphysical solutions at low densities, causing simulations to fail. As a result, very few of our calculations using the Tillotson EOS reached the prescribed runtime of 5000 s, which likely accounts for the slightly lower loss we calculate in some cases.

The Tillotson EOS is designed to model the behavior of material in the shocked state but does not provide a good description of material properties in lower-density, expanded states. This is particularly an issue in the multiphase liquid and vapor region where a minimum density cutoff is imposed which is typically a sizeable fraction of the reference density. In contrast, the water EOS from Senft & Stewart (2008) is designed for use in planetary collisions and includes a liquid-vapor phase region to more accurately describe expanded states. The efficiency of loss in the ocean case is dictated by the complex combination of multiple waves (see discussion below) and so the use of a high-quality EOS for water is critical to accurately determine loss. At high ground velocities, more of the water reaches expanded states and hits the minimum density cutoff, likely accounting for the lower loss at high ground velocities when using the Tillotson EOS. Given that the Senft & Stewart (2008) EOS provides a more accurate description of expanded states, our results are likely more realistic than those of Genda & Abe (2005) using the Tillotson EOS. For a discussion of the comparison between the Tillotson and more advanced EOSs, see Stewart et al. (2020).

4.2.2. Dependence on Ocean Depth and Atmospheric Pressure

To explore the effect of surface pressure and ocean depth, as well as planetary mass (Section 4.2.3), on the efficiency of loss as a function of ground velocity, we conducted simulations with every combination of six different planetary masses (0.107, 0.3, 0.5, 0.7, 0.9, and 1 M_Earth), seven surface pressures (p₀ = 1, 5, 10, 50, 100, 300, and 500 bar), and nine ocean depths (0.1, 0.5, 1, 2, 3, 5, 10, 20, and 30 km). We also conducted additional simulations for each planetary mass with 900 bar atmospheres and oceans of 0.1 km depth. Atmospheres were CO₂ (m_a = 44 g mol⁻¹, γ = 1.29) with surface temperatures of 300 K. Simulations were performed for ground velocities at 0.05 v_esc intervals between 0.05 and 0.95 v_esc.

Figure 7 shows the efficiency of loss from an Earth-mass planet for six example surface velocities for different atmospheric pressures (symbols) and ocean depths (x-axis). The amount of atmosphere and ocean lost is presented as the sum of the atmospheric and ocean loss (which we refer to as combined loss) with one being the total loss of atmosphere and two being the total loss of both ocean and atmosphere. For reference, the lower open gray triangle to the left of each panel shows the loss expected in the absence of an ocean. Note that the x-axis is the inverse of ocean height (1/H_oc).

Zoom In Zoom Out Reset image size

Variation in the atmospheric pressure and ocean depth can make the difference between almost zero and total loss of an atmosphere, and between zero and almost total loss of the ocean. Loss is more efficient from planets that initially have deeper oceans and/or lower-pressure atmospheres. For shallow oceans and high-pressure atmospheres, the efficiency of loss tends toward a low-loss limit where the efficiency of loss is the same as that in the no-ocean case (open gray triangle to the left of each panel in Figure 7). Increasing the ocean depth while keeping the atmospheric pressure constant leads to an increase in the efficiency of loss, until the loss plateaus at a high-loss limit for very deep oceans. Over the range of conditions we have considered only the lowest pressure atmospheres plateau in the high-loss limit. At lower ground velocities, where only the atmosphere is being lost, the value of the plateau is dependent on the initial atmospheric pressure. As we will describe below, this phenomenon is due to the fact that the lower the initial atmospheric pressure, the higher the velocity of the ocean surface upon release to the atmosphere and so the stronger the driver for atmospheric loss. When the atmosphere is totally lost, ocean loss in the high-loss limit is almost invariant of initial atmospheric pressure.

It is evident from Figure 7 that the physics of atmospheric loss in the presence of an ocean is more complicated than a simple impedance-match calculation (Section 2). Filled symbols to the left of each panel in Figure 7 show the loss determined by convolving the velocity of the ocean surface determined from an impedance-match solution (see Figure 12) with the parameterization for atmospheric loss in the case of no ocean (see Section C), i.e., the loss that would be expected if loss of the atmosphere was only controlled by the ocean surface driving a shock at the impedance-match velocity. The dark blue line on the left of each panel in Figure 7 shows a similar calculation except using the velocity of the ocean surface expected upon release of the ocean to very low pressure (1 Pa). The plateau in loss seen in our simulations is typically higher than that calculated assuming the impedance-match velocity as the driving velocity for loss, and additional processes must be at play.

To explain the dependence of loss on atmospheric pressure and the depth of the ocean, it is necessary to understand the dynamics of the system beyond the initial breakout of the shock (Section 2). Figure 8 shows examples of the early evolution of the ocean and atmosphere after breakout of the impact shock from the ground. Each column shows the evolution for planets with the same depth of ocean (3 km) and the same ground velocity (4 km s⁻¹), but with increasing initial atmospheric pressures going from left to right (1, 50, and 500 bar). Upon breakout from the planet, the system evolves as dictated by the impedance-match between the different layers, as described in Section 2. The shock propagates through the ocean, compressing the water and accelerating it to the ground velocity. When the shock front reaches the surface of the ocean the water releases to the impedance-match pressure between the water and the atmosphere (gray dashed lines in Figure 8), driving a shock wave into the atmosphere. The shock accelerates as it travels up the atmosphere and the initial evolution is similar to that seen in Figure 4 for the no-ocean case, but with the ocean surface taking the role of the ground.

Download figure:

Video Standard image High-resolution image

Lamentably, further evolution of the system after the initial transmission of the shock into the atmosphere complicates the simple impedance-match picture. As discussed in Section 2, the release wave propagating downwards into the ocean causes more of the ocean to accelerate. Furthermore, as the pressure at the base of the atmosphere falls (Figure 8), the ocean continues to decompress and accelerate and can reach velocities significantly above the impedance-match velocity at later times. The ocean is stretched between the rapidly moving ocean surface and the ground, and the pressure in the ocean decreases rapidly. In cases with low-pressure atmospheres and deep oceans (e.g., Figures 8(A) and (B)), the pressure in the ocean remains above that of the atmosphere for most, if not all, of the subsequent evolution and the atmospheric shock is well supported. The ocean surface continues to decompress, accelerating to velocities approaching those expected upon release of the ocean to a vacuum. This effect is responsible for the plateauing of loss in the high-loss limit.

However, in most cases the ballistic expansion of the ocean causes the ocean pressure to drop below that of the lower atmosphere after a time dependent on the initial surface conditions and strength of the shock (Figures 8(C)–(F)). The difference in pressure exerts a force to slow the ocean, and the ocean velocity is decreased by a series of pressure waves traversing the ocean. The shock is no longer fully supported in the atmosphere and a release wave from the bottom of the atmosphere can retard the acceleration of the upper fraction of the atmosphere, leading to decreased loss. In cases with either very thin oceans or very-high-pressure atmospheres (e.g., Figures 8(E) and (F)), pressure waves can equalize the pressure throughout the ocean before the shock wave has propagated far into the atmosphere. The velocity of the ocean surface slows to that of the ground, and the evolution of the atmosphere is very similar to that in the no-ocean case. This explains why the atmospheric loss tends to that in the no-ocean case in the low-loss limit (Figure 9).

Zoom In Zoom Out Reset image size

In all cases, for sufficiently high ground velocities continued expansion of the ocean, contributed to by a release wave from the top of the atmosphere, leads to a slow acceleration of the ocean surface and loss of the top of the ocean (e.g., Figure 8(A), top of the ocean on the right side of the yellow solid line). When the whole atmosphere is lost the ocean effectively releases to zero pressure and the initial pressure of the atmosphere becomes almost irrelevant. The high-loss limit for ocean loss is therefore relatively insensitive to atmospheric pressure.

The time at which the pressure in the ocean becomes less than that at the base of the atmosphere is critical in governing the efficiency of loss. The earlier in time this transition occurs, the earlier the ocean surface and atmosphere are slowed, and the lower the degree of atmospheric/ocean loss. The timing of this transition depends on two factors: the depth of the ocean, and the initial atmospheric pressure. For deeper oceans, the increase in the depth of the ocean layer as the ocean surface expands is a smaller fraction of the total depth of the ocean. The density, and hence pressure, of the ocean therefore falls more slowly. The dependence of loss on atmospheric pressure is more complicated as there are two competing effects. The atmospheric pressure dictates the impedance-match pressure and velocity of the ocean surface and base of the atmosphere. The lower the initial atmospheric pressure, the higher the impedance-match velocity of the ocean surface and the greater the driver for atmospheric loss. In addition, the lower the initial atmospheric pressure, the lower the pressure in the atmospheric shock. All else being the same, as the ocean expands it takes longer for the pressure in the ocean to fall below the pressure in the shocked lower atmosphere, leading to slowing of the ocean surface later in time. However, working against this effect is that the higher the surface velocity of the ocean, the more rapidly the ocean expands. The pressure in the ocean decreases more rapidly and so falls below the pressure in the lower atmosphere earlier in time, leading to an earlier slowing of the ocean surface by pressure waves. Determining the balance between these two effects of atmospheric pressure is nontrivial.

The efficiency of atmospheric loss is therefore controlled by the initial depth of the ocean and atmospheric pressure in two principle ways: the atmospheric pressure controls the impedance-match velocity between the ocean and atmosphere that drives enhanced loss; and a combination of the ocean depth and initial atmospheric pressure determine when the ocean surface begins to slow. Genda & Abe (2005) previously suggested that loss was dependent on the ratio of the mass of the atmosphere to the mass of the ocean (${ \mathcal R }={M}_{\mathrm{atm}}/{M}_{\mathrm{oc}}$). For a planet of a given mass, this mass ratio is proportional to the ratio of initial atmospheric pressure to ocean depth (p₀/H_oc):

$$\begin{eqnarray}&&{ \mathcal R }\sim \displaystyle \frac{{p}_{0}}{{H}_{\mathrm{oc}}}\displaystyle \frac{1}{g},\end{eqnarray} \tag{ 4 }$$

where g is the gravitational acceleration at the base of the atmosphere. We find that loss correlates well with both ${ \mathcal R }$ (Figure 9) and p₀/H_oc over a wide range of atmospheric pressures and ocean depths, including across the transition between the low- and high-loss regimes. This demonstrates the key role that atmospheric pressure and ocean depth play in loss.

The correlation of loss with ${ \mathcal R }$ provides an alternative way of understanding the limits on loss. In the high-${ \mathcal R }$ limit, the ocean is much less massive than the atmosphere and so the release of the ocean cannot provide sufficient momentum to drive any enhancement in atmospheric loss. The efficiency of loss is then the same as if it was just driven by the ground motion alone, and the whole atmosphere, or any amount of ocean, is not lost until the ground velocity is very close to the escape velocity. In the low-${ \mathcal R }$ limit, the ocean is much more massive than the atmosphere and so the atmosphere offers little impediment to the release of the ocean to very low pressures. The transition between the high- and low-loss regimes occurs when the mass of the atmosphere is comparable to that of the ocean (i.e., ${ \mathcal R }\sim 1$) and takes place over one to two orders of magnitude in ${ \mathcal R }$, dependent on the ground velocity. The transition occurs at higher ${ \mathcal R }$ as ground velocity increases as the ocean itself begins to be lost and the influence of the atmosphere diminishes.

We find that considering loss as a function of ${ \mathcal R }$ is more useful than p₀/H_oc when considering planets with different masses (Section 4.2.3). We therefore use ${ \mathcal R }$ as the principal control on loss from here on out. However, it is important to bear in mind the close relationship between the two ratios.

Scaling with ${ \mathcal R }$ does not capture all the factors that influence the relationship between ground velocity and loss. In Figures 9(A) and (B), the dependence of the upper limit of loss on atmospheric pressure is noticeable, but the scaling with M_atm ∼ p₀ of the x-axis means that the offset in loss at a given ${ \mathcal R }$ is smaller than at the same H_oc. In addition, loss in transition between the high- and low-${ \mathcal R }$ regimes does not scale perfectly with ${ \mathcal R }$, particularly in the regime where ocean is being lost, with loss being higher in cases with initially higher-pressure atmospheres at the same ${ \mathcal R }$. In such cases, the scaling with M_atm ∼ p₀ is overcorrecting for the effect of atmospheric pressure as the ocean is able to decompress to lower pressures with limited restriction from the atmosphere. Over the range of parameters we have considered, variation in atmospheric pressure leads to differences in loss that are much smaller than the overall ${ \mathcal R }$ effect, with the exception of in the ocean loss transition region when the variation due to atmospheric pressure can be ∼a third of the total variation in loss. If we were to consider cases with high-pressure atmospheres but much larger ocean depths, the simple scaling with ${ \mathcal R }$ would not well describe the value to which loss plateaus in the high-loss regime. The lower impedance-match velocity of higher-pressure atmospheres would result in loss plateauing at lower values than for lower-pressure atmospheres (Figure 7). This effect would cause deviation on the order of the total variation in loss for cases with a similar ${ \mathcal R }$. Considering lower-pressure atmospheres than those we examine here could compound this effect as they could plateau at higher loss fractions than the 1 bar minimum pressure we simulated. A wider range of parameters will be considered in future work but, for now, we advise caution when using the results of this work beyond the parameter regime simulated.

4.2.3. Dependence on Planetary Mass

The scaling of loss with ${ \mathcal R }$ holds well, with some caveats, when considering planets of different masses. Figure 10 shows the combined loss as a function of ${ \mathcal R }$ for six example velocities (normalized to v_esc) and six different mass planets (colors) between the mass of Mars (M_Mars) and the mass of Earth (M_Earth). The symbols are the same as in Figures 7 and 9, and colored lines show the results of our parameterization for atmospheric loss for each mass planet (Section C). The dependence on loss is generally well captured by scaling with ${ \mathcal R }$, with the exception that in the low-${ \mathcal R }$ regime when there is only partial atmospheric loss the plateau in loss is lower for lower-mass planets. For a ground velocity which is a given fraction of v_esc, the absolute ground velocity, and hence the strength of the shock, for a lower-mass planet is lower. The resulting impedance-match velocity for the ocean–atmosphere interface is a lower fraction of the escape velocity of the smaller planet, leading to less efficient loss (filled symbols to the left of each panel in Figure 10). The influence on loss is compounded by the fact that the atmospheric loss function is highly nonlinear at low velocities (Figure 5) while the impedance-match velocity is relatively linear with respect to absolute velocity. In the regime in which the entire atmosphere is lost, the loss of ocean is controlled by expansion of the ocean to low pressure and the maximum loss is relatively insensitive to planetary mass.

Zoom In Zoom Out Reset image size

There is also substantial deviation from a simple ${ \mathcal R }$ scaling in the transition region between the low- and high-loss regimes, with loss from more massive planets being more efficient. This variation is likely largely due to the sensitivity of the dynamics of loss to the shock and release path of water, which itself is a function of the absolute strength of the shock. The variation due to planetary mass is compounded by the variation due to initial atmospheric pressure/ocean depth (Section 4.2.2), and the difference in combined loss can be as much as 50% at the same ${ \mathcal R }$ in regions where the loss fraction is varying rapidly with ${ \mathcal R }$. These variations are the main cause of error in our parameterization of loss (Section C).

4.2.4. Effect of Atmospheric Composition

In the presence of an ocean, the effect of atmospheric composition on the relationship between ground velocity and loss is complex as the composition of the atmosphere can affect the efficiency of loss in a number of different ways. First, the composition of the atmosphere controls the compressibility of the shocked gas, and the impedance-match velocity and pressure of the ocean surface. For an ideal gas in the strong-shock limit

$$\begin{eqnarray}&&\displaystyle \frac{{\rho }_{{\rm{s}}}}{{\rho }_{0}}\approx \displaystyle \frac{\gamma +1}{\gamma -1},\end{eqnarray} \tag{ 5 }$$

where ρ_s and ρ₀ are the density of the shocked and unshocked material, respectively, and γ is the ratio of specific heat capacities of the gas. The compression of material due to the shock is entirely dependent on γ which varies from 1.25 to 1.4 for the gases we consider, resulting in a variation of density in the shocked gas of 6–9. The particle velocity, u_p, at a given shock pressure, p_s , in the strong-shock limit is

$$\begin{eqnarray}&&{u}_{{p}}\approx \sqrt{\displaystyle \frac{{{p}}_{{\rm{s}}}}{{\rho }_{0}}\left(\displaystyle \frac{2}{\gamma +1}\right)}.\end{eqnarray} \tag{ 6 }$$

The lower the initial density of the gas, the higher the particle velocity required to reach a given shock pressure. As a result, the impedance-match velocity at the ocean surface is higher for lighter gases while the impedance-match pressure is lower. In particular, H₂ is much less dense than gases such as N₂ and CO₂, and the initial velocity of the ocean surface can be tens of percent larger for H₂ than for the heavier gases. The lower pressure of the impedance-match solution in such cases means that the pressure at the base of the ocean can remain above that at the base of the atmosphere for longer, helping sustain the velocity of the ocean surface. The height of the atmosphere, and hence the time it takes for the shock to reach the top of the atmosphere, can vary significantly depending on its composition. For example, it can take the shock more than 10 times longer to reach the top of a H₂ atmosphere than a CO₂ atmosphere. As a result, the surface of the ocean and the ground slow earlier in the evolution relative to the progress of the shock through the atmosphere. So, although the velocity of the ocean surface is sustained for longer for H₂ atmospheres in absolute time, the ocean surface is typically slowed earlier relative to the evolution of the atmosphere, acting to reduce the efficiency of loss. Finally, the release wave from the top of the atmosphere reaches the ocean surface earlier during the loss of CO₂ atmospheres compared to H₂ atmospheres, and so the pressure in the ocean drops more rapidly and the ocean reaches higher velocities earlier in the evolution.

The net effect of these competing factors depends on the initial conditions and the regime of loss, and can be difficult to isolate. For example, Figure 11(A) shows the same cases as in Figure 6(A) for both H₂ and CO₂ atmospheres, which show some of the common tradeoffs. The loss for different heavy atmospheres (e.g., N₂ and CO₂) is very similar, but the loss of H₂ atmospheres can be significantly different. In the high-loss limit (e.g., yellow curves in Figure 11(A)), the shock in the atmosphere is supported late into the evolution in both the CO₂ and H₂ cases and the efficiency of loss is very similar. In the regime where only part of the atmosphere is lost, loss in the H₂ case is slightly greater due to the larger impedance-match velocity. However, this effect is muted by the continued decompression of the ocean surface, and the ocean surface in the CO₂ case reaches greater velocities than in the H₂ cases within a few seconds. In the high-loss regime when the whole atmosphere is lost, the loss of ocean in the CO₂ and H₂ cases is nearly identical as the low-mass atmosphere provides little impediment to the loss of the ocean, regardless of atmospheric composition.

Zoom In Zoom Out Reset image size

In the transition between the high- and low-loss limits (e.g., teal and blue lines in Figure 11(A)) the differences between CO₂ and H₂ atmospheres are more complex. At low ground velocities, loss of CO₂ atmospheres is more efficient largely because the ocean surface in the H₂ cases slows significantly before the shock reaches the top of the atmosphere. There is then an intermediate ground velocity regime where the situation is similar to that in the high-loss limit and loss of H₂ atmospheres is more efficient. Close to near-total loss of the atmosphere and, continuing into the ocean-loss regime, loss of CO₂ atmospheres and their oceans once again becomes more efficient. This effect is likely a result of acceleration caused by the more rapid drop in the pressure of the ocean in the CO₂ cases allowing the ocean to reach higher velocities. In cases where the ocean is almost entirely lost, the situation can reverse (e.g., teal line in Figure 11(A)). This may indicate that the higher impedance-match velocity in the H₂ case controls the loss of the deep ocean.

The difference in loss between CO₂ and H₂ atmospheres can be tens of percent, particularly in the ocean-loss regime. However, for the range of parameters we have explored, the loss of H₂ atmospheres is within, if sometimes at the extreme, of the variation we see between different H_oc, p₀, and M_p combinations for CO₂ atmospheres with similar ${ \mathcal R }$ (Figure 10). Therefore, although the effect of atmospheric composition can be significant, it is a second-order effect compared to the dominant ${ \mathcal R }$ scaling.

4.2.5. Effect of Surface Temperature

As discussed in Section 2, in the no-ocean case the ground velocity is set by an impedance match between the rock surface and the atmosphere. Figures 12(A)–(C) shows the ground velocity as a function of the particle velocity in the shock in the planet before breakout for atmospheres with different pressures, compositions, and temperatures (colored lines). For reference, the particle velocity of the ground on release to very low pressures (1 Pa, black dashed line), and the impedance-match velocity for release of the ground into an ocean (deeper blue line) are also shown. For high-temperature, low-pressure atmospheres the ground velocity is close to the low-pressure release velocity, as the impedance-match pressure is low enough that the forsterite release curve is very steep in u_p–p space by the time it intersects the shock Hugoniot of the gas. This can be seen in the example impedance-match calculations shown in u_p–p space in Figure 13(A) (similar to the schematics in Figures 2 and 3), where the black dashed lines show rock release curves and black symbols show the impedance-match solutions. For the same strength of shock in the planet, the impedance-match velocity decreases with increasing atmospheric pressure and decreasing temperature, and is also lower for heavier atmospheres (Figure 12(B)). Therefore, even though the dependence of loss on the ground velocity is largely insensitive to the atmospheric properties (Section 4.1 and Figure 5), it is easier to lose hotter, lighter, and lower-pressure atmospheres from terrestrial planets due to the higher ground velocity resulting from any given impact.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The effect of the varying impedance-match velocities on the magnitude of loss is shown in Figure 12. Panels (E)–(G) show the loss calculated by inputting the impedance-match velocities shown in Figures 12(A)–(C) into our parameterization of the relationship between ground velocity and atmospheric loss (Section C). Panel (H) shows similar results but for the four atmospheres considered in the 3D loss simulations of Kegerreis et al. (2019) utilizing the EOS of Hubbard & MacFarlane (1980). Panels (I)–(L) show the difference between the calculated loss from each of the atmospheres and the loss that would be driven by release of the surface to very low pressures (1 Pa, black dashed line in panels (A)–(H)). Figure 12 also shows an ocean line, which we will discuss in Section 5.2. The efficiency of loss for a given strength of impact shock can vary by 10 s % between different atmospheres (Figures 14(I)–(L)), purely as a function of the difference in the impedance-match velocity between the ground and the atmosphere. The ability of giant impacts to remove volatiles from planets is therefore dictated, in part, by the surface conditions on the colliding bodies before the impact.

Zoom In Zoom Out Reset image size

As a result of the dependence of ground velocity on atmospheric properties, care needs to be taken when combining u_G–loss scalings (Section C) with the results of 3D impact simulations and when applying the results of 3D impact loss simulations to accretion models. The peak ground motion in a 3D impact simulation will be dictated by the pressure of the atmosphere, or the lack of atmosphere, used in that simulation. When combining 1D simulations with the ground velocity from 3D impact simulations (Kegerreis et al. 2019) to calculate the loss expected from a planet with a different atmosphere than that used in the 3D simulation, the peak ground motion must be corrected to that dictated by the impedance match with the alternative atmosphere. Furthermore, results of directly calculated 3D loss simulations from a planet with an atmosphere of one composition, pressure, and temperature cannot necessarily be applied to loss from a planet with a different atmosphere without significant error. For example, the difference in the loss efficiency curves shown in Figures 12 (H) and (L) may be at least partly responsible for the difference in atmospheric loss between different mass atmospheres observed by Kegerreis et al. (2019, e.g., their Figure 7). The dependence of ground velocity on atmospheric properties complicates efforts to develop universal scaling laws for atmospheric loss.

Figure 14(A) shows the impedance-match velocity between the ocean and the atmosphere for a given strength of shock in the planet. The impedance-match velocities are much larger than that between the ground and atmosphere in the no-ocean case (Figures 12(A)–(D)), explaining the significant increase of loss achieved in the high-loss regime (Figure 10). However, as we discussed in Section 4.2, the loss efficiency can be significantly perturbed from that expected from an impedance-match calculation.

When there is an ocean on the preimpact body, the ground velocity is set, not by release of the ground to the atmosphere, but by an impedance match between the rock surface and the bottom of the ocean. The darker blue line in Figures 12(A)–(D) shows the ground velocity upon release to the ocean as a function of the shock particle velocity in the planet before breakout. Note that, although the Hugoniot depends on the initial pressure/density of the material, over the range of ocean depths and atmospheric pressures considered in this paper the forsterite and water Hugoniots are similar and there is little variation in the ground velocity upon release. The ground velocity in the ocean case can be tens of percent lower than the ground velocity in the no-ocean case for the same impact (Figures 12(I)–(L)).

The reduced ground velocity in the ocean case can lead to a previously unrecognized phenomenon: The presence of an ocean can actually reduce the efficiency of loss (see Genda & Abe 2005). In the large atmosphere/ocean mass ratio, low-loss regime, the loss efficiency is the same as the loss in the no-ocean case for the same ground velocity. However, the ground velocity in the ocean case, for a given impact, is lower than in the no-ocean case. Therefore, the amount of atmospheric loss in a given impact will be lower if an ocean is present than if it were absent. Figure 14(B) shows the fraction of atmosphere lost due to a given particle velocity depending on the ocean-to-atmosphere mass ratio. The black dashed line shows the equivalent relationship for low-pressure atmospheres in the absence of an ocean. The atmosphere-to-ocean mass ratio needs to be sufficiently small in order for the loss efficiency in the ocean case to exceed that in the no-ocean case. Figure 14(C) shows the critical mass ratio required for loss in the ocean case to equal that in the no-ocean case, assuming release of the ground to very low pressure (∼1 Pa). The mass of the ocean must be at least comparable to the mass of the atmosphere in order for the presence of an ocean to enhance loss, over that in the no-ocean case. The critical mass ratio is lower for higher-pressure atmospheres as the release velocity is comparatively lower in the no-ocean case.

We have shown that the surface conditions on colliding bodies before giant impacts can have a strong influence on the efficiency of atmospheric and ocean loss. In the no-ocean case, hotter, lower-pressure, and lower-molecular-weight atmospheres are more easily lost due to the effect of atmospheric properties on the impedance-match velocity between the ground and atmosphere. If terrestrial planets grow to a large enough size while the nebula is still present, they can accrete a primary H₂–He atmosphere (e.g., Ikoma & Genda 2006), as inferred in the case of many exoplanets (e.g., Jontof-Hutter 2019). After the nebula dissipates, such a light atmosphere would be comparatively susceptible to loss by giant impacts, in addition to more generally considered thermal or stellar-radiation-driven loss (e.g., Lammer et al. 2014). Over time, the atmospheres of planets can become dominated by heavier-molecular-weight species produced by degassing of accreted solids and from the planet's interior. These heavier-molecular-weight atmospheres, often called secondary atmospheres, would be less susceptible to loss by giant impacts compared to H₂–He atmospheres. For atmospheres of all compositions, depending on the efficiency of volatile accretion, repeat impact events could reduce the atmospheric pressure, providing a positive feedback for increasingly efficient loss through accretion. The ability of giant impacts to remove atmospheric volatiles thus can evolve substantially during planet formation as the atmospheres of planets change composition, mass, and oxidation state.

Toward the end of planet formation, the cadence of giant impacts and the flux of planetesimal accretion in our solar system are both low enough to allow significant cooling between giant impacts. Thus, water may condense on protoplanets, given a suitable oxidation state, in the time between giant impacts (Abe & Matsui 1988). This may also be the case in many exosystems, but will depend strongly on the timescale of accretion. Therefore, giant impacts that occur later in accretion are likely to be between planets with preimpact oceans. We have shown that, over the range of atmospheres we considered, loss from planets with oceans is only weakly dependent on the atmospheric composition, temperature, and pressure, and is largely dominated by the atmosphere-to-ocean mass ratio (as proposed by Genda & Abe 2005). We find a relatively rapid transition between a high- and a low-loss regime, over one to two orders of magnitude in the atmosphere-to-ocean mass ratio. Contrary to previous thinking, the presence of an ocean can reduce the efficiency of loss in the low-loss regime. Condensation of an ocean could then actually protect the surface volatile inventory from loss by giant impacts. However, if the ocean is sufficiently massive, with the critical mass being typically within a factor of a few of the atmospheric mass, the presence of an ocean can significantly increase the efficiency of atmospheric loss, in agreement with Genda & Abe (2005). Therefore, a critical factor in understanding the evolution of planetary atmospheres during accretion is determining whether they are in the high- or low-loss regimes.

In our solar system, the compositions of chondritic meteorites are often used as proxies for the compositions of the building blocks of Earth. For most chondrites, if their full budget of H, C, and N were converted into a CO₂ and N₂ atmosphere over a pure H₂O ocean, the resulting atmosphere-to-ocean mass ratio would be of order 10⁰–10¹ (see Table 1). For all but the highest ground velocities, such an atmosphere and ocean would be at the edge of the low-loss efficiency regime and only high-energy impacts would drive substantial atmospheric loss (Kegerreis et al. 2020a, 2020b). However, since loss of ocean is less efficient than loss of atmosphere (Section 7.2), there is a positive feedback where loss events push the atmosphere-to-ocean mass ratio lower, increasingly the susceptibility of the remaining atmosphere to loss in subsequent giant impacts. The rapid transition in the efficiency of loss with decreasing atmosphere-to-ocean mass ratio could result in planets reaching a tipping point in volatile evolution once their atmosphere-to-ocean mass ratio becomes sufficiently small to be in the high-loss regime. Such an effect relies on a planet receiving a sufficient number of giant impacts with a relative timing that allows an ocean to condense between impacts. The strong sensitivity to the order and timing of giant impacts could result in planets with similar accretionary histories obtaining substantially different final volatile budgets.

Table 1. The Likely Volatile Budgets of Accreting Protoplanets Span from the High- to the Low-loss Regimes

Reservoir	${ \mathcal R }$ (oxidized)	${ \mathcal R }$ (reduced)
Earth's surface today	3.6 × 10−3	⋯
Earth's exosphere	0.21	0.13
BSE (Marty 2012)	0.7 ± 0.5	0.4 ± 0.5
BSE (Halliday 2013)	0.3 ± 0.1	0.2 ± 0.1
CC (Marty 2012)	1.1 ± 0.3	0.7 ± 0.3
CI	8.85	5.71
H	12.8	8.24
L	17.2	11.0
LL	13.6	8.75
EH	15.9	10.4

Notes. Presented are the atmosphere-to-ocean mass ratios for the Earth today, and different reservoirs or meteorite groups assuming that all the carbon and nitrogen were partitioned into the atmosphere and the hydrogen was in the ocean. Estimates assuming an N₂ and CO₂ (oxidized) or N₂ and CO (reduced) atmospheres are given in different columns. Estimate of Earth's exosphere volatile budget from Marty (2012). Bulk silicate Earth (BSE) estimates from Marty (2012) and Halliday (2013). Average carbonaceous chondrite (CC) estimated by Marty (2012) with data from Kerridge (1985) and Robert (2003). CI values based on Orgueil Lodders (2003). Average composition of ordinary chondrites (H, L, LL) taken from compilation by Schaefer & Fegley (2007). Enstatite (EH) chondrite values based on Indarch with data from Wiik (1956), Moore & Lewis (1966), Moore & Gibson (1969), and Grady et al. (1986) as compiled by Schaefer & Fegley (2017).

Download table as:  ASCIITypeset image

Beyond just the total volatile budget, critical to the loss or retention of volatile elements is the partitioning of volatiles between reservoirs in forming planets. As an example, the present-day Earth has an atmosphere-to-ocean mass ratio of ∼3.6 × 10⁻³, firmly in the high-loss efficiency regime. However, if all the H, C, and N in the present-day atmosphere, ocean, and sedimentary rocks (exosphere) were present as a CO₂ and N₂ atmosphere and a H₂O ocean, the atmosphere-to-ocean mass ratio would be ∼0.21 (Table 1), which is in the transition between the high- and low-loss regimes. Similar calculations using two different estimates for the volatile content of the entire bulk silicate Earth (BSE) give atmosphere-to-ocean mass ratios of 0.7 ± 0.6 (Marty 2012) and 0.3 ± 0.15 (Halliday 2013). The oxidation state of a planet's surface also plays a role in dictating the efficiency of atmospheric loss. For example, if CO was the dominant carbon species in the atmospheres calculated above the atmospheric mass would be ∼a third lower (right column of Table 1).

How volatiles partition between reservoirs, and the oxidation state of the surface in the period between giant impacts depend on a number of factors, including the separation of volatiles from the silicate during condensation of the silicate vapor produced in the giant impact (Stewart et al. 2018; Lock & Stewart 2020; Caracas & Stewart 2023), degassing of the mantle during magma ocean solidification (see, e.g., Elkins-Tanton 2012; Bower et al. 2019), dissolution of carbon in the ocean and potentially precipitation and storage of abiogenic (or even potentially biogenic) carbonates, weathering of crust transferring volatiles into the crust or sediments, and the thermal state of the atmosphere and surface (e.g., Zahnle et al. 2010). Given the difficulty in understanding each of these processes and the interactions between them, it is likely not possible to precisely determine the mass and composition of ocean and atmosphere before each impact. We therefore advocate taking a statistical approach to exploring the effect of giant impacts on the volatile budgets of terrestrial planets, making use of the high- and low-loss efficiency regimes described here to bound the evolution of planetary volatile budgets.

One of the main outstanding questions in regard to Earth's chemistry is that the C/N and H/N ratios of the BSE are significantly fractionated from those of known chondrites (Halliday 2013). It has been proposed that preferential loss of atmosphere to ocean in impacts could fractionate N and C from H (Tucker & Mukhopadhyay 2014). If C was dissolved in the ocean or stored in the crust or sediments, more C could be retained, leading to N/C fractionation.

Our results confirm that giant impacts drive significantly more loss of atmospheric species than loss of water. In the high-loss efficiency regime, total atmospheric loss can be achieved from sections of the colliding bodies where ground motions reach only ∼0.35 v_esc for an Earth-mass planet, whereas substantial ocean loss requires ground velocities that are a large fraction of v_esc. The amount of ocean loss can also be reduced by slowing of the ground by release waves within the planet, a process which atmospheric loss is relatively insensitive to (Section 6). If the proto-Earth, or the planetary embryos that accreted to form it, underwent a giant impact when an ocean was present, there is significant potential for H/N and C/N fractionation. Such an effect would be particularly strong if the ratio of atmospheric to ocean mass was low enough to be in the high-loss efficiency regime. Smaller impacts can also lead to significant atmospheric loss (Schlichting et al. 2015), but the degree of ocean/atmosphere fractionation in such events has not been determined. If loss by giant impacts is required to explain the Earth's H/N and C/N ratios, this would imply that Earth experienced (potentially multiple) giant impacts late in formation when it had an ocean and potentially a low enough atmosphere-to-ocean mass ratio to make atmospheric loss efficient.

Rayleigh–Taylor instabilities develop when a fluid attempts to accelerate (i.e., push) another fluid of greater density than itself. The growth of the instability disrupts the material interface and significantly reduces the acceleration of the denser fluid due to the forcing of the lighter fluid. If such instabilities were to develop at the ground–ocean, ground–atmosphere, or ocean–atmosphere boundaries during the process of atmospheric loss the efficiency of loss could be significantly reduced. Our 1D simulations cannot capture instabilities and so it is necessary for us to compare the density of the different layers during our simulations to determine the likelihood for Rayleigh–Taylor instabilities.

The ground–ocean boundary is not at risk of experiencing Rayleigh–Taylor instabilities over the range of parameters considered in this work. As discussed in Section 6, we do not directly model the planet in our simulations, but we can calculate the expected density of the ground by using the pressure at the base of the ocean/atmosphere from our calculations combined with the calculated release curve of forsterite (Figure 14). We find that the density of the ground is never lower than the base of the ocean, even in cases where the ground begins to vaporize, and so we would expect the ground–ocean boundary to remain stable.

Determining the stability of the ground–atmosphere and ocean–atmosphere boundaries is more difficult due to the limitations of the available EOSs. High-quality EOSs for heavy gases such as N₂, CO₂, and air mixtures that cover the range of conditions considered here are not publicly available or easily attainable. In this work, we therefore decided to use an ideal gas EOS for atmospheres. For lower-pressure atmospheres (< ∼150 bar for CO₂, ∼4 kbar for H₂) the density of the ideal gas Hugoniot is lower than the density of forsterite and water on the release curve at the impedance-match point, even in cases where the ground/atmosphere begins to vaporize. However, for higher pressures the density of the ideal gas EOS exceeds the density of the water and even forsterite. This is likely due to the fact that the ideal gas EOS significantly overestimates the density of gases at high pressure as it neglects the volumes and interactions of particles. More sophisticated EOSs for gases (Span & Wagner 1996; Lemmon et al. 2000; Span et al. 2000) predict a much lower density for the Hugoniot of atmospheric gases, below that of the forsterite and water at the impedance-match pressure. However, the maximum pressure covered by these EOSs preclude using them to calculate shocks over much of the range considered here. It is therefore likely that the ocean–atmosphere and ground–atmosphere boundaries are indeed stable to Rayleigh–Taylor instabilities, but wider-ranging, high-quality gas EOSs are required to confirm this.

We have produced scaling laws for loss as a function of ground velocity (Section C) and scripts to calculate impedance-match velocities (available through GitHub and Zenodo), which allow for linking of our results with ground velocity distributions calculated from 3D impact simulations (e.g., Yalinewich & Schlichting 2018; Kegerreis et al. 2019). We hope that this will allow for investigations of the effect of surface conditions on global atmospheric loss which are currently not possible due to computational limitations or expense.

When combining our results with 3D simulations there are two particularly important factors to consider. First, it is important to correct the ground motion from the 3D simulation to the impedance-match velocity for the atmosphere/ocean under consideration from that used in the 3D simulations. This can be done with the scripts and functions made available with this manuscript.

Second, by considering a 1D system in our simulations, we have inherently assumed that the shock in the planet breaks out perpendicular to the ground. That is, that the ground velocity is always vertical relative to the local surface. In reality, the shock wave will be traveling at an angle relative to the surface. On breakout of the shock, the wave will be refracted depending on the impedance of the materials either side of the boundary (Henderson 1989). It is not clear what the best approach is to account for a nonperpendicular incidence angle, but we suggest that a correction to the local escape velocity due to the component of the shock particle velocity parallel to the surface may be adequate. However, different methods will need to be tested using full 3D simulations as ground truth before proceeding with hybrid 1D-3D models.

The core of the code is the solution of the fundamental fluid dynamics equations in 1D. Conservation of momentum can be written as

$$\begin{eqnarray}&&\rho a=-\displaystyle \frac{\partial p}{\partial x}-\displaystyle \frac{\partial q}{\partial x}+\rho g,\end{eqnarray} \tag{ A1 }$$

where ρ is density, a is acceleration, p is pressure, x is the Lagrangian spatial coordinate, g is gravitational acceleration, and q is the viscous stress. For our purposes, q is simply the artificial viscosity. The acceleration in a Lagrangian reference frame is given by

$$\begin{eqnarray}&&a=\displaystyle \frac{\partial u}{\partial t},\end{eqnarray} \tag{ A2 }$$

where u is the velocity:

$$\begin{eqnarray}&&u=\displaystyle \frac{\partial x}{\partial t}.\end{eqnarray} \tag{ A3 }$$

Conservation of energy is the balance between the rate of change of internal energy and the rate at which work is being done,

$$\begin{eqnarray}&&\rho \displaystyle \frac{\partial \epsilon }{\partial t}=(p+q)\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial t},\end{eqnarray} \tag{ A4 }$$

where is the specific internal energy.

In WONDY, the conservation equations are solved using a finite-difference method. The fluid is divided into Lagrangian zones (or cells) and the parameters of each cell are advanced in incremental time steps using the conservative equations. The bulk fluid properties (m, ρ, q, p, , c, T) are defined at the center of each zone and the kinematic variables (a, u, x) are defined at the boundaries of the zones. c is adiabatic sound speed, T is the temperature, m is the mass of a given Lagrangian zone, and the other variables are described above. For our purposes all kinematic variables are defined in reference to the center of the planet. Linear interpolation is used to find values between these points, e.g., to find the bulk fluid properties at the edge of a zone. Zones are labeled in order with an index, j, starting at 1, with the lower boundary referred to with the index 0. The notation we employ here, after Kipp & Lawrence (1982), for the fluid properties is ${\psi }_{j}^{n}$, which donates an arbitrary quantity, ψ, at the jth zone and the nth time step. Integer j denote values of the quantity at the upper boundary of a zone whereas half-integer values indicate quantities defined at the center of the zone. Similarly, integer n indicate values known at that time step and half-integer values indicate quantities defined halfway between time steps. Bulk fluid properties are defined at ${\psi }_{j-\tfrac{1}{2}}^{n}$ and the kinematic variables at ${\psi }_{j}^{n}$, with the exception of u, which is defined at ${u}_{j}^{n+\tfrac{1}{2}}$. These choices of when and where to define variables are for numerical convenience, which will become evident below. The second-order centered-difference method is used to calculate both time and spatial derivatives of all quantities.

Using this notation the fundamental equations can be defined in finite-difference form. The finite-difference momentum equation is

$$\begin{eqnarray}\begin{array}{rcl}{a}_{j}^{n} & = & 2\left[\displaystyle \frac{({p}_{j-\tfrac{1}{2}}^{n}+{q}_{j-\tfrac{1}{2}}^{n})-({p}_{j+\tfrac{1}{2}}^{n}+{q}_{j+\tfrac{1}{2}}^{n})}{{\rho }_{j+\tfrac{1}{2}}^{n}({x}_{j+1}^{n}-{x}_{j}^{n})+{\rho }_{j-\tfrac{1}{2}}^{n}({x}_{j}^{n}-{x}_{j-1}^{n})}\right]\\ & & -\,\displaystyle \frac{{{GM}}_{p}}{{\left({x}_{j}^{n}\right)}^{2}},\end{array}\end{eqnarray} \tag{ A5 }$$

where g_j ⁿ has been substituted for a radial gravity profile for a planet of mass M_p . G is the universal gravitational constant. Equation (A5) can be used to calculate the acceleration at the nth time step from quantities that are already known from that time step. The velocity at the next half time step can then be determined from the finite-difference form of Equation (A2):

$$\begin{eqnarray}&&{u}_{j}^{n+\tfrac{1}{2}}={u}_{j}^{n-\tfrac{1}{2}}+\displaystyle \frac{1}{2}\left({t}^{n+1}-{t}^{n-1}\right){a}_{j}^{n}.\end{eqnarray} \tag{ A6 }$$

Similarly, from Equation (A3) the position of zone j at the next full time step is

$$\begin{eqnarray}&&{x}_{j}^{n+1}={x}_{j}^{n}+\left({t}^{n+1}-{t}^{n}\right){u}_{j}^{n+\tfrac{1}{2}}.\end{eqnarray} \tag{ A7 }$$

Having used Equations (A5), (A6), and (A7) to determine the kinematic variables for the next time step, the code recalculates the thermodynamic properties for the zone. Conservation of mass in spherical coordinates provides the updated density:

$$\begin{eqnarray}&&{\rho }_{j-\tfrac{1}{2}}^{n+1}=\displaystyle \frac{3{m}_{j-\tfrac{1}{2}}}{4\pi \left[{\left({x}_{j}^{n+1}\right)}^{3}-{\left({x}_{j-1}^{n+1}\right)}^{3}\right]}.\end{eqnarray} \tag{ A8 }$$

To reduce round-off error, the code actually implements the following, equivalent, expression:

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{j-\tfrac{1}{2}}^{n+1} & = & \left\{\displaystyle \frac{1}{{\rho }_{j-\tfrac{1}{2}}^{n}}\right.\\ & & +\,{\left.\displaystyle \frac{4\pi {\rm{\Delta }}{t}^{n+\tfrac{1}{2}}}{3{m}_{j-\tfrac{1}{2}}}\left[{\xi }_{j}^{n+\tfrac{1}{2}}{u}_{j}^{n+\tfrac{1}{2}}-{\xi }_{j-1}^{n+\tfrac{1}{2}}{u}_{j-1}^{n+\tfrac{1}{2}}\right]\right\}}^{-1},\end{array}\end{eqnarray} \tag{ A9 }$$

where

$$\begin{eqnarray}&&{\xi }_{j}^{n+\tfrac{1}{2}}={\left({x}_{j}^{n+1}\right)}^{2}+{x}_{j}^{n+1}{x}_{j}^{n}+{\left({x}_{j}^{n}\right)}^{2},\end{eqnarray} \tag{ A10 }$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}={t}^{n+1}-{t}^{n}.\end{eqnarray} \tag{ A11 }$$

The internal energy and pressure for the next time step are found by a combination of the energy equation and the EOS of the material. The finite-difference form of the energy equation is

$$\begin{eqnarray}&&\begin{array}{l}{\epsilon }_{j-\tfrac{1}{2}}^{n+1}={\epsilon }_{j-\tfrac{1}{2}}^{n}\\ +\,2\left({p}_{j-\tfrac{1}{2}}^{n+1}+{p}_{j-\tfrac{1}{2}}^{n}+2{q}_{j-\tfrac{1}{2}}^{n+\tfrac{1}{2}}\right)\displaystyle \frac{\left({\rho }_{j-\tfrac{1}{2}}^{n+1}-{\rho }_{j-\tfrac{1}{2}}^{n}\right)}{{\left({\rho }_{j-\tfrac{1}{2}}^{n+1}-{\rho }_{j-\tfrac{1}{2}}^{n}\right)}^{2}}.\end{array}\end{eqnarray} \tag{ A12 }$$

We will examine the solution of this equation for each of the EOSs used in this study in Appendix A.2.

Resolving shocks in numerical codes is challenging due to the very rapid change in material properties and velocity across the shock front. The natural viscosity of the fluids we consider is low, leading to very thin shock fronts that cannot be feasibly resolved at planetary scales. Artificial viscosity is used in WONDY to smear the shock front over several zones and avoid discontinuities. WONDY includes both a quadratic viscosity,

$$\begin{eqnarray}&&{q}_{1}=\rho {b}_{1}^{2}{\left(\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial t}\right)}^{2},\end{eqnarray} \tag{ A13 }$$

and a linear viscosity,

$$\begin{eqnarray}&&{q}_{2}={b}_{2}c\left(\displaystyle \frac{\partial \rho }{\partial t}\right),\end{eqnarray} \tag{ A14 }$$

with the total viscosity given as a simple sum of q₁ and q₂ (i.e., q = q₁ + q₂). b₁ and b₂ are parameters that control the strength of the artificial viscosity. To make the artificial viscosity insensitive to the absolute spatial scales, b_i are scaled to the size of the zone,

$$\begin{eqnarray}&&{b}_{i}={B}_{i}{\rm{\Delta }}x,\end{eqnarray} \tag{ A15 }$$

so the true constants, B_i , are approximately equal to the number of zones that any disturbance will be smoothed over. The values for all constants used in the code are given in Table 4 and the sensitivity of our results to each parameter is examined in Appendix B.1.

q₁ is large only when the rate of change is large (i.e., in the shock front) and small elsewhere so it is used mainly to control gradients in the shock front and avoid discontinuities. q₂ is more effective at controlling spurious oscillations elsewhere where q₁ is negligible. Both forms of the viscosity are only used when the material is in compression, i.e., where

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}\lt 0.\end{eqnarray} \tag{ A16 }$$

The finite-difference form of the combined viscosity is

$$\begin{eqnarray}&&\begin{array}{l}{q}_{j-\tfrac{1}{2}}^{n+\tfrac{1}{2}}=\displaystyle \frac{{\rho }_{j-\tfrac{1}{2}}^{n+1}+{\rho }_{j-\tfrac{1}{2}}^{n}}{2}\\ \left[{B}_{2}\left(\displaystyle \frac{{x}_{j}^{n+1}-{x}_{j-1}^{n+1}+{x}_{j}^{n}-{x}_{j-1}^{n}}{2}\right){c}_{j-\tfrac{1}{2}}^{n}\left(\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial t}\right)\right.\\ \quad +{B}_{1}^{2}{\left(\displaystyle \frac{{x}_{j}^{n+1}-{x}_{j-1}^{n+1}+{x}_{j}^{n}-{x}_{j-1}^{n}}{2}\right)}^{2}\\ \qquad \times \,\left.\left(\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial t}\right)\left|\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial t}\right|\right],\end{array}\end{eqnarray} \tag{ A17 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial t}=\displaystyle \frac{2\left({\rho }_{j-\tfrac{1}{2}}^{n+1}-{\rho }_{j-\tfrac{1}{2}}^{n}\right)}{{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}\left({\rho }_{j-\tfrac{1}{2}}^{n+1}+{\rho }_{j-\tfrac{1}{2}}^{n}\right)},\end{eqnarray} \tag{ A18 }$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}={t}^{n+1}-{t}^{n}.\end{eqnarray} \tag{ A19 }$$

The time step is chosen to ensure numerical stability. The criterion for linear stability of the difference equations used in the WONDY code (Hicks 1978) is

$$\begin{eqnarray}&&\begin{array}{c}{\rm{\Delta }}t\lt {\rm{\Delta }}x\left[{B}_{2}c+2{B}_{1}^{2}\left|\frac{\dot{\rho }}{\rho }\right|{\rm{\Delta }}x\right.\\ \quad {\left.+\sqrt{{\left({B}_{2}c+2{B}_{1}^{2}\left|\frac{\dot{\rho }}{\rho }\right|{\rm{\Delta }}x\right)}^{2}+{c}^{2}}\,\right]}^{-1}.\end{array}\end{eqnarray} \tag{ A20 }$$

At the end of each cycle, the value of the right-hand side of Equation (A20) is evaluated for each cell to determine the maximum time step that can be used to ensure stability. Since the criteria in Equation (A20) are based on a linear stability analysis, and then assumed to apply generally, the critical value for each cell is scaled by a constant, ${K}_{{{\rm{t}}}_{1}}$, to account for any nonlinear effects and ensure stability, such that

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{t}_{j-\tfrac{1}{2}}^{n+\tfrac{3}{2}}={K}_{{{\rm{t}}}_{1}}{\rm{\Delta }}{x}_{j-\tfrac{1}{2}}^{n+1}\left[{B}_{2}{c}_{j-\tfrac{1}{2}}^{n+1}+2{B}_{1}^{2}\left|\displaystyle \frac{\dot{\rho }}{\rho }\right|{\rm{\Delta }}{x}_{j-\tfrac{1}{2}}^{n+1}\right.\\ {\left.+\,\sqrt{{\left({B}_{2}{c}_{j-\tfrac{1}{2}}^{n+1}+2{B}_{1}^{2}\left|\displaystyle \frac{\dot{\rho }}{\rho }\right|{\rm{\Delta }}{x}_{j-\tfrac{1}{2}}^{n+1}\right)}^{2}+{\left({c}_{j-\tfrac{1}{2}}^{n+1}\right)}^{2}}\,\right]}^{-1}.\end{array}\end{eqnarray} \tag{ A21 }$$

The next time step is then set as

$$\begin{eqnarray}&&{\rm{\Delta }}{t}^{n+\tfrac{3}{2}}=\min \left({\rm{\Delta }}{t}_{j-\tfrac{1}{2}}^{n+\tfrac{3}{2}},{K}_{{{\rm{t}}}_{2}}{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}\right),\end{eqnarray} \tag{ A22 }$$

to ensure stability and to limit the increase in the time step to a factor ${K}_{{{\rm{t}}}_{2}}$.

Although we will discuss the code in dimensional parameters here, these equations are implemented in the code in a dimensionless form to reduce numerical errors. All the parameters, with the exception of temperature, are nondimensionalised using combinations of the ground velocity, u_G, the radius of the planet, R_p, and p₀.

Note that the gravity field in the 1D calculation is assumed to be fixed and radial. Since the time for loss is small (∼100 s) compared to the timescale of the impact, any change to the field over this period is minimal, but the field could still be distorted by the deformation of the target and the presence of the impactor. As long as the gravity field does not vary significantly from a 1/r² dependence, a local escape velocity in the direction of the ground motion should be sufficient to scale our results to calculate the loss in 3D impact simulations (Kegerreis et al. 2020a, 2020b).

In this work, we have used a variety of EOSs, in particular for water. The EOS dictates the method used to solve the finite-difference form of the energy equation (Equation (A12)). We will now outline the EOSs used and their implementation.

A.2.1. Ideal Gas Equation of State

The ideal gas EOS is the simplest EOS for gases and its implementation is straightforward. The EOS can be written in the form

$$\begin{eqnarray}&&p=(\gamma -1)\rho \epsilon ,\end{eqnarray} \tag{ A23 }$$

where γ is the ratio c_p/c_V, and c_p and c_V are the specific heat capacities at constant pressure and volume, respectively. This linear relation between pressure and energy allows for an analytical solution to the energy equation,

$$\begin{eqnarray}&&\begin{array}{l}{\epsilon }_{j-\tfrac{1}{2}}^{n+1}=\left[{\epsilon }_{j-\tfrac{1}{2}}^{n}\right.\\ \ +\,\left.2\left({p}_{j-\tfrac{1}{2}}^{n}+2{q}_{j-\tfrac{1}{2}}^{n+\tfrac{1}{2}}\right)\displaystyle \frac{\left({\rho }_{j-\tfrac{1}{2}}^{n+1}-{\rho }_{j-\tfrac{1}{2}}^{n}\right)}{{\left({\rho }_{j-\tfrac{1}{2}}^{n+1}-{\rho }_{j-\tfrac{1}{2}}^{n}\right)}^{2}}\right]\\ \quad {\left[1-(\gamma -1){\rho }_{j-\tfrac{1}{2}}^{n+1}\right]}^{-1},\end{array}\end{eqnarray} \tag{ A24 }$$

which exclusively uses quantities that are known from the previous time step. The pressure at the next time step is given by

$$\begin{eqnarray}&&{p}_{j-\tfrac{1}{2}}^{n+1}=(\gamma -1){\rho }_{j-\tfrac{1}{2}}^{n+1}{\epsilon }_{j-\tfrac{1}{2}}^{n+1}.\end{eqnarray} \tag{ A25 }$$

The only parameter that needs to be specified for an ideal gas EOS is γ, although other parameters are needed to initialize the atmosphere (see Appendix A.4).

A.2.2. Tabulated Equations of State

We have also added to WONDY the ability to use a tabulated EOS in the form of a T – ρ table. In this work, we have used the water table form Senft & Stewart (2008), a high-quality water EOS based partially on the IAPWS water EOS (see Appendix A.2.3) but extended to the high pressures relevant for planetary impacts.

Use of a tabulated EOS requires a numerical solution to the energy equation. We solve for the temperature at each time step that solves the energy equation using a Newton–Raphson method. The solution is found iteratively using

$$\begin{eqnarray}&&{T}_{i+1}={T}_{i}-\displaystyle \frac{f({T}_{i})}{f^{\prime} ({T}_{i})},\end{eqnarray} \tag{ A26 }$$

where the subscripts mark the number of the iteration,

$$\begin{eqnarray}&&\begin{array}{l}f({T}_{i})=\epsilon ({\rho }^{n+1},{T}_{i})-{\epsilon }^{n}\\ \ -\,\displaystyle \frac{\left({\rho }^{n+1}-{\rho }^{n}\right)}{{\left({\rho }^{n+1}-{\rho }^{n}\right)}^{2}}\left(p({\rho }^{n+1},{T}_{i})-{p}^{n}-2{q}^{n+\tfrac{1}{2}}\right),\end{array}\end{eqnarray} \tag{ A27 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}f^{\prime} ({T}_{i})={c}_{V}({\rho }^{n+1},{T}_{i})\\ \ -\,\displaystyle \frac{\left({\rho }^{n+1}-{\rho }^{n}\right)}{{\left({\rho }^{n+1}-{\rho }^{n}\right)}^{2}}\left[{\left.\displaystyle \frac{\partial \,p}{\partial T}\right|}_{\rho }({\rho }^{n+1},{T}_{i})\right].\end{array}\end{eqnarray} \tag{ A28 }$$

Note we have dropped the spatial subscripts for clarity as all values required are for the zone under consideration. The values of , p, c_V , and ∂p/∂T∣_ρ are found by linear interpolation of the tabulated values. The derivatives are precalculated at each point in the table using a second-order modified central difference method for unevenly spaced points. The initial value for the iteration is determined by assuming an isentropic volume change from the previous time step. The change in temperature, ΔT, for a given density change, Δρ, is given by

$$\begin{eqnarray}&&{\rm{\Delta }}T={\left.\displaystyle \frac{\partial \,T}{\partial \rho }\right|}_{S}{\rm{\Delta }}\rho ,\end{eqnarray} \tag{ A29 }$$

which can be written

$$\begin{eqnarray}&&{\rm{\Delta }}T={\left.\displaystyle \frac{T}{{c}_{V}{\rho }^{2}}\displaystyle \frac{\partial \,p}{\partial T}\right|}_{\rho }{\rm{\Delta }}\rho ,\end{eqnarray} \tag{ A30 }$$

which is straightforward to calculate from the tabulated values. The solution is considered converged when

$$\begin{eqnarray}&&\left|\displaystyle \frac{{T}_{i+1}-{T}_{i}}{{T}_{i+1}}\right|\lt {10}^{-10}.\end{eqnarray} \tag{ A31 }$$

A maximum of 1000 iterations are allowed but the routine typically converges within two to three iterations.

A.2.3. IAPWS Water Equation of State

A.2.4. Tillotson Equation of State

The Tillotson EOS (Tillotson 1962) has been widely used in shock physics and the study of planetary collisions. The EOS is based on empirical fits along the Hugoniot forced to converge to the Thomas–Fermi limit at high pressures. Although the EOS is designed for use primarily in compression, it includes parameters to approximate the behavior of material upon release. For this work, we use the parameters from O'Keefe & Ahrens (1982) for water.

Following Melosh (1989), the EOS implemented here has four different regimes: compressed states (ρ/ρ₀ ≥ 1, where ρ₀ is the reference density), cold expanded states (ρ/ρ₀ < 1 and < _iv , where _iv is the energy of incipient vaporization), hot expanded states (ρ/ρ₀ ≤ 1 and > _cv , where _cv is the energy of complete vaporization), and a transition region between the hot and cold expanded states (where ρ/ρ₀ < 1 and _iv < < _cv ). First, for compressed states and for cold expanded states the pressure is given by

$$\begin{eqnarray}&&p=\left[a+\displaystyle \frac{b}{\tfrac{\epsilon }{{\epsilon }_{0}{\eta }^{2}}+1}\right]\rho \epsilon +A\mu +B{\mu }^{2},\end{eqnarray} \tag{ A32 }$$

where a, b, A, B, and ₀ are fitted parameters, η = ρ/ρ₀, μ = η − 1, and ρ₀ is the reference density of the material. In order to ensure numerical stability, a low-density cutoff is applied in the cold expanded state when ρ/ρ₀ < 0.8, where a nominal pressure of 10 Pa is given to a zone. In hot expanded states,

$$\begin{eqnarray}&&\begin{array}{l}p=a\rho \epsilon \\ \ +\,\left[\displaystyle \frac{b\rho \epsilon }{\tfrac{\epsilon }{{\epsilon }_{0}{\eta }^{2}}+1}+A\mu \exp \left\{-\beta \left(\displaystyle \frac{1}{\eta }-1\right)\right\}\right]\\ \ \exp \left\{-\alpha {\left(\displaystyle \frac{1}{\eta }-1\right)}^{2}\right\},\end{array}\end{eqnarray} \tag{ A33 }$$

where α and β are two more fitted parameters that control the rate at which the EOS tends to the ideal gas law. For expanded states with internal energies between that of incipient and complete vaporization a hybrid formula is used to transition smoothly between the cold and hot expanded states,

$$\begin{eqnarray}&&p=\displaystyle \frac{\left(\epsilon -{\epsilon }_{{iv}}\right){p}_{E}+\left({\epsilon }_{{cv}}-\epsilon \right){p}_{C}}{{\epsilon }_{{cv}}-{\epsilon }_{{iv}}},\end{eqnarray} \tag{ A34 }$$

where p_C is the pressure calculated using Equation (A32) and p_E is the pressure calculated using Equation (A33).

Using the Tillotson EOS the energy equation can be solved analytically. However, identifying the relevant regime to use in our code requires knowledge of the solution. We therefore find the possible solutions for in all four regimes and use the physically meaningful solution. Subbing Equations (A32), (A33), and (A34) into the energy equation in finite-difference form (Equation (A12)) produces three separate quadratic equations. These can then be solved analytically giving two solutions for each equation. In the case that the material is in compression, we choose the solution of the equation derived from Equation (A32) that (i) gives the correct sign of the change in for the density change, and (ii) is closest to the previous value of . If the material is in an expanded state, the solutions from all three equations are potentially valid and so we choose the solution that (i) is in the correct energy range for its regime, (ii) gives the correct sign of the change in for the density change, and (iii) is closest to the previous value of . The relevant equation can then be used to find the pressure. Note that care must be taken in averaging properties in the cross-over region.

In order to calculate artificial viscosity we also need to calculate the adiabatic sound speed. The sound speed is defined as

$$\begin{eqnarray}&&{c}^{2}={\left.\displaystyle \frac{\partial p}{\partial \rho }\right|}_{S},\end{eqnarray} \tag{ A35 }$$

where S is specific entropy. Using standard differential relations, this can be expressed as

$$\begin{eqnarray}&&{\left.{c}^{2}=\displaystyle \frac{\partial p}{\partial \rho }\right|}_{\epsilon }+{\left.\displaystyle \frac{\partial p}{\partial \epsilon }\right|}_{\rho }{\left.\displaystyle \frac{\partial \epsilon }{\partial \rho }\right|}_{S}.\end{eqnarray} \tag{ A36 }$$

Substituting for the final term from the fundamental thermodynamic relation,

$$\begin{eqnarray}&&{\left.{c}^{2}=\displaystyle \frac{\partial p}{\partial \rho }\right|}_{\epsilon }+\displaystyle \frac{p}{{\rho }^{2}}{\left.\displaystyle \frac{\partial p}{\partial \epsilon }\right|}_{\rho }.\end{eqnarray} \tag{ A37 }$$

All the terms in this equation are known or can be found by differentiating the relevant equation for pressure in each regime. Equation (A37) is used to calculate the sound speed in each regime with the exception of the low-density cutoff where a nominal sound speed of 0.5u_G is used.

Due to the limitations of the Tillotson EOS, WONDY can encounter unphysical solutions in expanded states. In such cases the code is unable to continue, and most of our calculations using Tillotson did not complete the full model run. However, unphysical states are usually encountered late in the simulation when the loss has plateaued, and so the total loss calculated is only slightly underestimated.

The breakout of the shock at the base of the atmosphere/ocean is simulated by giving the lower boundary an initial velocity, u_G, and then allowing the boundary to evolve ballistically. The lower boundary evolves following the momentum equation (Equation (A5)) ignoring forces other than gravity, which in finite-difference form is expressed as

$$\begin{eqnarray}&&{a}_{0}^{n}=-\displaystyle \frac{{{GM}}_{p}}{{\left({x}_{0}^{n}\right)}^{2}}.\end{eqnarray} \tag{ A38 }$$

In this framework, u_G is the particle velocity of the surface upon release of the shock in the planet to the base of the atmosphere or ocean. As for the case of the ocean–atmosphere interface, u_G can be determined using an impedance-match calculation. Due to the lower impedance of gases as opposed to water, release of the shock directly to an atmosphere results in higher ground velocities than in the presence of an ocean, and in both cases the ground velocity is lower than if the shock released to vacuum. When calculating loss by convolving 1D results with the ground velocity calculated in 2D or 3D impact simulations it is imperative that corrections are made to translate the ground velocity.

Driving the lower boundary in this way creates some numerical complications. The sudden acceleration of the boundary generates excessively high artificial viscosities in the lowermost zones. This causes these zones to have artificially high temperatures and low densities which can be outside the range of the EOS. We overcome this by limiting the artificial viscosity in the first four zones to be less than q_max. The exact value of q_max, within a reasonable range, has little effect on the amount of loss so we set it to 0.2 of that calculated for the first time step for the first zone using Equation (A17). We also initialize the first zone with ${q}_{1}^{0}={q}_{\max }$ and the second and third zones with ${q}_{j-\tfrac{1}{2}}^{0}=0.01{q}_{j-\tfrac{3}{2}}^{0}$. This helps to maintain reasonable values for variables in the lowermost zones with no effect on the final results.

In some simulations (see Section 6) the velocity of the ground is forced to increase linearly over some finite rise time, t_rise, until the specified ground velocity is reached. The specification of the artificial viscosity in cells just above the boundary is the same as in the standard case.

The lower boundary is assumed to stop when it returns to its initial position. However, it is necessary to slow the boundary gradually to avoid numerical errors. The boundary is prescribed to fall exponentially after it returns to a position ${x}_{0}\lt {R}_{t}+0.15\left({x}_{0}^{\max }-{R}_{t}\right)$ where x₀ ^max is the maximum height reached by the boundary. The rate is chosen so that the velocity of the boundary is continuous from before being slowed. In finite-difference form,

$$\begin{eqnarray}&&{u}_{0}^{n+\tfrac{1}{2}}={u}_{0}^{n-\tfrac{1}{2}}\exp \left\{\displaystyle \frac{{u}_{\mathrm{slow}}{\rm{\Delta }}{t}^{n-\tfrac{1}{2}}}{0.15\left({x}_{0}^{\max }-{R}_{t}\right)}\right\},\end{eqnarray} \tag{ A39 }$$

where u_slow is the velocity of the boundary at the point the boundary starts to be slowed. Since the slowdown happens late in the simulation the time step can be quite large by the time the slowdown begins. We therefore artificially reduce the time step when the boundary is ${R}_{t}+0.15\left({x}_{0}^{\max }-{R}_{t}\right)\lt {x}_{0}\lt {R}_{t}\,+0.3\left({x}_{0}^{\max }-{R}_{t}\right)$ to ≤0.1 s so that the slowdown can be properly calculated. After slowdown begins the time step is allowed to evolve normally.

In some simulations (see Section 6) the ground is stopped at some time, t_stop. This is achieved by immediately setting the velocity of the boundary to zero at the specified time. Once the boundary has stopped, it begins to fall under gravity and its evolution continues as described above.

The boundary at the top of the atmosphere is treated as a stress-free boundary, i.e., a vacuum. This is implemented by having an additional zone at the top of the atmosphere with (ρ, p) = 0 at all time steps. The finite-difference equations are then applied as normal.

The initial conditions for our simulations are a hydrostatic atmosphere and, optionally, an ocean. We set the extent of the atmosphere by setting the pressure and temperature at the base of the atmosphere. The atmosphere is then assumed to be polytropic,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\tfrac{\gamma -1}{\gamma }}={\left(\displaystyle \frac{\rho }{{\rho }_{0}}\right)}^{\gamma -1}=\displaystyle \frac{T}{{T}_{0}},\end{eqnarray} \tag{ A40 }$$

where p₀, ρ₀, and T₀ are the pressure, density, and temperature at the base of the atmosphere. We then integrate up from the base of the atmosphere using a 4th-order Runge–Kutta method to find the top of the atmosphere (in this model the pressure drops to zero at a finite height). The atmosphere is then divided into the specified number of zones, N_atm, each of equal thickness. For the ideal gas EOS the structure of the atmosphere is determined by p₀, T₀, γ, and m_a (the molar mass of the gas), which sets the density at a given p, T. The values for γ and m_a for each of the atmospheres used here are given in Table 2.

The extent of the ocean is set by specifying the depth of the ocean, H_oc. The temperature and pressure of the bottom of the atmosphere are assumed to be continuous with the ocean. The exact structure of the ocean, however, depends on the EOS being used. For the tabulated EOS an isothermal ocean is initialized with temperature T₀. In the case of the IAPWS EOS an isentropic ocean is used with the entropy set by p₀ and T₀. In order to compare directly with the work of Genda & Abe (2005), when using the Tillotson EOS we initialize a constant specific internal energy ocean with = 120 J kg⁻¹. In each case a 4th-order Runge–Kutta method is used to integrate down from the surface to find the properties for each of the N_oc ocean cells. Due to the limited compressibility of water, the ocean is close to isothermal in both the isentropic and isoenergetic initializations. Cells are of equal thickness.

The masses, radii, and escape velocities of the planets used in this study are given in Table 3. For Earth- and Mars-mass planets the radii were taken as the present-day equatorial radii of those planets. For intermediate-mass planets, the radii were calculated using the HERCULES planetary structure code (Lock & Stewart 2017; Lock 2019) using the present-day Earth's core mass fraction of 0.323 (Yoder 1995). In HERCULES, a body is modeled as consisting of a series of nested concentric layers of constant density and a potential field method is used to calculate the equilibrium structure of the body with a given thermal state, composition, mass, and angular momentum using a realistic EOS. For this work, the mantle and core of each planet were assumed to be forsterite and pure iron, respectively, and were modeled using the the ANEOS EOS model (Thompson & Lauson 1972; Melosh 2007; Canup 2012). The EOSs are documented as "aneos-T70" (iron) and "aneos-gadget" (forsterite), respectively, in two Zenodo repositories (Stewart et al. 2019; Stewart 2020). The mantle was assumed to be isentropic with a specific entropy of 3.2 kJ K⁻¹ kg⁻¹, corresponding to a mantle potential temperature of around 1900 K. This thermal state approximates that of the Hadean mantle. The core was also assumed to be isentropic and have a thermal state similar to the present day. The core specific entropy was set at 1.5 kJ K⁻¹ kg⁻¹, corresponding to a temperature of 3800 K at the pressure of the present-day core–mantle boundary. We used the same HERCULES parameters as in previous studies, which have been shown to accurately model the structure of Earth-like planets (Lock & Stewart 2017, 2019; Lock 2019; Lock & Stewart 2020). The planets were modeled by ${N}_{\mathrm{lay}}^{\mathrm{core}}=20$ evenly spaced layers in the core and ${N}_{\mathrm{lay}}^{\mathrm{mantle}}=80$ layers in the mantle, with N_μ = 1000 points describing the shape of each layer. The expression for the gravity field was truncated at order $2{k}_{\max }=12$. The minimum pressure at the surface of the planet was set to 10 bar. The tolerance for the convergence of the shape of equipotential layers was ${\xi }_{\mathrm{toll}}^{\mu }={10}^{-10}$ and the tolerance for the convergence of the mass of the planet was ξ_toll = 10⁻⁸. The step used for calculating gradients in the solution algorithm was δ ξ = 10⁻². For further details of the definitions of these parameters, the reader is referred to the HERCULES user manual (Lock 2019).

We have tested the sensitivity of our 1D models to the intrinsic code parameters (B₁, B₂, ${K}_{{{\rm{t}}}_{1}}$, and ${K}_{{{\rm{t}}}_{2}}$) as well as the number of zones in the ocean and atmosphere (N_atm and N_oc). To do this, we ran a series of calculations varying each of the parameters in turn and comparing the calculated loss. Default paramters are shown in Table 4. The final results showed little variation with any of the parameters over the ranges we explored (e.g., Figure 16), although some runs with lower B₁ values failed. For some combinations of parameters ringing was observed around the shock front, indicating insufficient numerical viscosity, but such effects did not affect the final loss. We are therefore confident that our results are not significantly affected by our choice of these parameters.

Zoom In Zoom Out Reset image size

Table 4. Values for Constants Used in All 1D Model Runs for Which Results Are Presented Here (Details in Text)

Parameter	Standard Value
B1	2
B2	0.1
${K}_{{{\rm{t}}}_{1}}$	0.75
${K}_{{{\rm{t}}}_{2}}$	1.05
Noc	500
Natm	500
tfinal	5 × 104 s

Notes. B₁ and B₂ are used to determine the strength of the artificial viscosity. ${K}_{{{\rm{t}}}_{1}}$ and ${K}_{{{\rm{t}}}_{2}}$ control the size of each time step. N_oc and N_atm are the number of zones used for the ocean and atmosphere, respectively. t_final is the total runtime for the simulations.

Download table as:  ASCIITypeset image

For the no-ocean case, we describe the loss using a modified logistics function,

$$\begin{eqnarray}&&{f}_{\mathrm{NO}}\left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right)={\alpha }_{1}{\left[1+\exp \left\{{\alpha }_{2}\left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right)+{\alpha }_{3}\right\}\right]}^{{\alpha }_{4}}+{\alpha }_{5},\end{eqnarray} \tag{ C1 }$$

where u_G is the ground velocity, v_esc is the escape velocity, and α_i are constants. By requiring loss to be zero in the case of zero ground velocity, i.e., f_NO = 0 when u_G = 0, we find

$$\begin{eqnarray}&&{\alpha }_{5}=-{\alpha }_{1}{[1+\exp ({\alpha }_{3})]}^{{\alpha }_{4}}.\end{eqnarray} \tag{ C2 }$$

Similarly, by requiring loss to be complete when u_G = v_esc, we find

$$\begin{eqnarray}&&{\alpha }_{1}=\displaystyle \frac{1}{{[1+\exp ({\alpha }_{2}+{\alpha }_{3})]}^{{\alpha }_{4}}-{[1+\exp ({\alpha }_{3})]}^{{\alpha }_{4}}}.\end{eqnarray} \tag{ C3 }$$

In addition, we force f_NO = 0 when u_G/v_esc < 0 and f_NO = 1 when u_G/v_esc > 1.

In the no-ocean case, we set α₃ = −1 and performed a least-squares fit to find α₂ = −4.32 and α₄ = −3.77. For the fit we used the calculations shown in Figure 5(A): 0.1, 0.5, 1, 5, 10, 50, 100, and 500 bar atmospheres on planets with mass 0.107 (M_Mars), 0.3, 0.5, 0.7, 0.9, and 1 M_Earth. The surface temperature was 300 K, and the atmosphere was CO₂ with a mean molecular weight of m_a = 44 and a ratio of specific heat capacities of γ = 1.29. Loss was calculated at 0.05 v_esc increments from 0.05 to 0.95 v_esc.

This parameterization offers a very accurate fit to our model runs. Figure 5(C) shows the misfit between our parameterization and the fitted calculations. The rms misfit for all runs is 0.0046 and the maximum misfit is 0.0089.

We chose to parameterize loss in the ${ \mathcal R }$ - u_G–loss space (Figure 10). Loss is described by a modified logisitics function,

$$\begin{eqnarray}&&f({u}_{{\rm{G}}},{v}_{\mathrm{esc}},{ \mathcal R })={\alpha }_{1}{[1+\exp ({\alpha }_{2}{\mathrm{log}}_{10}\left({ \mathcal R }\right)+{\alpha }_{3})]}^{{\alpha }_{4}}+{\alpha }_{5},\end{eqnarray} \tag{ C4 }$$

where α_i are functions of the ground velocity and the escape velocity, v_esc. We find that loss due to a given ground velocity tends to the no-ocean case as ${ \mathcal R }\to \infty$ (Figure 10) and so we enforce f → f_NO as ${ \mathcal R }\to \infty$ by setting

$$\begin{eqnarray}&&{\alpha }_{1}={f}_{\mathrm{NO}}-{\alpha }_{5},\end{eqnarray} \tag{ C5 }$$

and requiring that α₂ ≤ 0.

The velocity dependences of α_i are given by an arbitrary function chosen as they provided a good fit to our simulation results. α₂ and α₃ are described by 3rd-order polynomials:

$$\begin{eqnarray}&&{\alpha }_{i}={a}_{1}^{i}+{a}_{2}^{i}\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}+{a}_{3}^{i}{\left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right)}^{2}+{a}_{4}^{i}{\left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right)}^{3},\end{eqnarray} \tag{ C6 }$$

and α₄ is given by

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{4} & = & {a}_{1}^{4}\exp \left\{\sin \left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\pi \right)\right.\\ & & \left.\times \,\left[{a}_{2}^{4}\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}+{a}_{3}^{4}{\left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right)}^{2}+{a}_{4}^{4}{\left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right)}^{3}\right]\right\},\end{array}\end{eqnarray} \tag{ C7 }$$

where aⁱ _j are coefficients that can depend on v_esc. The functional form for the ground velocity dependence of α₅ is more complex. The value of α₅ is the asymptote for loss as ${ \mathcal R }\to 0$ and is close to the lowest ${ \mathcal R }$ loss curve we calculate (e.g., the yellow line in Figure 6). To describe a loss curve in u_G–loss space, we divide the functional form for α₅ into two parts for the atmospheric-loss (g₁) and ocean-loss (g₂) regimes, respectively:

$$\begin{eqnarray}{\alpha }_{5}=\left\{\begin{array}{ll}0 & {u}_{{\rm{G}}}\leqslant 0\\ {g}_{1} & 0\lt {u}_{{\rm{G}}}\lt {u}_{\mathrm{tal}}\\ 1 & {u}_{\mathrm{tal}}\leqslant {u}_{{\rm{G}}}\lt {v}_{\mathrm{esc}}\ \&amp;\ {g}_{2}\leqslant 1\\ {g}_{2} & {u}_{\mathrm{tal}}\leqslant {u}_{{\rm{G}}}\lt {v}_{\mathrm{esc}}\ \&amp;\ {g}_{2}\gt 1\\ 2 & {u}_{{\rm{G}}}\geqslant {v}_{\mathrm{esc}},\end{array}\right.\end{eqnarray} \tag{ C8 }$$

where

$$\begin{eqnarray}&&{g}_{1}={a}_{1}^{5}{\left(\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right)}^{{a}_{3}^{5}}\exp \left({a}_{2}^{5}\displaystyle \frac{{u}_{{\rm{G}}}}{{v}_{\mathrm{esc}}}\right),\end{eqnarray} \tag{ C9 }$$

à

$$\begin{eqnarray}\begin{array}{rcl}{g}_{2} & = & 2+{a}_{4}^{5}\left[\displaystyle \frac{{u}_{{\rm{G}}}-{u}_{\mathrm{tal}}}{{v}_{\mathrm{esc}}-{u}_{\mathrm{tal}}}-1\right]\\ & & +\,{a}_{5}^{5}\left[{\left(\displaystyle \frac{{u}_{{\rm{G}}}-{u}_{\mathrm{tal}}}{{v}_{\mathrm{esc}}-{u}_{\mathrm{tal}}}\right)}^{2}-1\right],\end{array}\end{eqnarray} \tag{ C10 }$$

and u_tal is the velocity at which total atmospheric loss is achieved (i.e., when g₁ = 1), given by

$$\begin{eqnarray}&&{u}_{\mathrm{tal}}={v}_{\mathrm{esc}}\left(\displaystyle \frac{{a}_{3}^{5}}{{a}_{2}^{5}}\right)\mathop{\min }\limits_{k=-1,0}{W}_{k}\left\{\displaystyle \frac{{a}_{3}^{5}}{{a}_{2}^{5}}{\left({a}_{1}^{5}\right)}^{-1/{a}_{3}^{5}}\right\},\end{eqnarray} \tag{ C11 }$$

where W_k is the Lambert W function of order k. The velocity regime between the atmospheric and ocean-loss regimes where α₅ = 1 loss emulates a feature seen in our low-${ \mathcal R }$ simulations where significant ocean loss is delayed for a short range of u_G after total atmospheric loss is achieved (Figure 6).

Dependence on planetary mass is captured by a linear dependence of the coefficients for α₂ (a² _j ), α₄ (a⁴ _j ), and α₅ (a⁵ _j ) on v_esc such that

$$\begin{eqnarray}&&{a}_{j}^{i}={a}_{j,1}^{i}+{a}_{j,2}^{i}\displaystyle \frac{{v}_{\mathrm{esc}}}{{v}_{\mathrm{esc}}^{\mathrm{Earth}}},\end{eqnarray} \tag{ C12 }$$

where ${v}_{\mathrm{esc}}^{\mathrm{Earth}}=11.2$ km s⁻¹ is the escape velocity of Earth. Note that it is possible to parameterize the v_esc dependence on α₅ by using a function for the absolute value of u_G convolved with the no-ocean loss function, but we found that, because of the complications of the varying pressure of release and the highly nonlinear no-ocean loss function, a simple escape velocity scaling of the parameters was able to give a more accurate fit.

To determine the parameters ${a}_{j,n}^{i}$ we performed a least-squares fit on the results of the simulations shown in Figure 10: all combinations of 1, 5, 10, 50, 100, 300, and 500 bar atmospheres and oceans of depths 0.1, 0.5, 1, 2, 3, 5, 10, 20, and 30 km, on planets with mass M_Mars (0.107 M_Earth) and 0.3, 0.5, 0.7, 0.9, and 1 M_Earth. Additional simulations with 900 bar atmospheres and oceans of 0.1 km depth were also included. The surface temperature was assumed to be 300 K, and the atmosphere was CO₂ (m_a = 44, γ = 1.29). Loss was calculated at 0.05 v_esc increments from 0.05 to 0.95 v_esc. The best-fit parameters are given in Table 5.

Table 5. Fitted Coefficients for Each of the Parameters ${a}_{j,k}^{i}$ (see Equation (C12))

i	j	${a}_{j,1}^{i}$	${a}_{j,2}^{i}$
2	1	−16.2584384	14.2828078
2	2	71.5029708	−86.4465281
2	3	−117.862430	164.871407
2	4	60.1395645	−95.6142191
3	1	−7.66878923	⋯
3	2	43.9524753	⋯
3	3	−105.616622	⋯
3	4	77.4557563	⋯
4	1	−0.767532594	−0.263620386
4	2	−12.0443695	−11.5024085
4	3	44.7091911e	39.8696345
4	4	−38.8346991	−28.7008315
5	1	10660.1668	13897.9517
5	2	−12.9490638	0.640100268
5	3	7.35306571	−0.384763832
5	4	2.12776033	−0.324256357
5	5	−1.14723912	0.418100341

Download table as:  ASCIITypeset image

We find that our parameterization offers an accurate representation of the dependence of loss on ${ \mathcal R }$, u_G, and M_p. Figure 10 shows the fit in ${ \mathcal R }$-loss space for example ground velocities. The global rms is 0.041, the rms at any given ground velocity does not exceed 0.062, and the maximum misfit of any simulation result is 0.25. Unsurprisingly, the misfit is largest in the transition between the low- and high-${ \mathcal R }$ regimes where the gradient of the function in the ${ \mathcal R }$–loss regime is greatest and the slight sensitivity to parameters such as absolute velocity, atmospheric pressure, ocean height, etc., is greatest.

The parameterization also reproduces realistic loss curves in u_G–loss space. Figure 17 shows our simulation results for the same four example loss curves as in Figure 6(A) but for CO₂ atmospheres (solid lines), with corresponding curves calculated using our parameterization (dashed lines). The gray band gives the full range of possible loss for an Earth-mass planet for any value of ${ \mathcal R }$ using our parameterization.

Zoom In Zoom Out Reset image size

As discussed in Section 5, the relationship between the strength of the shock in the planet and the ground velocity varies depending on the properties of the ocean and/or atmosphere in contact with the ground. It is also worth noting that the relationship could be affected somewhat by changing the shock properties of the ground. When calculating the loss due to a given shock, for example when combining 1D and 3D simulations (Section 7.4), it is therefore necessary to calculate the relationship between the properties of the shock in the planet (here we use the particle velocity, u_p, as the reference variable) and the ground velocity for the specific combination of ground and ocean and/or atmosphere that is under consideration. As this relationship can vary substantially, we do not attempt to provide a comprehensive set of u_p–u_G relations. Instead, we provide a Python package and a widget that workers can use to calculate individual or sets of u_p–u_G relations (see data availability) for different surface conditions and a range of EOSs.

We do make one exception. The u_p–u_G relationship for breakout of a shock into the base of the ocean from a forsterite mantle is relatively constant over the range of conditions we have considered (oceans of depths 0.1–30 km, and atmospheric pressure up to 900 bar). In this case, the u_p–u_G relationship can be well described by

$$\begin{eqnarray}&&{u}_{{\rm{G}}}={\beta }_{1}{u}_{p}^{{\beta }_{2}},\end{eqnarray} \tag{ C13 }$$

where β₁ = 2.6441 and β₂ = 0.9252.