SPH calculations of Mars-scale collisions: the role of the Equation of State, material rheologies, and numerical effects
SSPH calculations of Mars-scale collisions: the role of the Equationof State, material rheologies, and numerical effects
Alexandre Emsenhuber a, ∗ , Martin Jutzi a , Willy Benz a a Physikalisches Institut & Center for Space and Habitability, Universität Bern, Gesellschaftsstrasse 6,CH-3012 Bern, Switzerland
Abstract
We model large-scale ( ≈ ) impacts on a Mars-like planet using a Smoothed ParticleHydrodynamics code. The effects of material strength and of using different Equationsof State on the post-impact material and temperature distributions are investigated. Theproperties of the ejected material in terms of escaping and disc mass are analysed as well.We also study potential numerical effects in the context of density discontinuities and rigidbody rotation. We find that in the large-scale collision regime considered here (with impactvelocities of / s ), the effect of material strength is substantial for the post-impact dis-tribution of the temperature and the impactor material, while the influence of the Equationof State is more subtle and present only at very high temperatures. Keywords:
Terrestrial planets, impact processes, interiors
1. Introduction
Giant impacts occurring at the end stages of planet formation define the properties of thefinal planets and moons. Examples include Mercury’s anomalously thin silicate mantle (Benzet al., 1988, 2007; Asphaug and Reufer, 2014), the origin of the Earth’s Moon (Hartmann andDavis, 1975; Cameron and Ward, 1976; Canup and Asphaug, 2001; Canup, 2004, 2012; Ćukand Stewart, 2012; Reufer et al., 2012) or the formation of the Pluto-Charon system (Canup,2005, 2011). Smaller, but still planet-scale collisions were proposed for the formation of themartian dichotomy (Wilhelms and Squyres, 1984; Frey and Schultz, 1988; Andrews-Hannaet al., 2008; Nimmo et al., 2008; Marinova et al., 2008, 2011). Numerical simulations of suchgiant impact scenarios are performed using shock physics codes, which are often based onthe smoothed particle hydrodynamics (SPH) method. In such simulations of planet-scalecollisions it is typically assumed that material strength is negligible due to large overburdenpressures caused by self-gravitation. However, it was found in previous studies (e.g. Jutzi,2015) that for collisions at
100 km to scale, the effect of material strength can still be ∗ Corresponding author
Email addresses: [email protected] (Alexandre Emsenhuber), [email protected] (Martin Jutzi), [email protected] (Willy Benz) a r X i v : . [ a s t r o - ph . E P ] O c t ery important. The impactor size range for which the assumption of purely hydrodynamicbehaviour is justified has not yet been studied systematically. Other uncertainties in giantimpact simulations are related to the Equation of State (EoS) models. Often used EoSare the simple Tillotson (Tillotson, 1962) EoS (e.g. Marinova et al., 2008, 2011) or the moresophisticated (M-)ANEOS (Thompson and Lauson, 1972; Melosh, 2007). Last but not least,there are potential numerical issues that have to be added to the uncertainties. Concerningthe SPH method, known numerical issues are related to surface and contact discontinuities(e.g. Hosono et al., 2013; Reinhardt and Stadel, 2017) or rotational instabilities in rigid bodyrotations (e.g. Speith, 2006).In this study, we focus on large-scale ( ≈ ) impacts on a Mars-like planet. Theconditions studied here are chosen to be quite general, but at the same time they also covera range of possible dichotomy forming events (Marinova et al., 2008). We want to assess therelative importance of the effects mentioned above (i.e., material strength; choice of EoS;numerical issues) in collision simulations of such scale. This knowledge is important in orderbe able to make an educated choice of the material models and numerical methods, and toknow their limitations.Collisions are modelled using an updated SPH code (Reufer et al., 2012), which includesself-gravity, a newly implemented strength model (following Jutzi, 2015) and various EoS(Tilloston and M-ANEOS). We investigate the effects of the above mentioned material prop-erties (namely material strength and the EoS models) on the outcome of the collisions fordifferent impact geometries. We focus on the material and temperature distributions in thefinal planet but also analyse the properties of the ejected material (orbiting and unbound).Finally, potential numerical effects are studied as well. We investigate the known SPH issuein the case of rigid body rotation. We also consider different schemes to compute the density,in order to investigate potential numerical issues at discontinuities (such as the free surfaceand the core-mantle boundary).A subset of the SPH calculations presented here is coupled with thermochemical simula-tions of planetary interiors in a companion paper (Golabek et al., 2017). The effects of thedifferent models on the long-term evolution (over a time period of about . ) and thecrust formation will be discussed there.In section 2 we describe our modeling approach and discuss in detail the implementedstrength model and its coupling with the EoS models. The initial conditions (target andprojectile properties and impact geometry) are detailed in section 3. In section 4, we presentthe results of our simulations using different material models and numerical schemes, followedby a discussion and conclusions in section 5.
2. Modelling approach
We use a smoothed particle hydrodynamics (SPH) code (e.g. Benz and Asphaug, 1994;Reufer et al., 2012; Jutzi, 2015) to model the impact event. The code used here is based onthe SPHLATCH code developed by Reufer et al. (2012). It includes a newly implementedpressure and temperature dependent shear strength model, as described below, which isappropriate for the large scale collisions considered here. To study smaller scale collisions2NEOS Tillotsoniron silicate iron olivineG [GPa] 76 44.4 76 72 σ M [GPa] 0.68 3.5 0.68 3.5 u m [J/kg] - - . · . · T m [K] - - 920 3400 Table 1: Material parameters for the solid rheology model. G is the shear modulus. σ M is the von Misesplastic limit. Melting temperature T m is computed from u m using equation (11). in geological materials, we use a different code, which also includes a tensile fracture and aporosity model (see Jutzi, 2015, for details).SPH uses a Lagrangian representation where material is divided into particles. Quantitiesare interpolated (‘smoothed’) by summing over surrounding particles (called neighbours)according to B ( (cid:126)x ) = (cid:88) i B i W ( (cid:126)x − (cid:126)x i , h i ) V i (1)where B i represents the quantity (field variable) to be interpolated, and (cid:126)x i , h i and V i are theposition, smoothing length, and volume of particle i, respectively. W ( (cid:126)x, h ) is a smoothingkernel, which has the propriety to vanish if (cid:107) (cid:126)x (cid:107) > h , so that the hydrodynamic variables(pressure, density, stress tensor) are integrated over a local group of neighbouring particles.In our simulations, we use a 3D cubic spline kernel (Monaghan and Lattanzio, 1985).This interpolation scheme is used in SPH to solve the relevant differential equations. Instandard SPH, the density is computed by direct summation (equation 1). We refer to thismethod as density summation through this article. Alternatively, it can be computed byusing the changing rate of density, which is given by the continuity equation: DρDt + ρ ∂v x ∂x + ρ ∂v y ∂y + ρ ∂v z ∂z = 0 (2)where D/Dt stands for the substantive time derivative. We call it density integration . Asimilar equation gives the time derivative of the smoothing length:
DhDt + h ∂v x ∂x + h ∂v y ∂y + h ∂v z ∂z = 0 (3)Thus the smoothing length grows in sparse regions and shrinks in dense regions, so that thenumber of neighbours within h is ≈ .Body accelerations are the result of the pressure gradient, and for this we use the pressureas computed from an equation of state (EoS) (see below), which is a function P ( ρ, u ) ofdensity ρ (see equation 2) and internal energy u . For solid materials the pressure gradientis generalised into a stress tensor (Benz and Asphaug, 1994) with the resulting equation ofmotion given by: Dv a Dt = 1 ρ ∂τ ab ∂x b (4)3he stress tensor is defined as τ ab = − P δ ab + σ ab (5)where P is the pressure, δ ab is the Kronecker symbol and σ ab is the traceless deviatoricstress tensor. The time evolution of σ ab is computed adopting Hooke’s law as in Benz andAsphaug (1994, 1995). In the limit of zero deviatoric stress this reduces to the familiarpressure-gradient acceleration. Energy conservation now reads DuDt = 1 ρ τ αβ ˙ ε αβ (6)with ˙ ε αβ being the strain rate tensor.In standard SPH, the computation of the strain and rotation rate tensors fails to conserveangular momentum in rigid body rotations. Following the approach by Speith (2006), we usea correction tensor in the computation of the velocity derivatives, which allows conservationof angular momentum at the cost of an additional computation step.Equations (4) and (5) describe an entirely elastic material. To model plastic behaviourwe use a Drucker-Prager-like yield criterion (e.g. Collins et al., 2004; Jutzi, 2015). The modelhas one yield strength σ i , corresponding to intact material, and another σ d correspondingto completely fragmented material: σ i = C + µ i P µ i P/ ( σ M − C ) (7) σ d = µ d P (8)where C is the cohesion (yield strength as zero pressure), σ M is the von Mises plastic limit,and µ i and µ d are the coefficients of friction for intact and completely damaged material,respectively. We set them as µ i = 2 and µ d = 0 . as in Collins et al. (2004). Note that σ d islimited to σ d < σ i at high pressures. For the large-scale collisions considered here cohesionand tensile strength are negligible due to the large gravity-induced lithostatic stress (see e.g.Jutzi et al., 2015). We therefore set C = 0 . However, shear strength (limited by σ M ) is stillimportant and cannot be neglected, as we shall see below.The yield strength of intact material is a function of temperature. To take this intoaccount, we adopt the same relation as in Collins et al. (2004): σ → σ tanh (cid:18) ξ (cid:18) T m T − (cid:19)(cid:19) (9)where ξ is a material constant which we set to ξ = 1 . , T is the temperature and T m is thecorresponding liquidus temperature (see below). If T > T m the material is molten and weset σ = 0 . Finally, if the measure of the stress state, the second invariant of the deviatoricstress tensor σ ll = (cid:115) σ ab σ ab (10)exceeds σ , the components of the deviatoric stress tensor are reduced by a factor σ/σ ll . Theparameters for the solid rheology model are provided in table 1.4ame Description Value ρ Reference density / m T Reference temperature p Pressure at the reference point B Bulk modulus at the reference point . · Pa T Debye
Reference Debye temperature
464 K γ max Limiting value of the Gruneisen coefficient / E sep Zero temperature separation energy · J / kg T melt Melting temperature H f Heat of fusion . · J / kg ρ liq /ρ solid Ratio of liquid to solid density at melt point . Table 2: Material parameters for iron with ANEOS equation of state. Other parameters are set to zero.
To complete the set of equations, an equation of state providing a relation betweenpressure, temperature and density is needed. We use either Tillotson (Tillotson, 1962) or(M-)ANEOS (Thompson and Lauson, 1972; Melosh, 2007). The Tillotson equation of stateprovides both the pressure P and the speed of sound as output, but it does not give thetemperature directly, so as an approximation we use the internal energy as a proxy fortemperature by dividing by the heat capacity T = u/c p (11)In this case equation (9) is modified to use u and u m instead of T and T m respectively,and u m is treated as a material constant independent of pressure. Iron parameters are fromTillotson (1962) (also in Melosh, 1989) and the ones for olivine from Marinova et al. (2008).For a more physical computation of temperature, and to handle phase transitions in amore thermodynamically consistent way than is done with Tillotson EoS, we use ANEOS(Thompson and Lauson, 1972) for iron with the parameters provided in table 2, and M-ANEOS (Melosh, 2007) for silicates. These equations of state provide numerous outputvariables including the temperature and phase information (e.g. melt and vapor fraction).Depending on the parameters used in the equation of state, this phase information can begiven in different ways. For iron, it can be used to infer the melting temperature at a givenpressure. However, this is not the case for silicates, for which we apply the same procedureas used in Senft and Stewart (2009) to obtain the melting temperature at a given pressure.It is important to point out that neither the Tillotson nor the ANEOS equation of stateinclude the latent heat of melting and therefore, temperatures above the melt temperatureare overestimated.Shocks occur when the impact speed (or particle velocity) exceeds the sound speed in themedium. If not treated properly, this can lead to particle interpenetration and other non-physical effects. We use a standard artificial viscosity term (Monaghan, 1992) that depositsinternal energy and momentum behind the shock in a manner that is consistent with theHugoniot shock relations (c.f. Melosh, 1989) but spread out over several smoothing lengths.5n addition to the body force accelerations given by equation (4), we also compute theaccelerations due to self-gravity, using a tree-code method (Barnes and Hut, 1986).The set of equations described above provides the time dependence of the physical quant-ities, which are then integrated using a prediction-correction scheme (Press, 2002). For amore detailed description of the method see e.g. Jutzi et al. (2015).
3. Initial conditions
We perform a series of SPH collision models, where we use a Mars-mass target body( r = 3400 km ) and an impactor with radius. We use different EoS (ANEOS orTillotson), material rheologies (solid or fluid) and numerical schemes (density summation orintegration), in order to study the effects on the collision outcome and to asses their relativeimportance. For each combination, we use two impact angles: 0° (head-on) or 45° (oblique),which gives a total of 16 simulations. The mutual escape velocity (Asphaug, 2010) definesthe impact speed. As a nominal case, we define the simulation with ANEOS, solid rheology,integrated density and oblique impact angle.The target and impactor both start with a Mars-like internal structure with an iron coreradius half of the body radius, while SiO2 comprises the bodies’ mantles. The use of SiO2to represent Mars’ basaltic mantle is a simplified but reasonable assumption, given availableequation of state information. Note also that the crust of the target is not resolved at theselarge scales and is not explicitly included in the model. For simulations using Tillotson EoS,each layer is initially isothermal; iron core temperature is fixed to T Fe = 1800 K and mantleto T Si = 1500 K . For the simulations with ANEOS, an isentropic profile is used for whichcore’s central temperature is and for the silicate at the core-mantle boundary it is . The choice of these values is explained in Golabek et al. (2017).We begin by setting up self-gravitating, hydrostatically-equilibrated planets, startingwith one-dimensional (1D) spherically symmetric bodies modelled using a Lagrangian hy-drocode (Benz, 1991) with the same EoS as in the SPH model. This structure is divided in1 000 cells (100 for the impactor) of equal mass. Cell boundaries according to the force bal-ance between self-gravity and pressure (including a damping term). This profile is evolveduntil hydrostatic equilibrium is reached so that radial velocities are small (less than 1% ofthe escape velocity). Afterwards we transfer the 1D radial profile onto SPH particles thatare placed onto a 3D lattice. Parameters of each particle are copied from the profile ac-cording to the radius. As particles are equally spaced on the lattice, variation of densityis taken into account by adjusting the particle mass. The spherical SPH bodies are thenalso evolved in an initialising step to reach a hydrostatic equilibrium and negligible radialvelocities. Thus the SPH simulations start with two relaxed differentiated spherical planetsapproaching from several radii away; they get tidally deformed prior to the impact.The SPH simulations are performed with a resolution of about one million particlesfor the target. The number of particles for the impactor is scaled according to the massratio between the two objects, so that particle spacing h is approximately constant, for the6reatest numerical accuracy during the collision. The corresponding smoothing length isthen approximately
60 km for both bodies.
Initial radial profiles of density, pressure and temperature are shown in figures 1 and 2for the impactor and target respectively. For the cases with the Tillotson EoS, the internalenergy is being evolved during the equilibration phase, hence the spread in temperature.The comparison with a case where this value is kept fixed is discussed later in section 4.3.3.The 1D profiles for the two different approaches to compute the density (density sum anddensity integration) are identical as the distinction only happens during the SPH phaseitself.There is a noticeable difference in the density between the cases with a different EoS:iron has higher density with ANEOS and the opposite happens for silicate. The densitycontrast at the core-mantle boundary is thus higher when using ANEOS. Since the objectshave the same mass in the two cases, their radius is affected by this density variationand we obtain about using ANEOS and using Tillotson, which makes arelative difference of the radii around . It is known that standard SPH has problemshandling sharp density changes, at boundaries between materials and at the surface. We seeartificial numerical effects at these locations. With density summation, particles beyond theboundary enter in the SPH sum (equation 1) and thus smooth the discontinuity. Particlesin this region then have a slightly different density. At the surface, the effect is due to thelack of close-by particles (see Reinhardt and Stadel, 2017 for an improved scheme). Thedensity variation is reflected in pressure and temperature. Although in the case of densityintegration there is no such summation involved to compute the density, there are stillsome artificial effects (oscillations) at the boundaries. More sophisticated ways to deal withdiscontinuous boundaries in SPH have been developed recently (e.g. Hosono et al., 2013)and shall be implemented in our code for future studies.However, as we shall see below, the amplitude of the unphysical ’noise’ in the initialprofiles is small enough not to affect significantly the outcome of the collision simulations.
4. Results
We show a time series of the nominal case along with the corresponding fluid case infigure 3. Shown are the post-impact material and temperature distributions. The simu-lations start roughly two hours before the first snapshot. Due to the angular momentumtransferred to the target by the impact, the target begins to rotate after a few hours, hencethe impact location is shifted by ∼ ◦ between the two last snapshots. One can note theimpactor’s tidal deformation before the impact mainly in the fluid case. The effect of thematerial rheology leads to significant differences in terms of the post-impact material andtemperature distribution. Another notable difference is the degree of the impact inducedtarget oscillation, which is much stronger in the fluid case (see snapshots at 3 and ).7 igure 1: Initial profiles of the different impactors used in this study. The 1D profile is depicted by the blackline. Each point is one SPH particle and its color represents the material: red for iron and blue for silicate. igure 2: Initial profiles of the different targets used in this study. T = . h T = . h T = . h T = . h T = . h T = . h T = . h T = . h Figure 3: Time series showing material distribution for the nominal case ( ◦ , ANEOS, solid, integrateddensity; two left columns) and the corresponding fluid case (two right columns). Impactor’s trajectory isclockwise. This is a slice of depth inside the impact plane. Plots are centered on the center of massof the main body. On material plots, colour represents the type and origin: blue is target’s mantle, purpleimpactor’s mantle, red target’s core and yellow impactor’s core. × . . . . . . C u m u l a t i v e m a ss Figure 4: Cumulative mass distribution at the end of the simulation as function of radius. Colour representsmaterial type and origin: red is target’s core, yellow impactor’s core, blue target’s mantle and purpleimpactor’s mantle, as for material plots in figure 3. Solid lines depict the run with solid rheology and dashedlines are for fluid rheology.
The heat generated by the impact is spread out over a broader region with solid rheologycompared to the fluid one. Already at t = 0 .
25 h the higher temperature zone extend toseveral hundreds of kilometres away from the contact zone and expands until covering awide region at the end of our simulations. With fluid rheology only the contact zone hashigh temperature material at the beginning and towards the end only impactor materialretains these high temperatures. Note that the annulus of warm material ( ∼ ) thatforms at the target’s core-mantle boundary with solid rheology is not due to the impact butis an artificial feature caused by the set up routine. The small remaining velocities (on theorder of
10 m / s ) at the beginning of the impact simulation, when forces are enabled andthe damping term removed are sufficient to create instabilities that lead to internal energyincrease due to the artificial viscosity.At early times, the impactor does not penetrate as deep inside the target with solidrheology as in the fluid case (see snapshots at .
50 h and .
92 h ).To get a better overview of the location where impactor material is deposited insidethe target at the end of the simulation, we show a cumulative mass plot against radius infigure 4. The curves are normalised by the total mass of each material.Note that some of the impactor’s silicate material is still orbiting around the body orhas been ejected. This accounts for
13 wt% and
18 wt% for the fluid and solid rheologiesrespectively. In both cases, about one sixth of this material is still bound whereas theremaining is on an escaping orbit.There are significant differences for impactor’s material distribution between the two12heologies.
27 wt% of impactor’s core remains inside the mantle with solid rheology whereasthis value falls to for fluid rheology. The remaining lies at the target’s core-mantleboundary. The impactor’s silicate shows a different behaviour. For fluid rheology, thereis a small part close to the core-mantle boundary, some inside the target’s mantle and themajority close to the surface. For solid rheology, there is one half distributed inside themantle and the other half close to the surface.
The results from our series of 16 simulations, using different material rheologies, EoSand numerical schemes are shown in figure 5 for oblique cases and in figure 6 for head-onones. The differences between density integration and summation will be discussed in section4.3.2. Temperature and energy distributions for four of these cases (oblique impacts; densityintegration) are provided in figure 7. Note that our nominal case and its fluid counterpartthat were discussed in the previous section are shown again in the first row of figure 5(ANEOS; integrated density) and with the black lines on figure 7.At first look, we note that the differences (concerning material and temperature distri-bution) discussed in the previous section between the runs with rheologies are also presentin the simulations using different EoS and numerical techniques (density computation). Themore pronounced temperature increase around the impact zone with solid rheology is no-ticeable on figure 7 between 2000 and ; it is about −
20 wt% for ANEOS and abit less for Tillotson.The differences between the runs with different EoS are more subtle. For instance, usingthe Tillotson EoS, the impactor’s particles close to the surface have higher temperatures.This effect appears both in figures 5 and 6 as well as in 7. In figure 7, the curves for the twoEoS diverge above about , whereas the energy distribution remains similar. ANEOSincludes the latent heat of vaporisation which our simple conversion formula (eq. 11) forTillotson neglects. Hence the high temperatures are overestimated with Tillotson. It shouldbe noted that for silicate, neither EoS take s the latent heat of melting into account. One canalso see in figure 7 that the behaviour at low temperature depends on the initial profile, andthat almost all mantle gets a small temperature increase (of the order of hundred kelvins).When looking at head-on cases, again the same general features appear. Hot surfacematerial is now located directly above the impact location. The fluid cases however showsome specific results. The impact itself produces surface waves that propagate and arefocused at the antipode. At this location, these waves lead to some spallation of mantlematerial, which gets heated when falling back. We thus see an increase in temperature atthe antipode; here there are some notable differences between the Tillotson and the ANEOScases.Comparing the influence of both effect, we find that rheology plays an important rolefor both material and heat distribution, even at this scale. The EoS, despite changing theradius of the bodies by ∼ , has little influence on these effects. Compared to Tillotson, themore sophisticated temperature computation by ANEOS, which includes the latent heat ofvaporisation, leads to differences only at relatively high temperatures (above − ).13olid, material Solid, temperature Fluid, material Fluid, temperature AN E O S d e n s i t y i n t e g r a t i o n AN E O S d e n s i t y s u mm a t i o n T ill o t s o n d e n s i t y i n t e g r a t i o n T ill o t s o n d e n s i t y s u mm a t i o n Figure 5: Material and temperature plots for the oblique impact cases at the end of the simulations. Thetop line is the same as the final one of figure 3. AN E O S d e n s i t y i n t e g r a t i o n AN E O S d e n s i t y s u mm a t i o n T ill o t s o n d e n s i t y i n t e g r a t i o n T ill o t s o n d e n s i t y s u mm a t i o n Figure 6: Same as figure 5, but for the corresponding head-on cases. Initial contact happens on the rightside. T [K]10 − − − − C u m u l a t i v e m a ss u [J / kg]10 − − − − C u m u l a t i v e m a ss Figure 7: Cumulative mass fraction of the mantle with temperature (left panel) or specific internal energy(right panel) greater than a given value. Black lines denote ANEOS whereas red is for Tillotson. Solid linesare for solid rheology, dashed lines for fluid rheology and the dotted lines show the initial target profile.Only density integration with oblique geometry is considered. The dashed grey line on the temperature plotshows the upper bound of the colour scale on figures 3, 5 and 6.
In figure 8 we show the effect of applying the correction tensor in the computation ofthe stress tensor to avoid rotational instability effects in the case of a solid rheology (e.g.Speith, 2006). As it can be seen, the post-impact material and temperature distributionsare very similar in the cases with and without correction tensor, except for the temperatureartefact due to the set up routine at the core-mantle boundary. We note that the amplitudeof the artificial temperature increase at the boundary is much smaller in the simulationwith the correction tensor (about ∼ increase) than in the one without (about ∼ increase). Furthermore, the longer term target rotation is slower in the case withoutcorrection tensor, because the angular momentum is not conserved.The angular momentum is plotted as a function of time in figure 9 for the two cases. Wealso include the corresponding fluid case for comparison. For the solid rheology with correc-tion and the fluid rheology, the variation of angular momentum over the whole simulations(about
20 h ) is lower than whereas for the solid rheology without correction, angularmomentum is reduced by almost
25 % . The decrease starting a few hours after time impactreflects the slow down of the rigid body rotation which is not treated correctly without thecorrection tensor. Ejected material is not subject to this problem as it does not encountersolid forces. The rotation rate decrease implies that target’s position is not the same duringlate stages of simulations. At 18 hours after impact, the target rotated around 20 degreesless than in the case where angular momentum is conserved (figure 8). Here, the targetrotates about ◦ between the merging of both cores at two hours after the impact and theend at 18 hours. 16orrected Non-correctedmaterial temperature material temperature Figure 8: End result comparison for case with angular conservation (two left panels) and without (two rightpanels). Colours for material plots are the same as in figure 3. . . . . . . . L / L i n i t i a l correctednon − correctedfluid Figure 9: Angular momentum as function of time for the two cases shown above ( ◦ , ANEOS, solid,integrated density), with/without correction tensor and the corresponding fluid case for comparison. T [K]10 − − − − C u m u l a t i v e m a ss T [K]10 − − − − C u m u l a t i v e m a ss Figure 10: Cumulative mass fraction of mantle with temperature greater than a given value. Same as leftpanel of figure 7, but comparing density computation instead. Left panel is for ANEOS, where blue linesare now density summation. Right panel is Tillotson with green lines for density summation.
These results confirm that standard SPH does not handle well problems which involverigid body rotation (Speith, 2006) and point out the importance of applying the correctiontensor to deal with such cases, in particular for situations where the rotation timescale iscomparable to the simulated time.
As discussed in section 2, we test two different methods to compute the densities of theSPH particles. We use density integration for our nominal cases.Major differences arise at the boundaries between materials or at the surface where steepdensity gradients are located. For density summation, a group of particles separates in allquantities shown in radial profiles (figures 1 and 2). These groups contain particles closerthan h from the boundary and where the adjacent material enters in the SPH sum (orwhere there is a lack of neighbouring particles when close to the surface). Their densityis shifted towards the one beyond the boundary. For integrated density, as no such sumis involved, particles close to the boundary do not separate from the remaining. However,there are still some oscillations close to the discontinuity.Despite the issues mentioned above, which result in different initial density profiles (atthe boundaries), the results at end of the simulations do not show any major differencesbetween the two schemes (figures 5 and 6). However, there are some subtle effects. Forinstance, we observe a different behaviour of hot surface particles which were ejected duringthe early stage of impact and reaccreted later on. With ANEOS and density summation,they form a layer separated from the planet’s surface. In the other cases, this feature isless pronounced and not present at all in the cases with density integration. In the ANEOScase, the different densities of the surface particles also affect the temperature (figure 10,left panel) although the energy distribution remains essentially identical. As surface particle18ntegrated energy Constant energymaterial temperature material temperature Figure 11: End result comparison for case with evolved internal energy during setup phase (two left panels)and constant (two right panels). Colours for material plots are the same as in figure 3. density is lower in the density summation case than for density integration, ANEOS doesnot find them in the same phase: mixed liquid-vapour for the former and liquid in the latter.Density summation simulations exhibit lower temperature since part of the internal energyis used for the phase transition rather than temperature increase. This explains the absenceof material with temperature above with ANEOS density summation on the figure.The second small bump located between and in the case with solid rheologyand density summation (solid blue line on figure 10) is due to the core-mantle boundary.As phase transition does not influence the temperature in the case with Tillotson EoS, wedon’t see this feature on the right panel of figure 10.We note that the differences between the numerical schemes discussed above result fromdifferences at small scales, representing only a small fraction of the total mass involved inthe simulation. Even at the relatively high resolution (one million SPH particles) used inthe simulations performed here, some of these small-scale features are under-resolved. Thisis certainly the case for the layer of reaccreted material discussed above. Also the peakpressures and temperatures produced in the initial stages of the collision may not be fullyresolved.With the resolution used in the simulations presented here, the numerical scheme withdensity integration appears to be the better choice as it leads to a better defined surface.It is clear, however, that a more consistent treatment of the boundaries, as well as veryhigh resolution simulations, are required to obtain more accurate physical properties at suchsmall scales.
Bodies computed using the Tillotson EoS were evolved during set up phase with theinternal energy allowed to vary. This leads to a spread in temperature for those bodiesshown in figures 1 and 2 as well as the density decrease close to the centre for the impactorcomputed using density summation.For test purposes, a second series of bodies was generated with internal energy kept19oS Rheology Density Period [Sidereal day]ANEOS Solid Integrated 1.22Solid Summed 1.19Fluid Integrated 1.11Fluid Summed 1.14Tillotson Solid Integrated 1.18Solid Summed 0.95Fluid Integrated 1.13Fluid Summed 1.13
Table 3: Impact-induced rotation periods for oblique geometry runs constant. For the latter case we don’t observe a density decrease close to body centres.We don’t observe any notable difference between the results from the two different setupschemes. A comparison is shown in figure 11. The main differences are located at the core-mantle boundary where the effect of the set up method discussed here and the instabilitywith solid rheology explained in section 4.1 overlap. This effect is small compared with allthe previously discussed parameters.
The rotation periods resulting from the oblique impacts considered in this study, calcu-lated by applying a correction tensor to the angular momentum, are presented in table 3.They are comparable to Mars’ present rotation period. All fluid rheology runs have a verysimilar rotation periods ( . − .
14 Martian sideral days ). The solid rheology usually leadsto a slightly slower rotation, except for the case with Tillotson EoS and density summationwhich does not follow this trend and is rotating faster.
We provide in table 4 components of non-accreted material. The amount of unboundmaterial ("Escaping" column) remains almost constant from several hours after impactonward. At the end of the simulations, it represents the biggest component of the ejecta.We find that runs with fluid rheology lead to roughly 30% less escaping mass than the runswith solid rheology. This is visible in figure 4 as the impactor’s mantle (purple lines) isless present in the planet with solid rheology. The impactor penetrate less deeply insidethe target during the initial phase for the latter case. As consequence, there is a slightincrease in the amount of material not directly accreted by the planet, which is reflected inthe escaping mass.The amount of material which is still bound but not (yet) reaccreted is lower than theamount of escaping material by 6 hours after impact. It decreases as material fall back ontothe planet and eventually reaches a value about three times lower at the end. For a furtheranalysis of that part, the reader is referred to Golabek et al. (2017).20jecta Disc EscapingEoS Rheology Density [ M Mars ] [ M Mars ] [ M Mars ]ANEOS Solid Integrated . · − . · − . · − Solid Summed . · − . · − . · − Fluid Integrated . · − . · − Fluid Summed . · − . · − . · − Tillotson Solid Integrated . · − . · − . · − Solid Summed . · − . · − . · − Fluid Integrated . · − . · − Fluid Summed . · − . · − . · − Table 4: Non accreted masses for oblique geometry runs.
Ejecta column gives the amount of bound ejecta,computed at 6 hours after impact.
Disc is the amount of material that remains orbiting around the sim-ulations’ end (18 hours after impact) and
Escaping represents the unbound material, also at simulations’end.
Also shown in table 4 are the disc’s masses. These values fluctuate and tend to increasewith time. Particles composing the disc come from a small subset of the bound ejecta. Wefind large variations of the disc mass, depending mostly on the numerical scheme. However,the discs masses are very small (a few 10 − M Mars at most) and are under-resolved (theycontain between 138 and 332 SPH particles for density summation and up to 74 for densityintegration). Much higher resolution would be required to obtain meaningful results interms of the properties of such small discs. We note that recent studies of the impact-formation of the Martian moons (Citron et al., 2015; Canup and Salmon, 2016), whichconsidered comparable disc masses, used similar or lower resolution and did not includematerial strength.
5. Discussion and Conclusions
In this paper, we model large-scale ( ≈ ) impacts on a Mars-like planet using anupdated SPH code (Reufer et al., 2012), which includes self-gravity, a newly implementedstrength model (Jutzi, 2015) and various EoS (Tilloston and M-ANEOS). A subset of thesesimulations is used as initial conditions in a thermochemical code to study the long-terminterior evolution of Mars (Golabek et al., 2017). Here, we investigate the effects of materialstrength and of using different EoS. Some potential numerical effects (in the context ofdiscontinuities) which may result from the different ways to compute the density in SPHare studied as well. Finally, an advanced SPH scheme to avoid rotational instability in rigidbody rotation (Speith, 2006) is applied and tested in the regime of large-scale collisions.We find that in the collision regime considered here (with impact velocities of / s ), theTillotson and ANEOS equations of state lead to post-impact temperature distributions whichare quite similar in the temperature range below ∼ with only subtle differences. Thisindicates that the Tillotson EoS works reasonably well in this regime (where vaporisation21s not significant) and the simple estimation of the temperature from the specific energy byusing a constant heat capacity T = u/c p is justified.On the other hand, our results strongly suggest that the effect of material strength issubstantial even for the large-scale collisions considered here. When strength is taken intoaccount, the post-impact distributions of the impactor material and temperature are verydifferent. For instance, in this case, the heat generated by the impact is spread out overa much broader region of mantle material localised around the impact point. We expectthis finding to be generally important for collisions in this regime. We note that previousSPH simulations at this scale (e.g. Marinova et al., 2008; Canup, 2005; Citron et al., 2015;Canup and Salmon, 2016) did not include material strength. Our results also suggested thatPluto-Charon type collisions (e.g. Canup, 2005) might well be in a regime where materialstrength is not negligible.The different density computation schemes used here do not show large differences inthe final global outcome of the simulations. However, we observe some differences at smallscales, and a more consistent treatment of density discontinuities (e.g. Hosono et al., 2013)is needed. Our study confirms that the known issue of standard SPH in the case of rigidbody rotation is solved by using an appropriate correction tensor which increases the SPHconsistency (Speith, 2006).In conclusion, our results show that in the modelling of collisions at scales of ≈ ,it is essential to use a material strength model. The details of the EOS are less important inthis regime. Finally, an improved treatment of boundaries as well as a very high resolutionare required to obtain accurate physical properties at small scales. Acknowledgements
We thank Gregor J. Golabek and Taras V. Gerya for helpful discussion in the study’sdesign, and anonymous referee for suggestions that improved this paper. A.E. acknowledgesthe financial support of the Swiss National Science Foundation under grant 200020_17246.This work has been carried out within the frame of the National Centre for Competencein Research PlanetS supported by the Swiss National Science Foundation. The authorsacknowledge the financial support of the SNSF.
References
J. C. Andrews-Hanna, M. T. Zuber, and W. B. Banerdt. The Borealis basin and the origin of the martiancrustal dichotomy.
Nature , 453:1212–1215, June 2008. doi: 10.1038/nature07011.E. Asphaug. Similar-sized collisions and the diversity of planets.
Chemie der Erde / Geochemistry , 70:199–219, 2010. doi: 10.1016/j.chemer.2010.01.004.E. I. Asphaug and A. Reufer. Making an Iron Planet: The Case for Repeated Hit and Run Collisions.
AGUFall Meeting Abstracts , page A3, December 2014.J. Barnes and P. Hut. A hierarchical O(N log N) force-calculation algorithm.
Nature , 324:446–449, December1986. doi: 10.1038/324446a0.W. Benz. An introduction to computation methods in hydrodynamics. In C. B. de Loore, editor,
Late stagesof stellar evolution, computational methods in astrophysical hydrodynamics , volume 373, pages 258–312.Springer, 1991. ISBN 0387536205. . Benz and E. Asphaug. Impact simulations with fracture. I - Method and tests. Icarus , 107:98, January1994. doi: 10.1006/icar.1994.1009.W. Benz and E. Asphaug. Simulations of brittle solids using smooth particle hydrodynamics.
Comput. Phys.Commun. , 1-2:253–265, 1995. doi: 10.1016/0010-4655(94)00176-3.W. Benz, W. L. Slattery, and A. G. W. Cameron. Collisional stripping of Mercury’s mantle.
Icarus , 74:516–528, June 1988. doi: 10.1016/0019-1035(88)90118-2.W. Benz, A. Anic, J. Horner, and J. A. Whitby. The Origin of Mercury.
Space Sci. Rev. , 132:189–202,October 2007. doi: 10.1007/s11214-007-9284-1.A. G. W. Cameron and W. R. Ward. The Origin of the Moon. In
Lunar and Planetary Science Conference ,volume 7 of
Lunar and Planetary Science Conference , page 120, March 1976.R. M. Canup. Simulations of a late lunar-forming impact.
Icarus , 168:433–456, April 2004. doi: 10.1016/j.icarus.2003.09.028.R. M. Canup. A Giant Impact Origin of Pluto-Charon.
Science , 307:546–550, January 2005. doi: 10.1126/science.1106818.R. M. Canup. On a Giant Impact Origin of Charon, Nix, and Hydra.
Astron. J. , 141:35, February 2011.doi: 10.1088/0004-6256/141/2/35.R. M. Canup. Forming a Moon with an Earth-like Composition via a Giant Impact.
Science , 338:1052–,November 2012. doi: 10.1126/science.1226073.R. M. Canup and E. Asphaug. Origin of the Moon in a giant impact near the end of the Earth’s formation.
Nature , 412:708–712, August 2001.R. M. Canup and J. Salmon. On an Origin of Phobos-Deimos by Giant Impact. In
Lunar and PlanetaryScience Conference , volume 47 of
Lunar and Planetary Science Conference , page 2598, March 2016.R. I. Citron, H. Genda, and S. Ida. Formation of Phobos and Deimos via a giant impact.
Icarus , 252:334–338, May 2015. arXiv:1503.05623. doi: 10.1016/j.icarus.2015.02.011.G. S. Collins, H. J. Melosh, and B. A. Ivanov. Modeling damage and deformation in impact simulations.
Meteorit. Planet. Sci. , 39:217–231, February 2004. doi: 10.1111/j.1945-5100.2004.tb00337.x.M. Ćuk and S. T. Stewart. Making the Moon from a Fast-Spinning Earth: A Giant Impact Followed byResonant Despinning.
Science , 338:1047–, November 2012. doi: 10.1126/science.1225542.H. Frey and R. A. Schultz. Large impact basins and the mega-impact origin for the crustal dichotomy onMars.
Geophys. Res. Lett. , 15:229–232, March 1988. doi: 10.1029/GL015i003p00229.Gregor J. Golabek, Martin Jutzi, Alexandre Emsenhuber, Taras V. Gerya, and Erik I. Asphaug. CouplingSPH and thermochemical models of planets: Methodology and example of a Mars-sized body.
Icarus ,2017.W. K. Hartmann and D. R. Davis. Satellite-sized planetesimals and lunar origin.
Icarus , 24:504–514, April1975. doi: 10.1016/0019-1035(75)90070-6.N. Hosono, T. R. Saitoh, and J. Makino. Density-Independent Smoothed Particle Hydrodynamics for aNon-Ideal Equation of State.
Publ. Astron. Soc. Jpn , 65:108, October 2013. arXiv:1307.0916. doi:10.1093/pasj/65.5.108.Martin Jutzi. SPH calculations of asteroid disruptions: The role of pressure dependent failure models.
Planet. Space Sci. , 107:3–9, March 2015. arXiv:1502.01860. doi: 10.1016/j.pss.2014.09.012.Martin Jutzi, K. Holsapple, K. Wünneman, and Patrick Michel. Modeling Asteroid Collisions and ImpactProcesse. In Patrick Michel, Francesca E. DeMeo, and William F. Bottke Jr., editors,
Asteroids IV , page679. The University of Arizona Press, 2015. arXiv:1502.01844.M. M. Marinova, O. Aharonson, and E. Asphaug. Mega-impact formation of the Mars hemispheric dicho-tomy.
Nature , 453:1216–1219, June 2008. doi: 10.1038/nature07070.M. M. Marinova, O. Aharonson, and E. Asphaug. Geophysical consequences of planetary-scale impacts intoa Mars-like planet.
Icarus , 211:960–985, February 2011. doi: 10.1016/j.icarus.2010.10.032.H. J. Melosh.
Impact cratering: A geologic process . Oxford University Press, 1989.H. J. Melosh. A hydrocode equation of state for SiO2.
Meteorit. Planet. Sci. , 42:2079–2098, 2007. doi:10.1111/j.1945-5100.2007.tb01009.x.J. J. Monaghan. Smoothed particle hydrodynamics.
Annu. Rev. Astron. Astrophys. , 30:543–574, 1992. doi: Astron. Astro-phys. , 149:135–143, August 1985.F. Nimmo, S. D. Hart, D. G. Korycansky, and C. B. Agnor. Implications of an impact origin for the martianhemispheric dichotomy.
Nature , 453:1220–1223, June 2008. doi: 10.1038/nature07025.William H. Press.
Numerical recipes in C++ : the art of scientific computing . Cambridge University Press,2002. ISBN: 0-521-75033-4.C. Reinhardt and J. Stadel. Numerical aspects of giant impact simulations.
Mon. Not. R. Astron. Soc. , 467:4252–4263, June 2017. arXiv:1701.08296. doi: 10.1093/mnras/stx322.A. Reufer, M. M. M. Meier, W. Benz, and R. Wieler. A hit-and-run giant impact scenario.
Icarus , 221:296–299, September 2012. arXiv:1207.5224. doi: 10.1016/j.icarus.2012.07.021.L. E. Senft and S. T. Stewart. Dynamic fault weakening and the formation of large impact craters.
EarthPlanet. Sci. Lett. , 287:471–482, October 2009. doi: 10.1016/j.epsl.2009.08.033.Roland Speith. Improvements of the numerical method Smoothed Particle Hydrodynamics. Habilitationss-chrift, 2006.S. L. Thompson and H. S. Lauson. Improvements in the CHART-D Radiation-hydrodynamic code III:Revised analytic equations of state. Technical Report SC-RR-71 0714, Sandia National Laboratories,1972.J. H. Tillotson. Metallic Equations of State for Hypervelocity Impact. Technical Report GA-3216, GeneralAtomic, July 1962.D. E. Wilhelms and S. W. Squyres. The martian hemispheric dichotomy may be due to a giant impact.
Nature , 309:138–140, May 1984. doi: 10.1038/309138a0., 309:138–140, May 1984. doi: 10.1038/309138a0.