où on rappelle que $q$ et $Q$ sont respectivement la plus courte et la plus grande
distance de l'objet au Soleil, et $a_0$ le demi grand axe (ici le rayon) de l'orbite de la
planète. Si cette condition est respectée, cela veut dire d'après l'équation (2.13) :

$$r = \frac{a(1 - e^2)}{1 + e \cos f},$$

qu'il existe deux valeurs de l'anomalie vraie $f$ (ambiguïté du cosinus) telles que
$r = a_0$. En posant $C = \cos f$ (de la même manière, on notera $S = \sin f$), le
projectile et la cible sont à la même distance du Soleil pour

$$C = \frac{\frac{a}{a_0}(1 - e^2) - 1}{e}. \qquad (5.13)$$

Puisque les deux orbites possèdent une inclinaison relative non nulle, elles ne
présenteront une réelle intersection ($\mathbf{r} = \mathbf{a}_0$) que si l'égalité $r = a_0$ se produit au
niveau d'un noeud commun. Comme $C = \cos f = \cos(\theta - \omega)$ et qu'ici, les noeuds
communs correspondent à $\theta = 0$ ou $\theta = \pi$, une telle intersection n'est possible
que pour quatre valeurs de l'argument du périhélie $\omega$, qui sont

$$\omega_1 = \cos^{-1} C; \quad \omega_2 = -\cos^{-1} C; \quad \omega_3 = \pi - \cos^{-1} C; \quad \omega_4 = \pi + \cos^{-1} C. \quad (5.14)$$

L'hypothèse étant faite que l'argument du périhélie varie uniformément dans
le temps, ces quatre configurations sont toujours possibles tant que la condition
donnée par l'équation (5.12) est respectée, et elles sont équiprobables (on note que
dans le cas où l'orbite de la cible n'est pas circulaire, l'existence d'une intersection
est également conditionnée aux valeurs de $\Omega$ et $\Omega_0$; voir la section 5.2 qui traite
du cas général). Un impact, s'il se produit, aura lieu pour une valeur de $\omega$ proche
de $\omega_1, \omega_2, \omega_3$ ou $\omega_4$. La probabilité totale d'impact sur une période orbitale du
projectile peut s'écrire

$$P_{orb} = \sum_{n=1}^{4} P_{\omega_n} = 4P_{\omega}, \qquad (5.15)$$

où $P_{\omega_n}$ est la probabilité d'impact au voisinage de $\omega_n$, égale à $P_{\omega}$ pour les quatre
$\omega_n$.

Il s'agit maintenant d'exprimer $P_{\omega}$. On peut la décomposer ainsi

$$P_{\omega} = \frac{\Delta\omega}{2\pi} \times \frac{\Delta t}{T_0}, \qquad (5.16)$$

où $\frac{\Delta\omega}{2\pi}$ est la fraction angulaire sur l'espace des valeurs de $\omega$ qui permet à un
impact de se produire (potentiellement) autour de $\omega_n$, et où $\frac{\Delta t}{T_0}$ est la fraction de
temps pendant laquelle les deux objets sont suffisamment proche l'un de l'autre
pour que la collision ait réellement lieu, $T_0$ étant la période orbitale de la cible.
Les composantes de la vitesse relative au noeud commun, $U_∞$, se calculent dans le repère orbital (OXYZ) comme

$$ 
\left\{ 
\begin{aligned}
U_X^2 &= v_0^2 \left( 2 - \frac{a_0}{a} - \frac{a}{a_0} (1 - e^2) \right), \\
U_Y &= -v_0 \left( \sqrt{\frac{a}{a_0}} (1 - e^2) \cos i - 1 \right), \\
U_Z^2 &= v_0^2 \frac{a}{a_0} (1 - e^2) \sin^2 i,
\end{aligned}
\right. 
\qquad (5.17) 
$$

et donc

$$ U_{\infty}^{2} = v_{0}^{2} \left( 3 - \frac{a_{0}}{a} - 2\sqrt{\frac{a}{a_{0}}(1 - e^{2})} \cos i \right), \qquad (5.18) $$

avec $v_0$ la vitesse orbitale de la planète telle que $v_0 = \sqrt{\frac{(GM_S+GM_0)}{a_0}} \approx \sqrt{\frac{GM_S}{a_0}}$. On considère que $U_∞$ est la vitesse relative au vrai point de collision, celui-ci se trouvant dans le voisinage du noeud commun, voisinage quasi ponctuel au regard de la distance planète-Soleil. On rappelle que la vitesse $U_∞$ n'inclut pas la perturbation gravitationnelle de la planète. On s'autorise à la qualifier de vitesse à l'infini en considérant que la distance au-delà de laquelle la perturbation gravitationnelle liée à la cible est négligeable (rayon de la sphère de Hill) est elle aussi très petite comparée à la distance planète-Soleil. Notons que par commodité pour la suite de notre étude, nous n'utilisons pas le repère traditionnel utilisé par Opik et ses successeurs. Dans ce dernier (appelons-le (Oxyz)), $U_x = -U_X$ et $U_y = -U_Y$.

Il y a collision lorsque la distance qui sépare les deux corps est inférieure à une certaine valeur $d_{min}$. Si on ne prend pas en compte l'attraction gravitationnelle, et si on considère que le projectile est suffisamment petit au regard de la cible pour être assimilé à un objet ponctuel, on a $d_{min} = R_0$, où $R_0$ est le rayon de la planète cible. En réalité, la planète, massique, attire à elle par gravité un certain nombre de projectiles qui passent à proximité. L'impact aura lieu si $d_{min} \le \tau$, où $\tau$ est le rayon de la section gravitationnelle de la cible, donné par l'équation (5.10).

La fraction angulaire sur l'espace $\omega$ qui permet à la distance minimale entre les deux orbites d'être inférieure à $\tau$ s'exprime

$$ \frac{\Delta\omega}{2\pi} = \frac{\tau}{a_0} \frac{1}{\pi \sin i} \frac{\sqrt{U_{\infty}^2 - U_Y^2}}{|U_X|} . \qquad (5.19) $$

La fraction temporelle est quant à elle donnée par

$$ \frac{\Delta t}{T_0} = \frac{\tau}{4a_0} \frac{U_\infty}{\sqrt{U_\infty^2 - U_Y^2}}. \qquad (5.20) $$

La probabilité totale d'impact par unité de temps est donnée par $P = P_{orb}/T$, où $T = 2\pi a^{3/2}/\sqrt{GM_S}$ est la période orbitale de l'astéroïde. La probabilité d'impact annuelle, puisque 1 an $\approx \frac{\sqrt{GM_S}}{2\pi}$, est donc
$$P = \frac{\tau^2 U_{\infty}}{a_0^2 |U_X|} \frac{1}{\pi \sin i} a^{-3/2}. \qquad (5.21)$$

On le voit, la probabilité d'impact est proportionnelle à $\tau^2$: plus la planète est grande et massique, et plus elle est susceptible de se faire bombarder (pour une géométrie orbitale identique).

On note que sur un cycle de précession de l'orbite du projectile, il y a quatre géométries possibles de collision. En effet, celle-ci peut se faire au niveau du noeud ascendant ou descendant, et quand le projectile s'éloigne du Soleil ou s'en approche. Selon que le noeud est ascendant ou descendant, l'inclinaison mutuelle entre les deux objets (ici, simplement l'inclinaison du projectile $i$), est par convention positive ou négative. $U_Z$ est donc donné par l'équation (5.17) comme $U_Z = \pm\sqrt{U_Z^2}$. Il en va de même pour $U_Z$ selon que l'approche se fait dans le sens planète-Soleil ou inversement ($U_X = \pm\sqrt{U_X^2}$). Chacune des quatre géométries résultantes est équiprobable. Ainsi, par exemple, une collision de probabilité $P$ aura lieu dans 25 % des cas alors que le projectile plonge vers la planète ($U_Z > 0$) avant d'avoir atteint son périhélie ($U_X > 0$).

Il faut enfin mettre en évidence le fait que la probabilité d'impact donnée par l'équation (5.21) présente deux singularités. Elle tend vers l'infini quand $i$ ou bien $U_x$ tendent vers zéro. Le premier cas est soluble en posant $\sin i = \tau/a_0$ quand $\sin i < \tau/a_0$. Le deuxième, qui correspond à une collision au niveau du périhélie ou de l'aphélie du projectile, se résout en posant $U_x = U_{min}$ quand $U_x < U_{min}$, avec $U_{min} = \frac{\tau}{a_0} v_0 (\frac{a}{a_0}(1 - e^2))$. Cela permet de conserver la validité du calcul avec une précision du troisième ordre [Dones et al., 1999].

### 5.2.1.2 Orbites hyperboliques

Lorsque les astéroïdes ou comètes possèdent une orbite non bornée, l'expression de la probabilité d'Opik reste similaire, à quelques différences près. Typiquement, ce cas de figure survient lorsque l'on s'intéresse au bombardement d'un satellite gravitant autour de sa planète mère. Le rôle du corps primaire (le Soleil dans la section précédente) est alors joué par la planète. Le corps cible (précédemment la planète) devient le satellite. Le projectile, capturé par la sphère d'influence de la planète, décrit une orbite hyperbolique autour de celle-ci, et non plus une orbite elliptique.

Sous les mêmes approximations que précédemment, la probabilité d'impact s'exprime exactement comme donnée par l'équation (5.21), à ceci près que les
grandeurs affublées de l'indice *l* font référence au satellite, et que l'indice *p* se substitue à l'indice *s*, indiquant que la planète prend ici la place du Soleil. Puisque l'orbite du projectile n'est pas bornée, il convient d'exprimer la probabilité d'impact par orbite, et non plus par unité de temps. On a donc

$$P = \frac{\tau^2 U'_{\infty}}{a_l^2 |U'_{x}|} \frac{1}{\pi \sin i}, \qquad (5.22)$$

où *P* est définie pour

$$q \le a_l, \qquad (5.23)$$

avec $q = a(e - 1)$ le péricentre de l'hyperbole. On note que $a,e$ et $i$ font ici référence aux éléments orbitaux de l'hyperbole, $i$ étant l'inclinaison du projectile par rapport au plan orbital du satellite.

La différence principale réside dans l'expression de la vitesse relative d'approche, qui fait intervenir le péricentre $q$. Comme indiqué dans Shoemaker et Wolfe [1982], celle-ci est donnée par

$$\left\{ \begin{aligned} U'_{\infty} &= v_l^2 \left( 3 - \frac{a_l(1-e)}{q} - 2\sqrt{\frac{q}{a_l}(1+e)} \cos i \right), \\ U'_{X} &= v_l^2 \left( 2 - \frac{a_l(1-e)}{q} - \frac{q}{a_l}(1+e) \right), \\ U'_{Y} &= -v_l \left( \sqrt{\frac{q}{a_l}(1+e)} \cos i - 1 \right), \\ U'_{Z} &= v_l^2 \frac{q}{a_l}(1+e) \sin^2 i, \end{aligned} \right. \qquad (5.24)$$

avec $v_l = \sqrt{\frac{GM_p}{a_l}}$.

Rappelons que la probabilité d'impact donnée par l'équation (5.22) n'est valide que sous l'hypothèse que l'argument du péricentre, $\omega$, varie uniformément. Les travaux de Zahnle et al. [1998] et Zahnle et al. [2001] utilisent cette formule pour décrire le bombardement des satellites joviens. Pourtant, si l'hypothèse apparaît valide dans le cas d'orbites elliptiques qui précessent autour du soleil, voire de comètes hyperboliques du système solaire externe sous la très forte domination de Jupiter, il n'en est rien dans le cas qui nous occupe principalement dans cette étude, à savoir le bombardement de la Lune. En effet, les orbites elliptiques qui pénètrent dans le système Terre-Lune pour se changer en orbites hyperboliques ne sont plus alors animées d'un mouvement de précession. Faire l'hypothèse de l'uniformité de $\omega$ reviendrait à perdre l'information sur les conditions d'approche avec le système Terre-Lune ($i_{\infty}$), calculées plus loin. Or, nous verrons que ces conditions d'approche sont déterminantes, en particulier pour mettre en évidence un effet de variation en latitude du flux d'impacts lunaire.
### 5.2.2 Cas général

Nous nous plaçons ici dans le cas où la planète cible possède une excentricité $e_0$ et une inclinaison $i_0$ non nulles. Cela complique sensiblement l'expression et le calcul de la probabilité d'impact, qui dépend alors de deux variables supplémentaires, $\omega_0$ et $\Omega_0$. Nous attestons ici des principales formules et de la démarche générale qui mènent à l'obtention de cette probabilité, le détail de celles-ci pouvant être trouvé dans le travail original de Greenberg [1982] et Bottke et Greenberg [1993]. Notons que Öpik [1951] incluait une simple facteur correctif à ses équations pour tenir compte de l'excentricité de la cible (voir aussi Shoemaker et Wolfe [1982]), ce qui fut démontré comme étant une source importante d'erreurs [Wetherill, 1967; Greenberg, 1982].

Avant toute chose, une collision n'est possible que si

$$q < Q_0 \text{ et } Q > q_0. \qquad (5.25)$$

Si ces conditions sont respectées, il faut encore que l'orientation des deux orbites soit favorable.

On rappelle que la distance qui sépare le projectile du Soleil est donnée par

$$r = \frac{a(1 - e^2)}{1 + eC},$$

où C est le cosinus de l'anomalie vraie. Le rayon orbital de la cible, $r_0$, est donné de la même manière, en substituant $a, e, i$ et $C$ par $a_0, e_0, i_0$ et $C_0$. On rappelle également que $C, C_0, S$ et $S_0$ sont une notation abrégée pour $\cos(\theta - \omega)$, $\cos(\theta_0 - \omega_0)$, $\sin(\theta - \omega)$ et $\sin(\theta_0 - \omega_0)$, respectivement.

Puisque chacun des deux objets possède une inclinaison non nulle, l'angle qui sépare les deux plans orbitaux varie selon la valeur des longitudes du noeud ascendant, $\Omega$ et $\Omega_0$, comme indiqué sur la figure 5.4. Cet angle, que nous appellerons l'inclinaison mutuelle $i_m$, est donné par

$$\cos i_m = \cos i \cos i_0 + \sin i \sin i_0 \cos \Delta\Omega, \qquad (5.26)$$

où $\Delta\Omega = \Omega_0 - \Omega$. On rappelle que $\Omega$ et $\Omega_0$ varient uniformément. Le signe de $i_m$ est choisi positif si le noeud de rencontre est le noeud ascendant du projectile, c'est-à-dire si le projectile coupe le plan orbital de la planète de bas en haut. Dans le cas contraire, $i_m$ est négatif. La symétrie du problème rend l'un ou l'autre cas équiprobables.

En se référant à nouveau à la figure 5.4, on en déduit que la valeur de $\theta_0$ au noeud commun est

$$(\theta_0)_{\text{n.c.}} = \Theta_0 = \tan^{-1} \left( \frac{-\sin \Delta\Omega}{(\cot i \sin i_0 - \cos \Delta\Omega \cos i_0)} \right). \qquad (5.27)$$
<figure><img src="image_29.png" /><figcaption>FIG. 5.4 – Inclinaison mutuelle et noeud commun.</figcaption></figure>

Une autre valeur de $\Theta_0$, distante de $\pi$, est possible. La symétrie du problème permet pourtant de se cantonner à la valeur retournée par l'équation (5.27).

Une condition nécessaire à ce que la distance entre les deux orbites soit nulle ($\mathbf{r} = \mathbf{r}_0$) est $r = r_0$, c'est-à-dire

$$C = \frac{1}{e} \left[ (1 + e_0 C_0) \frac{a(1 - e^2)}{a_0(1 - e_0^2)} - 1 \right]. \qquad (5.28)$$

Là encore, quatre géométries d'impact équiprobables sont possibles. Ce sont les mêmes que pour le cas circulaire traité plus haut. Ici, elles se traduisent par le signe des sinus de l'anomalie vraie du projectile et de la cible, $S$ et $S_0$, indifféremment (et indépendamment) positifs ou négatifs.

Plaçons nous dans le repère (OXYZ) au moment où $\mathbf{r} = \mathbf{r}_0$, comme illustré par la figure B.1. Notons que pour rester aussi fidèle que possible à l'article de Greenberg [1982], les axes ont été inversés. Les vecteurs vitesse de l'astéroïde et de la cible font un angle $\alpha$ et $\alpha_0$ avec l'axe ($-OX$). L'angle $\alpha$ est donné par

$$\cos \alpha = \frac{eS}{\sqrt{1 + 2eC + e^2}} \quad \text{et} \quad \sin \alpha = \frac{1 + eC}{\sqrt{1 + 2eC + e^2}}. \qquad (5.29)$$

De la même façon, $\alpha_0$ s'exprime en substituant $e$, $C$ et $S$ par $e_0$, $C_0$ et $S_0$ dans l'équation précédente.
<figure><img src="image_30.png" /><figcaption>FIG. 5.5 – Géométrie du mouvement autour du noeud commun.</figcaption></figure>

La vitesse de la particule au noeud commun est

$$
\begin{align}
v^2 &= \frac{GM_s}{a(1-e^2)}(1 + 2eC + e^2), \nonumber \\
v_X &= -v \cos \alpha, \tag{5.30} \\
v_Y &= -v \sin \alpha \cos i_m, \nonumber \\
v_Z &= v \sin \alpha \sin i_m, \nonumber
\end{align}
$$

tandis que la vitesse de la cible est

$$
\begin{align}
v_0^2 &= \frac{GM_s}{a_0(1-e_0^2)}(1 + 2e_0C_0 + e_0^2), \nonumber \\
v_{0,X} &= -v_0 \cos \alpha_0, \tag{5.31} \\
v_{0,Y} &= -v_0 \sin \alpha_0, \nonumber \\
v_{0,Z} &= 0. \nonumber
\end{align}
$$

La vitesse relative d'approche est donnée par

$$
\mathbf{U}_{\infty} = \mathbf{v} - \mathbf{v}_{0}. \tag{5.32}
$$

L'expression de la probabilité de collision se base sur le cas **r** = **r**₀, et inclus
la distance dans l'espace (C, C₀) pour laquelle ||**r** − **r**₀|| < τ, ainsi que la fraction
de temps qui permet à cette condition d’être respectée. En reprenant les mêmes notations que dans le cas circulaire,

$$ \frac{\Delta\omega}{2\pi} = \frac{\tau}{\pi} \frac{A(1 + eC)^2}{a(1 - e^2)e|S| \sin i_m}, \qquad (5.33) $$

avec

$$ A = [\cot^2 \alpha + \cot^2 \alpha_0 + \sin^2 i_m - 2 \cot \alpha_0 \cot \alpha \cos i_m]^{1/2}, \qquad (5.34) $$

et

$$ \frac{\Delta t}{T_0} = \frac{\pi \tau}{2 T_0} \frac{U_\infty / v_0}{[U_\infty^2 - (U_{\infty,x} \cos \alpha_0 + U_{\infty,y} \sin \alpha_0)^2]^{1/2}}. \qquad (5.35) $$

Pour un couple $(\omega_0, \Delta\Omega)$ donné, et pour une des quatre géométries possibles, la probabilité de collision s'exprime

$$ P_{\omega_0, \Delta\Omega} = \frac{1}{(2\pi)^2} \times \frac{\Delta\omega}{2\pi} \times \frac{\Delta t}{T_0} \times \frac{1}{T}, \qquad (5.36) $$

où le facteur $\frac{1}{(2\pi)^2}$ tient compte de l'espace total sur $\omega_0$ et $\Delta\Omega$.

Afin d'obtenir la probabilité de collision totale selon toute les orientations orbitales possibles, $P_{\omega_0,\Delta\Omega}$ doit être évaluée aux quatre points correspondant à un couple $(C, C_0)$ donné (c'est-à-dire à chaque quadruplet $C, C_0, \pm S, \pm S_0$) comme il a été mis en évidence par Bottke et Greenberg [1993], puis intégrée (numériquement) sur l'espace $\omega_0$ and $\Delta\Omega$. Pourtant, comme il est noté dans Bottke et al. [1994], la vitesse et l'inclinaison d'approche dépendent de $\omega_0$ et $\Delta\Omega$. Puisque nous souhaitons obtenir la distribution réelle des conditions d'approche et non la distribution de leurs valeurs moyennes sur un cycle de precession, nous évaluerons la probabilité de collision et les conditions associées pour chaque orientation orbitale donnée par le couple $(\omega_0, \Delta\Omega)$. Notons que si l'orbite de la planète cible est parfaitement circulaire, cela n'est pas nécessaire.

## 5.3 Formules applicables au satellite

Dans le cas d'une lune, une autre étape intervient. Les projectiles pénètrent dans la sphère d'influence planétaire; leur trajectoire elliptique (par rapport au Soleil) est alors hyperbolique par rapport à la planète, et ses paramètres orbitaux sont modifiés par la masse de la planète. On note $U'$ les conditions d'approche avec le système planète-lune. Il s'agit de déterminer les conditions d'approche avec la lune, $U_\infty$, ainsi que leur probabilité, en fonction de $U'$. On se place dans le cas d'une orbite lunaire circulaire d'inclinaison nulle. Cette hypothèse simplifie
radicalement les calculs présentés plus loin, et s'avère excellente pour la grande majorité des satellites.

<figure><img src="image_31.png" /><figcaption>FIG. 5.6 – Bombardement d'un système planète-lune généré lors d'une phase d'approche décrite par $\mathbf{U}'$. Les projectiles sont uniformément distribués sur le disque ($\sigma$)</figcaption></figure>

La géométrie du problème est illustrée en figure 5.6. Là encore, chaque approche $\mathbf{U}'$ génère une nuée de projectiles qui pénètrent dans le système avec un vecteur vitesse parallèle et de même norme. La révolution du satellite engendre une symétrie qui permet de s'affranchir de l'angle $j'$. Posons $j' = 0$. La position initiale sur le disque est donnée par $b$ et $\delta$, représentés en figure 5.3. L'expression des éléments orbitaux d'une trajectoire hyperbolique ($q, e, i, \dots$) dans notre géométrie est détaillée en annexe A. Les grandeurs qui correspondent au satellite ont pour indice la lettre $l$. L'indice $p$ fait référence à la planète.
L'impact est possible pour $q \le a_l$, ce qui correspond à une valeur de $b_{max}$ donnée en annexe par l'équation (A.22) telle que

$$b_{max} = a_l \sqrt{1 + 2 \left( \frac{v_l}{U'} \right)^2}$$

<figure><img src="image_32.png" /><figcaption>FIG. 5.7 – Noeuds communs de l'hyperbole qui pénètre dans le système planète-lune. Pour une valeur de $\delta$, il peut exister deux valeurs de $b$ qui autorisent l'existence d'un noeud commun, ascendant ($b_a$), ou descendant ($b_d$).</figcaption></figure>

Au sein de la nuée de projectiles distribués sur le disque, il existe potentiellement, pour une valeur de $\delta$, deux valeurs de $b$ qui permettent à certains projectiles d'intersecter exactement l'orbite lunaire, correspondant respectivement à un noeud ascendant ou descendant, et notées respectivement $b_a$ et $b_d$ (voir la figure 5.7). La symétrie du problème, engendrée par le fait que $i' \in [-\pi/2, \pi/2]$ est positif ou négatif avec la même probabilité, permet de se cantonner au noeud descendant. Cherchons donc $b_l$ tel que projectile et satellites se croisent en un même point, le noeud descendant. Notons que l'existence de $b_l$ n'est pas systématiquement permise.

La distance qui sépare le projectile sur sa trajectoire hyperbolique du centre du système (centre de la planète) s'écrit, en posant $C = \cos f$ (Eq (2.17))

$$r = \frac{q(1 + e)}{1 + eC} \qquad (5.37)$$
La valeur de C telle que $r = a_l$ est donc

$$C_l = \frac{q/a_l(1+e) - 1}{e} . \qquad (5.38)$$

Les orbites ne peuvent s'intersecter qu'au noeud commun. Puisque l'inclinaison du satellite est nulle dans le plan de référence, cela correspond au moment où le projectile coupe ce plan. Pour le noeud descendant, on a $\theta = \pi$ et donc $C = \cos(\pi - \omega) = -\cos \omega$. En identifiant $C_l$ avec $-\cos \omega$, et en exprimant $e, q$ et $\cos \omega$ en fonction de $b$ à l'aide des équations (A.19), (A.21) et (A.32), on trouve

$$\left(\frac{U'}{v_l}\right)^2 \left(\frac{b_l}{a_l}\right)^2 - \left(\frac{U'}{v_l}\right)^2 \frac{b_l \sin i'}{a_l \sin i} - \sin \delta \frac{\cos i'}{\sin i} - 1 = 0, \qquad (5.39)$$

où $i$ est l'inclinaison de l'orbite du projectile par rapport au plan orbital lunaire, donnée par les équations (A.28) et (A.29) telle que

$$\cos i = \cos \delta \cos i' \quad \text{et} \quad \sin i = \sqrt{1 - \cos^2 \delta \cos^2 i'}. \qquad (5.40)$$

La solution (physiquement acceptable) de l'équation (5.39) est

$$b_l = \frac{a_l}{2} \left( \frac{\sin i'}{\sin i} + \sqrt{\left(\frac{\sin i'}{\sin i}\right)^2 + 4 \left(1 + \frac{\tan \delta}{\tan i}\right) \left(\frac{v_l}{U'}\right)^2} \right) . \qquad (5.41)$$

Pour une approche $U'$ et une valeur donnée de $\delta$, on connaît donc la valeur du paramètre d'impact qui permet l'intersection exacte des orbites du projectile et de la lune, sa cible (quand le projectile est en phase descendante). A partir de la valeur de $b_l$, il est possible d'établir la probabilité d'impact $P_l(U', i', \delta)$ avec la lune, ainsi que les conditions associées d'approche ($U_\infty, i_\infty, j_\infty$). Pour ce faire, nous nous inspirons des formules établies par J. Öpik.

Connaissant $b_l$, la vitesse d'approche $U_\infty$ et ses trois composantes sont données par l'équation (5.24) correspondant au cas hyperbolique. En exprimant $e$ et $q$ en fonction de $b_l$ à l'aide des équations (A.19) et (A.21), on a

$$\left\{ \begin{aligned} U_{\infty}^{2} &= U'^{2} + 3v_{l}^{2} - 2U'v_{l}\frac{b_{l}}{a_{l}}\cos i, \\ U_{X}^{2} &= U'^{2}\left(1 - \left(\frac{b_{l}}{a_{l}}\right)^{2}\right) + 2v_{l}^{2}, \\ U_{Y} &= v_{l} - U'\frac{b_{l}}{a_{l}}\cos i, \\ U_{Z}^{2} &= U'^{2}\left(\frac{b_{l}}{a_{l}}\right)^{2}\sin^{2}i, \end{aligned} \right. \qquad (5.42)$$

Les angles d'approche $i_\infty$ et $j_\infty$ en découlent via l'équation (5.7).

La différence essentielle avec les probabilités d'Öpik est que ces dernières s'appliquent en considérant que l'argument du péricentre $\omega$ varie uniformément. Cela ne saurait être vrai ici, comme il est expliqué au début de cette section.
On rappelle que les probabilités d'Opik s'obtiennent en faisant intervenir le produit

$$ \frac{\Delta\omega}{2\pi} \times \frac{\Delta t}{T_0}, \qquad (5.43) $$

où $\frac{\Delta\omega}{2\pi}$ est la fraction angulaire sur l'espace des valeurs de $\omega$ qui permet à un impact de se produire autour de $\omega_n$, $\omega_n$ étant l'argument du périhélie pour lequel les deux orbites s'intersectent exactement. Le second terme du produit, $\frac{\Delta t}{T_0}$, représente la fraction de temps pendant laquelle les deux objets sont suffisamment proche l'un de l'autre pour que la collision ait réellement lieu, $T_0$ étant la période orbitale de la cible. De manière analogue, nous définissons notre probabilité comme

$$ P_l(U', i', \delta) = \frac{\Delta(b^2)}{b_{max}^2} \times \frac{\Delta t}{T_l}, \qquad (5.44) $$

où $\Delta(b^2)$ est l'intervalle autour de $b_l^2$ tel que la distance minimale qui sépare les deux orbites est inférieure à $\tau_l$, le rayon de la section gravitationnelle de la lune. Ce sont les variations du carré de $b$ qui nous intéressent ici, puisqu'il peut prendre une valeur équiprobable entre 0 et $b_{max}^2$, alors que $b$ n'est pas uniforme. Le terme $\frac{\Delta t}{T_l}$, où $T_l$ est la période orbitale de la lune, s'exprime exactement comme dans le cas des probabilités d'Opik, c'est-à-dire, avec les notations de cette section,

$$ \frac{\Delta t}{T_l} = \frac{\tau_l}{4a_l} \frac{U_\infty}{\sqrt{U_\infty^2 - U_Y^2}}. \qquad (5.45) $$

Reste à exprimer $\Delta(b^2)$. Au noeud, on rappelle que $\theta = \pi$. En notant $\ell_{max}$ l'écart maximal de $r_{\theta=\pi}$ à $a_l$ tel que la distance minimale entre les deux orbites est égale à $\tau_l$, $\Delta(b^2)$ vérifie la relation

$$ \frac{2\ell_{max}}{\Delta(b^2)} \approx \left( \frac{\partial r_{\theta=\pi}}{\partial b^2} \right)_{b=b_l}, \qquad (5.46) $$

tel qu'il est démontré pour un cas similaire dans l'annexe B.

Or, d'après l'équation (B.4),

$$ \ell_{max} = \tau_l \sqrt{1 + U_X^2 / U_Z^2}. \qquad (5.47) $$

L'équation (2.17) donne $r_{\theta=\pi}$ en fonction de $q$, $e$ et $\cos\omega$. On exprime $q$, $e$ et $\cos\omega$ en fonction de la géométrie du problème avec les équations (A.19), (A.21) et (A.32). La distance $r$ n'est plus alors fonction que de $b$, $\delta$, $U'$, $i'$ et de l'angle $\theta$. En dérivant par rapport à $b^2$, puis en posant $b = b_l$, on obtient

$$ \left( \frac{\partial r_{\theta=\pi}}{\partial b^2} \right)_{b=b_l} = \frac{1}{2a_l} \left( \frac{U'}{v_l} \right)^2 \sin i \frac{2 \sin i + 2 \cos i' \sin \delta + \left( \frac{U'}{v_l} \right)^2 \frac{b_l}{a_l} \sin i'}{\left( \sin i + \cos i' \sin \delta + \left( \frac{U'}{v_l} \right)^2 \frac{b_l}{a_l} \sin i' \right)^2}. \qquad (5.48) $$
La probabilité s'exprime enfin

$$P_l(U', i', \delta) = \frac{\tau_l^2 U_\infty v_l^2}{a_l^2 |U_Z| 2v_l^2 + U'^2} \frac{1}{\sin i} \frac{\left(\sin i + \cos i' \sin \delta + \frac{U'|U_Z|}{v_l^2} \frac{\sin i'}{\sin i}\right)^2}{2 \sin i + 2 \cos i' \sin \delta + \frac{U'|U_Z|}{v_l^2} \frac{\sin i'}{\sin i}} \quad (5.49)$$

Pour respecter la symétrie du problème, on calcule cette probabilité pour $+i'$ et $-i'$. On rappelle que $b_l$ n'est pas toujours défini. Dans ce cas, on pose $P_l = 0$. Le signe de la composante verticale de la vitesse d'approche avec la lune est choisi de manière équitable : $U_Z = \pm\sqrt{U_Z^2}$.

Le signe de $U_X$ se détermine comme suit. Comme le montre la figure B.1, $U_X < 0$ pour $\cos \alpha > 0$. D'après l'équation (B.1), cela correspond à $S = \sin f = \sin(\theta - \omega) > 0$, soit $\omega \in [0, \pi]$ pour le noeud descendant ($\theta = \pi$). $\omega$ est donné par l'équation (A.34), et l'on a

$$U_X = -\sqrt{U_X^2} \quad \text{pour} \quad \omega \in [0, \pi] \qquad (5.50)$$

$$U_X = \sqrt{U_X^2} \quad \text{pour} \quad \omega \in ]\pi, 2\pi[ \qquad (5.51)$$

$$\qquad (5.52)$$

Si $U_X$ est négatif (le projectile s'éloigne de la planète) alors que $q \le R_p$, une collision avec la planète a nécessairement déjà eu lieu, et l'on pose $P_l(U', i', \delta) = 0$. La valeur de $b$ correspondant à $q = R_p$ est

$$b_p = \sqrt{2R_p a_l \left(\frac{v_l}{U'}\right)^2 + R_p^2} \qquad (5.53)$$

La probabilité d'impact avec la planète est donc

$$P_p(U', i', \delta) = \frac{1}{2} \frac{b_p^2}{b_{max}^2} - P_l(U', i', \delta) = \frac{1}{2} \frac{R_p^2 + 2R_p a_l \left(\frac{v_l}{U'}\right)^2}{a_l^2 + 2a_l^2 \left(\frac{v_l}{U'}\right)^2} - P_l(U', i', \delta), \quad (5.54)$$

où le facteur 1/2 vient du fait que nous nous sommes restreints à l'un des deux noeuds possibles.

# 5.4 Méthodologie

## 5.4.1 Planètes

Nous résumons ici étape par étape la démarche qui permet la construction de la distribution de probabilité des conditions d'approche pour chaque planète ou
satellite, $p(\mathbf{U}_\infty)$ (soit $p(U_\infty, i_\infty, j_\infty)$), ainsi que le calcul du flux d'impacts relatif entre celles-ci.

Pour les planètes, la démarche est relativement simple. Dans le cas général (orbite planétaire excentrique), on procède comme suit. Les projectiles sont choisis selon une triade (a, e, i) aléatoire. Si le périhélie et l'aphélie de l'orbite correspondante vérifient l'équation (5.25), la collision est possible, et le processus continue en choisissant une valeur aléatoire de l'orientation relative $\Delta\Omega$ entre 0 et $2\pi$, à partir de laquelle l'inclinaison mutuelle $i_m$ est déterminée par l'équation (5.26). L'angle $\Delta\Omega$ donne également la position du noeud commun sur l'orbite de la cible, $\Theta_0$, via l'équation (5.27). L'argument du périhélie, $\omega_0$, est choisi aléatoirement dans le domaine de valeurs permises telles que $C_0 = \cos(\Theta_0 - \omega_0)$ retourne une valeur définie de C (entre -1 et 1) via l'équation (5.28). A partir de C, on trouve $S = \pm\sqrt{1-C^2}$, dont le signe est choisi de manière équitable. $S_0$ est calculé de la même manière. Les couples $(C_0, S_0)$ et $(C, S)$ déterminent une des quatre géométries de collision possibles, et permettent de calculer $U_\infty$, $i_\infty$, $j_\infty$ et $P_{\omega_0,\Delta\Omega}$ (en utilisant successivement les équations (5.29), (5.31), (5.30), (5.32), (5.7), (5.33), (5.35) et (5.36)). Puisque la valeur de $P_{\omega_0,\Delta\Omega}$ est singulière quand $i_m = 0$, on pose $\sin i_m = \tau/r_0$ quand $\sin i_m < \tau/r_0$. Qui plus est, quand la collision a lieu près de l'aphélie ou du périhélie du projectile, S tend vers zéro et la probabilité vers l'infini. D'après Greenberg [1982], cela est évité en permutant de manière appropriée la définition du projectile et de la cible (voir la référence pour plus de détails). La probabilité de collision est finalement pondérée par le nombre relatif d'objets dont l'orbite est décrite par la triade (a, e, i), telle qu'elle a été déterminée au chapitre 3 par la quantité $n(a, e, i)$. La procédure est répétée de sorte que $p(\mathbf{U}_\infty)$ converge, et le nombre d'itérations est le même pour chaque planète. En sommant toutes les probabilités de collision calculées en construisant $p(\mathbf{U}_\infty)$, puis en divisant par l'aire de la surface planétaire ($4\pi R_0^2$) nous obtenons le flux d'impact relatif entre les planètes. On note que si les variations séculaires de l'orbite de la cible sont incluses dans la simulation, une autre étape qui consiste à échantillonner $e_0$ et $i_0$ selon leur distribution de probabilité (cf. chapitre 4) est nécessaire à chaque itération. Si on considère que l'orbite de la planète est circulaire (approximation excellente pour la Terre et Venus), le calcul est très simplifié : on évalue les conditions d'approche ainsi que la probabilité de collision de chaque triade (a, e, i) avec les équations (5.17), (5.7) et (5.21).
## 5.4.2 Satellites

On commence par calculer la distribution de probabilité des conditions d'approche avec le système planète-lune global, que l'on note $p'(U', i')$. Notons que cette distribution de probabilité est intégrée sur toutes les valeurs possibles de l'angle $j'$, ce dernier ne jouant aucun rôle en raison de la symétrie du problème. $p'$ est calculée exactement de la même façon que dans le cas des planètes décrit juste au-dessus, à ceci près que la section gravitationnelle est ici posée comme étant

$$ \tau' = a_l \sqrt{1 + 2 \left( \frac{v_l}{U'} \right)^2}, \qquad (5.55) $$

et se substitue à $\tau$ dans l'expression de $P_{\omega_0, \Delta\Omega}$ (équation (5.36)).

Pour chaque couple $(U', i')$, on tire $N$ fois une valeur de $\delta$ entre 0 et $2\pi$, $N$ étant aussi grand que possible. On rappelle que $i'$ est positif ou négatif de manière équiprobable, et l'on choisit aléatoirement son signe à chaque itération. On détermine l'inclinaison $i$ de l'hyperbole qui traverse le système via l'équation (5.40), puis la valeur du paramètre d'impact $b_l$ selon l'équation (5.41). Si $b_l$ n'est pas défini, on pose la probabilité d'impact avec la lune comme nulle. A partir de $b_l$, on calcule les composantes de la vitesse d'approche par l'équation (5.42), où $U_Z = \pm\sqrt{U_Z^2}$ de manière équiprobable et où le signe de $U_X$ est donné par l'équation (5.52) selon la valeur de $\omega$, elle-même donnée par l'équation (A.34). Les angles d'approche sont donnés par l'équation (5.7). On en déduit la probabilité de collision $P_l$ donnée par l'équation (5.49) si $U_X > 0$ et $b_l > b_p$, avec $b_p$ donné par l'équation (5.53). Dans le cas contraire, on pose $P_l = 0$. On calcule ensuite la probabilité d'impact avec la planète $P_p$ (équation (5.54)). On pondère enfin la contribution de chaque projectile en multipliant $P_l$ et $P_p$ par $p'(U', i')$. En répétant cette procédure sur l'ensemble des conditions d'approche avec le système planète-satellite, on construit la distribution des conditions d'approche avec le satellite $p(U_\infty, i_\infty)$. Le bombardement relatif entre la planète et sa lune s'obtient en prenant le ratio de la somme des $P_p$ sur la somme des $P_l$.

## 5.4.3 Traitement simplifié du cas du satellite

Pour traiter le cas du satellite (on parlera du système Terre-Lune ici), une approche plus simpliste est possible, qui consiste à considérer que la Lune se situe en lieu et place de la Terre, et que cette dernière n'existe pas. La pseudo-orbite lunaire est donc l'orbite terrestre, tandis que la Lune conserve sa masse et son rayon. On applique alors certains termes correctifs pour reproduire les
conditions réelles. On suppose donc que le bombardement relatif entre la Terre et la Lune ne dépend que de leur rayon et masse respectifs, et n'est pas sensible à deux effets opposés :

- un effet de bouclier terrestre qui réduit le nombre des impacts lunaires.

- un effet de lentille gravitationnelle qui concentre une partie des projectiles vers la Lune.

Nous verrons que, pour notre modèle de projectiles et pour la distance Terre-Lune actuelle, cette approche s'avère être une approximation raisonnable. La probabilité d'impact est calculée à partir de l'équation (5.21) (en supposant ici par souci de simplicité que l'orbite terrestre est parfaitement circulaire), avec $R_0 = R_l, M_0 = M_l, a_0 = a_p$. La vitesse d'approche calculée dans ce cas reste exactement la même que si elle avait été calculée pour la Terre; notons-la par conséquent $U_\infty^p$ (par contre, la probabilité d'impact diffère, car elle dépend de $\tau$, donc de $R_0$ et $M_0$; la Terre engendre de plus forte probabilité d'impact pour les faibles vitesses que la Lune, car le pouvoir d'attraction terrestre s'exerce alors à plein). Il est possible d'inclure certains termes correctifs. La vitesse d'approche lunaire peut être recalculée pour tenir compte de la distance Terre-Lune : par conservation de l'énergie mécanique, il peut aisément être montré que la vitesse recalculée à une distance correspondant à celle de l'orbite réelle du satellite ($a_l$) est $v_\infty^l = \sqrt{U_\infty^p 2 + 2GM_p/a_l}$. La probabilité d'impact s'en trouve alors modifiée, en substituant $v_\infty^l$ à $U_\infty^p$ dans l'expression de $\tau$ (et uniquement à cet endroit). Cet accroissement de la vitesse reste toutefois très modeste : pour $U_\infty^p=10$ et $20 \text{ km.s}^{-1}$, on trouve $v_\infty^l=10.1$ and $20.05 \text{ km.s}^{-1}$, respectivement. Shoemaker et Wolfe [1982] appliquent cette correction dans le cas des satellites galiléens, et ajoutent un terme qui corrige de l'effet de concentration gravitationnelle que Jupiter exerce dans le voisinage de ses satellites, en multipliant l'expression (5.21) par un facteur $1 + \frac{2GM_p/a_l}{U_\infty^p}$. Ce facteur corrige la valeur du flux relatif entre la Terre et la Lune, mais pas la distribution des angles d'approche. D'ailleurs, si $v_\infty^l$ correspond à une vitesse relative dans le système Terre-Lune, ce n'est plus vrai dans le référentiel lunaire, puisque la vitesse orbitale de celle-ci s'ajoute ou se retranche selon l'orientation de l'approche, pour donner la vraie vitesse relative, $U_\infty$. Cela se traduit par un accroissement du flux d'impact autour de l'apex du satellite. Bon nombre d'auteurs [Neukum, G., 1983; Shoemaker et Wolfe, 1982; Zahnle et al., 1998, 2001] considèrent que l'orientation du vecteur $v_\infty^l$ est isotrope dans le référentiel du satellite, et en déduisent une expression analytique traduisant l'accroissement du flux selon la distance angulaire à l'apex. Notons
qu'ils ignorent de ce fait les autres asymétries potentielles (effet de latitude ou
effet face visible / face cachée). Il reste possible de traiter l'effet de latitude en
considérant que la Terre ne modifie pas substantiellement la distribution des i∞,
mais l'étude de l'effet face visible / face cachée, qui ne dépend que de la présence
de la Terre, nécessite une autre approche. Nous comparerons, au chapitre suivant,
ce genre d'approximations avec nos résultats, plus rigoureux.

5.5 Résultats

5.5.1 Système solaire actuel

Les distributions de probabilité marginalisées pour la vitesse d'approche sont représentés pour chaque cible sur la figure 5.8. Comme il a été mentionné au chapitre 3, les croiseurs planétaires possèdent en général une vitesse d'autant plus élevée qu'ils se trouvent à proximité du Soleil. Par conséquent, la vitesse d'approche moyenne est plus élevée pour Mercure qu'elle ne l'est pour Venus, et plus élevée pour Venus que pour Mars. Comme les planètes elles-mêmes ont une vitesse orbitale qui décroît avec la distance au Soleil, et que cette vitesse orbitale s'ajoute ou se retranche à celle des croiseurs selon la géométrie de l'approche, la distribution martienne est considérablement plus étroite que son homologue mercurienne. Les vitesses d'approche moyennes calculées sont : 42.2 km.s⁻¹ pour Mercure, 22.4 km.s⁻¹ pour Venus, 16.4 km.s⁻¹ pour la Terre, 19.2 km.s⁻¹ pour la Lune, et 9.1 km.s⁻¹ pour Mars. Les vitesses d'impact correspondantes sont 42.5, 25.2, 20.4, 19.4 et 10.6 km.s⁻¹, respectivement. Celles-ci sont d'autant plus élevées par rapport à la vitesse d'approche que la vitesse de libération à la surface de la planète est grande. Si la vitesse d'approche moyenne calculée pour la Lune est plus élevée qu'elle ne l'est pour la Terre – alors que la population de croiseurs est sensiblement la même : la Terre, on l'a vu, ne modifie pas significativement leur vitesse à une distance correspondant à la distance Terre-Lune – c'est parce que la probabilité d'impact dépend du pouvoir d'attraction de la cible, qui lui-même dépend de la vitesse d'approche. Les projectiles de faible vitesse sont beaucoup mieux captés par la Terre que par la Lune.

On rappelle que ces distributions sont construites en calculant les conditions
d'approche pour chaque orientation relative des orbites du projectile et de la cible.
La figure 5.9 compare la distribution des vitesses calculées de cette façon pour
Mercure et la distribution des inclinaisons calculées pour Mars, avec le résultat
<figure><img src="image_33.png" /><figcaption>FIG. 5.8 – Distribution des vitesses d’approche avec les planètes telluriques et la Lune.</figcaption></figure>

obtenu en faisant pour chaque croiseur la moyenne des conditions d'approche sur un cycle de précession (probabilité intégrée sur $\omega_0$ et $\Delta\Omega$). On constate que, si l'excentricité de la planète cible est conséquente, comme c'est le cas pour ces deux planètes, cette dernière approche peut conduire à une vue sensiblement biaisée des conditions du bombardement. Or il n'est pas rare qu'elle soit utilisée (par exemple dans Marchi et al. [2005] pour Mercure ou Ivanov et al. [2002] pour Mars).

La figure 5.10 montre la distribution des inclinaisons d'approche par rapport au plan orbital des planètes, $i_∞$, comparées au cas isotrope, où $i_∞$ est distribué avec une probabilité proportionnelle à $\cos i_∞$. Notons qu'ici, comme dans la suite de cette section, $i_∞$ est restreint aux valeurs positives, puisque les valeurs négatives sont équiprobables. On constate que Mars, et dans une moindre mesure Mercure, présentent une distribution qui favorise les hautes inclinaisons, tandis que la Terre et la Lune se verront préférentiellement bombardées par des projectiles approchant avec une inclinaison modérée. Venus ne se trouve pas non plus
<figure><img src="image_34.png" /><figcaption>FIG. 5.9 – Comparaison entre les distributions des conditions d'approche calculées en faisant ou non la moyenne sur un cycle de précession. A gauche. Distributions des vitesses d'approche calculées pour Mercure. A droite. Distributions des inclinaisons d'approche calculées pour Mars.</figcaption></figure>

dans le cas isotrope, mais d'une certaine façon, les hautes et basses inclinaisons se compensent.

Une telle différence de comportement dans le bombardement des planètes peut surprendre, en particulier dans le cas de Mars qui montre un fort excès de hautes inclinaisons. D'après les équations (5.7), (5.32), (5.30) et (5.31), l'inclinaison du vecteur relatif d'approche par rapport au plan orbital est donnée par

$$ \sin i_{\infty} = \frac{v}{U_{\infty}} \sin \alpha \sin i_{m}, $$

où *v* est la vitesse du projectile, *α* est l'angle entre **v** et le vecteur qui relie le Soleil à la planète, et *i*<sub>m</sub> est l'inclinaison entre les deux orbites. Dans le cas de Mars, dont la distance au Soleil peut atteindre 1.8 AU, la figure 3.5, qui donne la proportion d'objets en fonction de la valeur de leur périhélie, montre que la planète croise un grand nombre d'astéroïdes lorsque ceux-ci sont proches de leur périhélie, et qu'elle-même est proche de son aphélie (les projectiles sont beaucoup plus nombreux que lorsque Mars se situe à un autre point de son orbite). Qui plus est, la probabilité d'impact dans ce genre de configuration est particulièrement élevée. La valeur de *i*<sub>∞</sub> est par conséquent élevée pour deux raisons. Premièrement, *α* est proche de 90° et sin *α* est donc proche de son maximum. Deuxièmement, *U*<sub>∞</sub> est minimisée car l'angle entre les deux vecteurs vitesse est au plus faible. Cela explique également, en partie, le léger surplus de fortes inclinaisons d'approche pour Mercure, même si l'augmentation de la proportion d'objets avec la valeur de leur périhélie n'est pas aussi importante dans cette zone (~ 0.4 AU). La figure 3.9
<figure><img src="image_35.png" /><figcaption>FIG. 5.10 – Distribution des inclinaisons d'approche avec les planètes telluriques, comparées à l'isotropie.</figcaption></figure>

montre également que, pour Mercure, les objets possèdent en moyenne une plus grande inclinaison. Dans la mesure où $i_m$ est à peu près égale à l'inclinaison des projectiles (l'inclinaison de Mercure est de 7°), l'inclinaison d'approche $i_∞$ sera en moyenne plus élevée que pour Venus et la Terre.

On appelle $R$ le flux relatif de projectiles par rapport à la Lune, par unité de surface. En considérant que les éléments orbitaux des projectiles ne sont pas corrélés à leur taille, $R$ compare le bombardement entre les planètes pour un certain diamètre de projectiles, et est valide à tout diamètre. Pour la Terre, $R = 1.58$, ce qui est proche de la valeur déterminée par Ivanov [2006] à partir des observations télescopiques à cette date. Pour Mercure, $R = 1.81$, comme déterminé par Marchi et al. [2005] (d'après Ivanov et Hartmann, W. K. [2007]). On calcule $R = 1.73$ pour Venus, ce qui correspond à un bombardement relatif à la Terre de 1.09, qui diffère de la valeur de 0.95 estimée par Shoemaker et al. [1991] à partir d'observations limitées de croiseurs terrestres et vénusiens. Mars présente le plus fort bombardement par unité de surface, puisque $R = 3.13$. A
partir du flux terrestre estimé au chapitre 3, le bombardement lunaire absolu est
$\Phi(d > 1\text{km}) = 2.1 \times 10^{-15} \text{ km}^{-2} \text{ an}^{-1}$. Il correspond à environ 950 NEO tels que
$d > 1\text{km}$.

<figure><img src="image_36.png" /><figcaption>FIG. 5.11 – Distribution de probabilité des conditions d’approche $U_∞$ et $i_∞$ pour les planètes telluriques et la Lune. Eléments orbitaux actuels. Valeurs normalisées par la probabilité maximale pour chaque planète.</figcaption></figure>

La figure 5.11 compare les distributions de probabilité en fonction de $U_∞$ et $i_∞$ pour les planètes telluriques et la Lune. En général, la rotation des planètes autour de leur axe rend la détermination de la distribution des $j_∞$ inutile pour le calcul du bombardement en fonction des cordonnées « géographiques », puisque celui-ci ne dépendra pas de la longitude. Deux exceptions cependant : la Lune, en rotation synchrone autour de la Terre, et la Terre, pour laquelle nous souhaitons déterminer l'intensité du bombardement en fonction de l'heure locale.

La figure 5.12 donne l'ensemble des distributions d'approche qui caractérisent
le bombardement terrestre. On constate notamment que la distribution des $j_∞$ dé-
pend fortement de la vitesse d'approche considérée. La distribution des coordon-
<figure><img src="image_37.png" /><figcaption>FIG. 5.12 – Distribution de probabilité des conditions d’approche terrestres. A gauche. Probabilité de chaque couple $(U_\infty, i_\infty)$. Au milieu. Probabilité de chaque couple $(U_\infty, j_\infty)$. A droite. Probabilité de chaque couple $(i_\infty, j_\infty)$.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>FIG. 5.13 – Comparaison de la distribution terrestre des coordonnées du radiant $(i_\infty, j_\infty)$ avec les données radar de Campbell-Brown [2008]. La figure de droite est directement extraite de Campbell-Brown [2008].</figcaption></figure>

nées du radiant, $i_\infty$ et $j_\infty$, est comparée aux données radar de Campbell-Brown [2008] en figure 5.13. Celles-ci ont été obtenues à partir des observations effectuées par CMOR (Canadian Meteor Orbit Radar) sur une durée de 5 ans. Les observations concernent des projectiles se situant dans une gamme de taille millimétrique. L'excellent accord entre notre simulation et les résultats de Campbell-Brown [2008] suggèrent fortement que les météoroïdes reproduisent le comportement orbital de leurs astéroïdes parents, plus gros.

La figure 5.14 donne l'ensemble des distributions d'approche qui caractérisent le bombardement lunaire, et l'on constate d'ores et déjà que le flux de projectiles est concentré à l'apex.

Les valeurs moyennes des grandeurs qui caractérisent le bombardement sont listées en table 5.1.
<figure><img src="image_39.png" /></figure>

FIG. 5.14 – Distribution de probabilité des conditions d'approche lunaires. Distance Terre-Lune actuelle. A gauche. Probabilité de chaque couple ($U_\infty, i_\infty$). Au milieu. Probabilité de chaque couple ($U_\infty, j_\infty$). A droite. Probabilité de chaque couple ($i_\infty, j_\infty$).

<table><caption>TAB. 5.1 – Conditions moyennes du bombardement des planètes telluriques.</caption><thead><tr><th>Planète</th><th>Mercure</th><th>Vénus</th><th>Terre</th><th>Lune</th><th>Mars</th><th>Mars (moyenne séculaire)</th></tr></thead><tbody><tr><td>Vitesse d'approche (km.s<sup>-1</sup>)</td><td>42.2</td><td>22.4</td><td>16.4</td><td>19.2</td><td>9.1</td><td>9.3</td></tr><tr><td>Déflexion gravitationnelle Δ</td><td>0.02</td><td>0.96</td><td>1.88</td><td>0.04</td><td>0.58</td><td>0.58</td></tr><tr><td>Vitesse d'impact (km.s<sup>-1</sup>)</td><td>42.5</td><td>25.2</td><td>20.4</td><td>19.4</td><td>10.6</td><td>10.8</td></tr><tr><td>Inclinaison d'approche i<sub>∞</sub> (°)</td><td>34.1</td><td>32.7</td><td>30.2</td><td>30.0</td><td>41.5</td><td>41.0</td></tr><tr><td>Bombardement relatif à la Lune R</td><td>1.81</td><td>1.73</td><td>1.58</td><td>1.00</td><td>3.13</td><td>2.75</td></tr></tbody></table>

### 5.5.2 Influence des variations séculaires

Comme on a pu le voir au chapitre 4, les éléments orbitaux des planètes sont soumis à des variations séculaires. En incluant ces variations dans le calcul des probabilités de collision détaillées tout au long de ce chapitre, et en supposant que la distribution orbitale des projectiles est restée inchangée, on quantifie dans quelle mesure la distribution des conditions d'approche et le bombardement relatif entre les planètes sont affectés.

Nous constatons que le flux d'impacts absolu pour le système Terre-Lune, Venus et Mercure ne diffère jamais de plus de 2%, et ce, même si les variations orbitales de Mercure, en particulier, sont grandes. La distribution des conditions d'approche sur le long terme est quand à elle difficilement distinguable de son homologue actuel. Ce dernier point est aussi vérifié pour Mars, comme le montre la figure 5.15. Par contre, dans la mesure où Mars se trouve à un moment de son
cycle séculaire où elle possède une forte excentricité, elle s'aventure actuellement dans des zones particulièrement proches de la ceinture d'astéroïdes, et donc très peuplées. Sur le long terme (3 milliards d'années), le bombardement absolu est réduit et l'on trouve $R = 2.75$. Cette valeur est équivalente à celle déterminée par Ivanov et al. [2002].

<figure><img src="image_40.png" /><figcaption>FIG. 5.15 – Comparaison des distributions de probabilité ($U_∞, j_∞$) entre l'actuelle Mars et Mars sur l'ensemble de son histoire orbitale.</figcaption></figure>

On l'a vu, la distance Terre-Lune est fonction du temps, même si la relation qui lie l'un à l'autre est encore mal contrainte. Si la distribution des $U_∞$ et $j_∞$ ne semble pas varier pour une distance Terre-Lune comprise entre 20 et 60 rayons terrestres, ce n'est pas le cas pour $j_∞$ : la vitesse croissante du satellite, à mesure qu'il s'approche de la planète, concentre d'autant plus le flux d'impact à l'apex (voir le chapitre suivant). On le vérifie sur la figure 5.16. Le bombardement relatif entre la Terre et la Lune change quand à lui tellement modérément (pour la Terre, $R$ est réduit d'environ 1 % entre 60 et 20 rayons terrestres) que nous considérons qu'il peut être tenu pour constant dans cette gamme de distances Terre-Lune.

Sous l'hypothèse que la distribution orbitale des astéroïdes et comètes est restée inchangée dans les derniers milliards d'années, on peut donc considérer que le bombardement des planètes, à l'exception de Mars, n'est pas sensible (en moyenne) aux variations séculaires, et que le bombardement actuel peut être pris pour représentatif du bombardement moyen sur le long-terme. Pour Mars,
<figure><img src="image_41.png" /><figcaption>FIG. 5.16 – Distribution de probabilité de l'angle $j_∞$ pour la Lune, pour différentes distances Terre-Lune.</figcaption></figure>

en revanche, si la distribution des conditions d'approche n'est pas perturbée, le flux absolu est réduit par les variations d'excentricité. Enfin, il faudra considérer les variations d'obliquité dans la détermination des variations géographiques du flux d'impacts (chapitre suivant). Nous verrons encore une fois que seule Mars verra les variations spatiales du flux d'impact grandement réduites par de forts changements d'obliquité. Les variations d'obliquité très modérées de l'ensemble des autres planètes n'auront pas d'incidence. Pour cette raison, nous distinguerons la Mars actuelle et la Mars sur le long terme, alors que nous ne ferons pas cette distinction pour les autres planètes. Pour la Lune, il sera toujours précisé quelle distance Terre-Lune est considérée.
# Chapitre 6

## Simuler le bombardement

Dans ce chapitre, nous nous attachons à décrire les variations du flux d'impacts à la surface des planètes et des lunes, en fonction de la latitude et de la longitude géographiques. La plupart des travaux qui traitent du bombardement des planètes considèrent que celui-ci est isotrope. Nous verrons dans une première partie ce que cela implique comme approximations, et comment simuler un tel bombardement. En réalité, si la source de projectiles n'est pas elle-même isotrope, ou si l'on traite du cas d'un satellite en rotation synchrone, des effets de latitude et de longitude sont à attendre. La deuxième partie les détaille de manière schématique. Connaissant les conditions d'approche avec la cible, nous développons dans la troisième partie une expression analytique du flux d'impacts en fonction de la position à la surface. Nous discutons dans une quatrième partie des paramètres qui influent sur l'amplitude des variations spatiales, avant de présenter les résultats obtenus pour les planètes telluriques et la Lune. Dans une sixième et dernière partie, nous comparons les résultats obtenus pour la Lune avec les données sismiques des stations Apollo.

### 6.1 Hypothèse d'isotropie

Considérer que le bombardement est isotrope revient à dire que le flux d'impacts (nombre de projectiles qui heurtent la planète par unité de temps et de surface) est constant quelle que soit la valeur du couple $(\lambda, \varphi)$, la latitude et la longitude géographiques.

En un point donné, l'élément surfacique s'exprime $dS = R_0^2 \cos \lambda \, d\lambda \, d\varphi$. Une densité d'impacts constante est donc vérifiée si la probabilité d'impact à une latitude donnée est proportionnelle à $\cos \lambda$, et si la longitude du site d'impact est uniformément distribuée entre 0 et $2\pi$. Outre la position de l'impact, il est important de connaître la valeur de l'angle d'incidence $\theta$ (que l'on choisit ici de mesurer par rapport à l'horizontale). En effet, cet angle joue un rôle important
dans la détermination de la taille du cratère d'impact (chapitre suivant). Il peut
aisément être montré [Pierazzo et Melosh, 2000] que cette probabilité est sin 2θ.
Notons que cela implique que l'angle d'incidence moyen à la surface de la planète
est π/4.

Echantillonner une variable revient à identifier sa probabilité cumulée (et normalisée) avec un nombre aléatoire compris entre 0 et 1. Pour une variable donnée $v$, on appelle « probabilité cumulée et normalisée » la probabilité d'avoir une valeur de $v'$ supérieure (ou inférieure) à $v$, donnée par $P(> v) = \text{Const.} \times \int_v^{v_{\text{max}}} P(v')dv'$, et normalisée telle que $\int_v P(> v)dv = 1$.

On peut ainsi simuler un bombardement isotrope en tirant pour N impacts
la position et l'angle d'incidence tels que

$$
\begin{equation}
\begin{cases}
\lambda = \sin^{-1}(2x_1 - 1) \\
\varphi = 2\pi x_2 \\
\theta = \frac{1}{2} \cos^{-1}(2x_3 - 1),
\end{cases}
\tag{6.1}
\end{equation}
$$

où les $x_i \in ]0, 1[$ sont des nombres aléatoires uniformément distribués.

Considérer que le bombardement est isotrope revient à dire que les distribu-
tions des angles d'approche $i_\infty$ et $j_\infty$ sont indépendantes, et que la probabilité du
premier est proportionnelle à $\cos i_\infty$, tandis que $j_\infty$ est uniformément distribué
entre 0 et $2\pi$. Nous avons vu au chapitre précédent qu'il n'en est rien. Une autre
possibilité pour justifier un bombardement isotrope est que la masse de la planète
est suffisamment grande pour homogénéiser par gravité les conditions d'impact
à sa surface, quelle que soit la distribution des conditions d'approche. Cela sera
discuté dans les sections suivantes.

6.2 Cas réel : asymétries

6.2.1 Planètes et lunes : effet de latitude

On illustre schématiquement l'effet de latitude provoqué par une source ani-
sotrope à l'aide de la figure 6.1. Sur ce schéma, l'inclinaison d'approche i∞ est
nulle, et les projectiles possèdent une vitesse infinie. Ils ne sont par conséquent
pas perturbés par la masse de la planète et parcourent un trajet rectiligne (un
cas équivalent reviendrait à poser une vitesse d'approche finie et une masse pla-
nétaire nulle). Par simple projection géométrique sur la sphère planétaire, on
déduit directement que le flux relatif d'impacts sera distribué selon le cosinus de
la latitude, c'est-à-dire maximal à l'équateur et nul au pôle.
Bien sûr, les conditions réelles sont autrement plus complexes : la source n’est jamais aussi anisotrope, et la masse de la planète, en courbant la trajectoire des projectiles par gravité, réduit l’amplitude des variations. Néanmoins, on doit par principe s’attendre à ce que le flux présente une dépendance en latitude sur l’ensemble des planètes et des lunes. Dans le cas des planètes, ce sera d’ailleurs le seul effet observé, étant donné que leur rotation propre annule toute dépendance longitudinale. On ne peut pas en dire autant des lunes en rotation synchrone, comme il est détaillé par la suite.

<figure><img src="image_42.png" /><figcaption>FIG. 6.1 – Illustration schématique de l’effet en latitude.</figcaption></figure>

### 6.2.2 Satellites en rotation synchrone : effet de longitude

#### 6.2.2.1 Effet apex / antapex

La géométrie du problème est schématisée en figure 6.2. On observe que, quelle que soit la direction d'où proviennent les projectiles, le flux d'impacts sera plus grand à l'apex (le centre de l'hémisphère qui est en tête du mouvement) qu'à l'antapex (le centre de l'hémisphère opposé). L'apex et l'antapex jouent respectivement le rôle de proue et de poupe. Autour de l'apex, la vitesse lunaire s'additionne à la vitesse des projectiles, alors qu'elle se soustrait autour de l'antapex. Le nombre d'impacts par unité de temps est donc nécessairement plus élevé autour du premier que du second. L'amplitude de cet effet dépend bien sûr de la différence de vitesse qui existe entre la lune et les projectiles, entre autres choses.

#### 6.2.2.2 Effet face visible / face cachée

Considérons un flux de projectiles parallèle au plan orbital lunaire ($i_∞ = 0$). La planète centrale, si sa masse est suffisante, concentre les projectiles (comme
<figure><img src="image_43.png" /><figcaption>FIG. 6.2 – Illustration schématique de l’effet apex / antapex.</figcaption></figure>

le ferait une lentille avec la lumière) le long d'une ligne « focale ». Lorsque la lune croise cette ligne, le bombardement s'intensifie. Ce phénomène est illustré en figure 6.3. On s'attend donc à observer une asymétrie face visible / face cachée. Pour cela, il faut cependant (entre autres) que le vecteur vitesse d'approche ne soit pas trop éloigné de la colinéarité avec le vecteur planète-lune.

<figure><img src="image_44.png" /><figcaption>FIG. 6.3 – Illustration schématique de l’effet face visible / face cachée.</figcaption></figure>

## 6.3 Expression analytique

La probabilité d'impact pour un type d'approche particulier peut s'écrire

$$p(\mathbf{U}_{\infty}) = \int_{b=0}^{\tau} \int_{\delta=0}^{2\pi} f(\mathbf{U}_{\infty}) \, db \, d\delta = f(\mathbf{U}_{\infty})\sigma, \qquad (6.2)$$
où $f$, constant sur $(\sigma)$, est le flux surfacique de projectiles qui traversent $(\sigma)$, c'est-à-dire leur nombre par unité de temps et de surface.

<figure><img src="image_45.png" /><figcaption>FIG. 6.4 – Relation entre le paramètre d'impact, $b$, et la distance angulaire au radiant, $\beta$.</figcaption></figure>

Sur la planète, le flux d'impacts global par unité de surface, généré par une certaine triade $(U_\infty, i_\infty, j_\infty)$, est donné par

$$ \Phi(\mathbf{U}_\infty) = \frac{p(\mathbf{U}_\infty)}{S_0} = f \frac{\sigma}{S_0}, \qquad (6.3) $$

où l'indice 0 indique que le flux est moyenné sur toute la surface $S_0$ de la planète, telle que $S_0 = 4\pi R_0^2$. Comme indiqué sur la figure 6.4, une certaine valeur de $b$ implique une certaine valeur de $\beta$, où $\beta$ est l'angle que fait le vecteur position du site d'impact avec le radiant, dans un repère centré sur la planète. On voit sur la figure que le problème est symétrique autour du radiant. Si l'on note $d\sigma$ l'aire du l'anneau d'épaisseur $db$ sur $(\sigma)$, telle que $d\sigma = 2\pi b \, db$, alors l'aire de l'anneau correspondant sur la planète s'exprime $dS = 2\pi R_0^2 \sin\beta \, d\beta$. Le flux surfacique sur la planète à une certaine distance angulaire $\beta$ s'exprime alors

$$ \Phi(\beta, \mathbf{U}_\infty) = f \frac{d\sigma}{dS}, \qquad (6.4) $$

et le flux relatif, normalisé par le flux global, est

$$ \phi(\beta, U_\infty) = \frac{\Phi(\beta, U_\infty)}{\Phi(U_\infty)} = \frac{S_0}{\sigma} \frac{d\sigma}{dS} = \frac{4}{\tau^2} \frac{b \, db}{\sin\beta \, d\beta}. \qquad (6.5) $$

Il nous faut à présent exprimer $b \, db$ en fonction de $\beta$ si l'on veut définir le flux d'impact comme une fonction de la seule position sur la planète. La distance
angulaire $\beta$ est donnée par $\beta = f_i - f_\infty$, où $f_i$ est l'anomalie vraie à l'impact, et $f_\infty$ est l'anomalie vraie à l'infini (annexe A). On a donc

$$ \left\{ \begin{aligned} \cos \beta &= C_i C_\infty + S_i S_\infty , \\ \sin \beta &= S_i C_\infty - C_i S_\infty . \end{aligned} \right. \qquad (6.6) $$

En faisant usage des équations (A.7)–(A.10), (A.11) et (6.6), on vérifiera que le paramètre d'impact $b$ s'exprime en fonction de la distance angulaire d'impact au radiant $\beta$ comme

$$ b = R_0 \frac{\tan \beta/2}{1 + \tan^2 \beta/2} \left( 1 + \sqrt{1 + \Delta(1 + \tan^2 \beta/2)} \right), \qquad (6.7) $$

et que la valeur maximale de $\beta$, qui correspond à $b = \tau$ est donnée par

$$ \cos \beta_{max} = -\frac{\Delta}{2 + \Delta}. \qquad (6.8) $$

On constate que $\lim_{\Delta \to 0} \beta_{max} = \pi/2$ et $\lim_{\Delta \to \infty} \beta_{max} = \pi$; quand l'attraction gravitationnelle est infinie au regard de l'énergie cinétique des projectiles, le flux d'impacts englobe la totalité de la sphère, tandis que dans le cas opposé, où la vitesse des projectiles est très grande, seul l'hémisphère apparent dans la direction du radiant est soumis au bombardement.

Le flux relatif en $\beta$, par unité de surface, une fois $b$ db exprimé et inséré dans l'équation (6.5), est

$$ \phi(\beta, U_{\infty}) = \frac{1 + \cos \beta}{2(1 + \Delta)} \frac{1 + \Xi}{\Xi} \left[ \Xi^2 + \Xi - (1 + \Xi) \frac{1 - \cos \beta}{1 + \cos \beta} \right], \qquad (6.9) $$

avec

$$ \Xi = \sqrt{1 + \Delta + \Delta \frac{1 - \cos \beta}{1 + \cos \beta}}. \qquad (6.10) $$

En l'absence de force gravitationnelle (planète « non-massique »), $\Delta = 0$ et le flux est simplement $\phi(\beta) = 4 \cos \beta$. À l'inverse, $\lim_{\Delta \to \infty} \phi(\beta, U_\infty) \sim \frac{1+\cos \beta}{2\Delta} \frac{\sqrt{\Delta}}{\sqrt{\Delta}} \left[ \Delta \left( 1 + \frac{1-\cos \beta}{1+\cos \beta} \right) \right] = 1$; le flux relatif d'impact est dans ce cas le même sur toute la surface de la planète.

Il nous faut maintenant relier $\beta$ à la position absolue sur la planète. Dans un premier temps, nous faisons intervenir la latitude et la longitude orbitales, $\lambda_0$ et $\varphi_0$. Dans le repère orbital (OXYZ), $\lambda_0$ est l'angle qui sépare le vecteur position du site d'impact et le plan orbital, et $\varphi_0$ est l'angle entre la projection de ce vecteur sur (OXY) et l'axe (OX). La différence avec $i_\infty$ et $j_\infty$ est qu'il s'agit des réelles cordonnées de l'impact, et non des cordonnées du radiant.
Dans (OXYZ), le vecteur unitaire dans la direction du radiant, **rad**, s'écrit

$$ \mathbf{rad} = \begin{pmatrix} C_{i\infty} & C_{j\infty} \\ C_{i\infty} & S_{j\infty} \\ S_{i\infty} & C_{j\infty} \end{pmatrix}, \qquad (6.11) $$

où $C_{i\infty}$ et $S_{i\infty}$ sont le cosinus et le sinus de $i\infty$ ; $C_{j\infty}$ et $S_{j\infty}$ le cosinus et le sinus de $j\infty$. En calculant le produit scalaire du vecteur unitaire correspondant à la position d'un impact quelconque avec **rad**, on en déduit

$$ \cos \beta = (\cos \lambda_0 \cos \varphi_0, \cos \lambda_0 \sin \varphi_0, \sin \lambda_0) \times \begin{pmatrix} C_{i\infty} & C_{j\infty} \\ C_{i\infty} & S_{j\infty} \\ S_{i\infty} & C_{j\infty} \end{pmatrix}, \qquad (6.12) $$

et il suffit de choisir la latitude et la longitude orbitales pour en déduire $\cos \beta$. Si $\beta > \beta_{max}$ (Éq. (6.8)), $\phi$ est nul, sinon, il est donné par l'équation (6.9). On écrit

$$ \phi(\lambda_0, \varphi_0, \mathbf{U}_\infty) = \phi(\beta(\lambda_0, \varphi_0, i_\infty, j_\infty), U_\infty) \quad \begin{array}{l} \text{si } \beta(\lambda_0, \varphi_0, i_\infty, j_\infty) \le \beta_{max}, \\ \text{si } \beta(\lambda_0, \varphi_0, i_\infty, j_\infty) > \beta_{max} \end{array} \qquad (6.13) $$

L'angle d'incidence mesuré à partir de l'horizontale, $\theta$, est donné par (cf. équation (A.18)),

$$ \theta(\lambda_0, \varphi_0, \mathbf{U}_\infty) = \cos^{-1} \left( \frac{1 + \Xi}{2\sqrt{1 + \Delta}} \sqrt{1 - \cos^2 \beta} \right) . \qquad (6.14) $$

tandis que la vitesse d'impact est (cf. équation (A.16))

$$ u(\lambda_0, \varphi_0, \mathbf{U}_\infty) = u(U_\infty) = U_\infty \sqrt{1 + \Delta} . \qquad (6.15) $$

On pourra vérifier en intégrant l'équation (6.14) entre $\beta = 0$ et $\beta = \beta_{max}$ que l'angle d'impact moyen sur l'ensemble de la surface bombardée est toujours égal à $\pi/4$, indépendamment des conditions d'approche et de la masse de la planète.

En intégrant sur l'ensemble des conditions d'approche, on obtient

$$ \phi(\lambda_0, \varphi_0) = \iint \int \phi(\lambda_0, \varphi_0, \mathbf{U}_\infty) p(\mathbf{U}_\infty) dU_\infty d\lambda_\infty d\varphi_\infty . \qquad (6.16) $$

Afin d'exprimer le flux dans un repère géographique et non plus orbital, il est nécessaire de prendre en compte une éventuelle obliquité $\epsilon$. On appelle (Oxyz) le repère géographique, $\lambda$ et $\varphi$ la latitude et la longitude associées. L'axe de spin précessse, et nous souhaitons obtenir le flux sur la moyenne de cette précession. La position de cet axe dans (OXYZ) est donnée par $\lambda_\epsilon = \pi/2 - \epsilon$ et la longitude $\varphi_\epsilon$ varie uniformément entre 0 et $2\pi$. La géométrie du problème est représentée
<figure><img src="image_46.png" /><figcaption>FIG. 6.5 – Passage du repère orbital au repère géographique par rotation d'un angle $\epsilon$ autour de l'axe (A).</figcaption></figure>

en figure 6.5. Si $\epsilon = 0$, (Oxyz) est confondu avec (OXYZ), et $(\lambda, \varphi) = (\lambda_0, \varphi_0)$. Le passage de (OXYZ) à (Oxyz) se fait par rotation d'un angle $\epsilon$ autour de l'axe (A) définie par $\varphi_\epsilon$ comme indiqué sur la figure. La longitude $\varphi$ ainsi obtenue correspond au cas d'une rotation synchrone, par exemple celui de la Lune. L'asynchronicité de la rotation et de la révolution des planètes rend le flux indépendant de $\varphi$, et sa définition triviale.

La matrice de passage de (OXYZ) à (Oxyz) s'écrit

$$ M = \begin{bmatrix} C_{\varphi\epsilon}^2 + S_{\varphi\epsilon}^2 C_{\epsilon} & C_{\varphi\epsilon} S_{\varphi\epsilon} (1 - C_{\epsilon}) & -S_{\varphi\epsilon} S_{\epsilon} \\ C_{\varphi\epsilon} S_{\varphi\epsilon} (1 - C_{\epsilon}) & S_{\varphi\epsilon}^2 + C_{\varphi\epsilon}^2 C_{\epsilon} & C_{\varphi\epsilon} S_{\epsilon} \\ S_{\varphi\epsilon} S_{\epsilon} & -C_{\varphi\epsilon} S_{\epsilon} & C_{\epsilon} \end{bmatrix}, \qquad (6.17) $$

où C et S désignent le cosinus et le sinus de la variable indiquée en indice. Dans ces conditions, en notant $W = M^{-1}$ la matrice inverse, $\lambda_0(\lambda, \varphi, \epsilon, \varphi_\epsilon)$ et $\varphi_0(\lambda, \varphi, \epsilon, \varphi_\epsilon)$ sont données comme

$$ \sin \lambda_0 = W_{31} \cos \lambda \cos \varphi + W_{32} \cos \lambda \sin \varphi + W_{33} \sin \lambda, \qquad (6.18) $$

et

$$ \begin{align} \cos \varphi_0 &= \frac{W_{11} \cos \lambda \cos \varphi + W_{12} \cos \lambda \sin \varphi + W_{13} \sin \lambda}{\sqrt{1 - \sin^2 \lambda_0}}, \notag \\ \sin \varphi_0 &= \frac{W_{21} \cos \lambda \cos \varphi + W_{22} \cos \lambda \sin \varphi + W_{23} \sin \lambda}{\sqrt{1 - \sin^2 \lambda_0}}. \tag{6.19} \end{align} $$

Le flux est finalement donné en fonction de la position géographique comme

$$ \phi(\lambda, \varphi) = \frac{1}{2\pi} \int_{0}^{2\pi} \phi(\lambda_0(\lambda, \varphi, \epsilon, \varphi_\epsilon), \varphi_0(\lambda, \varphi, \epsilon, \varphi_\epsilon)) \, d\varphi_\epsilon. \qquad (6.20) $$
Dans le cas des planètes, où la définition de $\varphi$ est artificielle, le flux, uniquement fonction de la latitude, est

$$ \phi(\lambda) = \frac{1}{(2\pi)^2} \int_{0}^{2\pi} \int_{0}^{2\pi} \phi(\lambda_0(\lambda, \varphi = 0, \epsilon, \varphi_\epsilon), \varphi_0) \, d\varphi_0 \, d\varphi_\epsilon . \quad (6.21) $$

Le flux absolu à une position donnée s'exprime à partir du flux total (donné au chapitre 6 par l'équation (5.1)) comme

$$ \Phi(> d, \lambda, \varphi) = \phi(\lambda, \varphi)\Phi(> d) . \quad (6.22) $$

Puisque $\phi$ est le flux d'impacts relatif, on rappelle que $\frac{1}{2\pi^2} \int_{\lambda} \int_{\varphi} \phi(\lambda, \varphi) = 1$.

Notons enfin que, sur une planète, le flux est homogène quand la distribution des inclinaisons par rapport au plan orbital vérifie l'isotropie, indépendamment de la vitesse d'approche. En d'autres termes, $\phi(\lambda) = 1 \forall \lambda$, quand $p(U_\infty, i_\infty) = p(U_\infty) p(i_\infty) = p(U_\infty) \frac{1}{2} \cos i_\infty$. C'est, avec le cas $\Delta = 0$, la deuxième cause possible à l'absence de variations latitudinales du flux d'impacts à la surface d'une planète. Nous verrons dans quelle mesure le fait que ces conditions ne sont pas respectées produit une dépendance spatiale.

## 6.4 Discussion

Dans cette section, nous explorons l'influence des paramètres qui gouvernent le bombardement, en particulier la variable $\Delta = \frac{GM_0}{R_0U_\infty^2}$. Nous avons pu observer précédemment que cette grandeur est un outil clef dans la prédiction de possibles variations spatiales du flux d'impacts. Lorsque $\Delta$ est grand, le flux est rendu homogène par la gravité de la planète. Au contraire, lorsque $\Delta$ est petit, les variations spatiales du flux sont gouvernées par les conditions d'approche.

La figure 6.6 présente quelques exemples de variations en latitude induites par différentes valeurs de $\Delta$. Nous considérons deux géométries opposées : (1) la planète est soumise à un bombardement d'objets qui s'approchent parallèlement à son plan équatorial ($i_\infty = 0$) (2) le bombardement est perpendiculaire au plan équatorial ($i_\infty = \pm \pi/2$). L'obliquité de la cible est nulle. Quand la vitesse d'approche est grande par rapport à la masse de la planète, $\Delta$ est proche de zéro, et la trajectoire des projectiles n'est quasiment pas influencée par la gravité de la planète. Dans ce cas, le flux d'impacts présente ses plus fortes variations : on a vu que dans le cas limite où $\Delta = 0$, $\phi = 4 \cos \beta$. Cela correspond à $\phi = \frac{4}{\pi} \cos \lambda$ dans le cas (1), et $\phi = 2 \sin \lambda$ dans le cas (2). On note que les variations spatiales provoquées par des projectiles fortement inclinés sont les plus importantes. A
<figure><img src="image_47.png" /><figcaption>FIG. 6.6 – Variations latitudinales du flux d'impacts pour différentes valeurs de $\Delta$ et $i_{\infty}$.</figcaption></figure>

mesure que la masse de la planète croît (ou que la vitesse des projectiles décroît), $\Delta$ augmente, et les trajectoires sont de plus en plus déviées par la planète. En conséquence, le flux d'impacts devient de plus en plus homogène. Dans le cas limité où $\Delta = \infty$, le flux sera identique en tout point de la planète, quelle que soit la valeur de $i_{\infty}$. Notons qu'en pratique, il suffit que $\Delta$ vaille environ 5 pour que les variations soient très faibles (inférieures au % sur l'ensemble de la surface.)

<figure><img src="image_48.png" /><figcaption>FIG. 6.7 – Variations latitudinales de l'angle moyen d'impact pour différentes valeurs de $\Delta$ et $i_{\infty}$.</figcaption></figure>
Pour les mêmes raisons, l'angle d'impact moyen présente les plus fortes variations latitudinales quand $\Delta$ est faible, comme le montre la figure 6.7. Il sera plus élevé à l'équateur dans le cas (1), et plus élevé aux pôles dans le cas (2). L'amplitude des variations se réduit à mesure que $\Delta$ augmente, jusqu'à ce que l'angle d'impact soit égal à 45° en tout point. On note qu'on retrouve ce dernier résultat lorsque la distribution de projectiles est isotrope.

L'hypothèse répandue selon laquelle le bombardement est isotrope est correcte, on l'a vu, si la source de projectile est elle-même isotrope, ou si $\Delta = \infty$. Pourtant, la première condition n'est en général pas respectée pour les planètes telluriques (voir les résultats du chapitre précédent). La deuxième condition non plus, comme le montrent les valeurs moyennes de $\Delta$ listées en table 5.1.

<figure><img src="image_49.png" /><figcaption>FIG. 6.8 – Variations du flux d'impacts en fonction de la distance angulaire à l'apex, pour différentes valeurs de $\Delta$ et $v_l$. La source de projectiles qui bombarde le système planète-lune est isotrope, et la vitesse d'approche est de 15 km/s.</figcaption></figure>

Pour finir, la figure 6.8 représente les variations du flux d'impacts en fonction de la distance angulaire à l'apex, $\gamma$, dans le cas d'un satellite en rotation synchrone, pour différentes valeurs de $\Delta$ et différentes vitesses orbitales $v_l$. Toutes les courbes ont été calculées pour une source isotrope et une vitesse d'approche de 15 km/s. On constate bien sûr que, plus la vitesse du satellite est grande, plus l'effet est important. D'autre part, une valeur croissante de $\Delta$ (donc de la masse du satellite ici) réduit l'amplitude des variations. La figure 6.9 poursuit la comparaison en montrant les résultats qui correspondent au conditions d'approche lunaires, avec certaines approximations – ou non – et anticipe les résultats de la section
<figure><img src="image_50.png" /><figcaption>FIG. 6.9 – Variations du flux d'impacts en fonction de la distance angulaire à l'apex. Vitesse d'approche moyenne de 19 km/s et vitesse orbitale lunaire de 1 km/s.</figcaption></figure>

suivante. Le bombardement du système Terre-Lune est caractérisé par une distribution de probabilité des conditions d'approche complexes. Les résultats que nous obtenons sont représentés par la courbe rouge. La courbe noire correspond à une simulation qui considère que la distribution de projectiles vérifie l'isotropie, et que chaque projectile possède une vitesse d'approche avec le système qui correspond à la vitesse moyenne du premier cas (soit $U' \simeq 19$ km/s). On constate que ces approximations mènent à un comportement différent du flux d'impacts, en particulier, la courbe est plus symétrique, et sous-estime le bombardement autour de l'apex. Nous comparons cette dernière aux expressions proposées par Zahnle et al. [1998] et Zahnle et al. [2001]. De manière analytique, Zahnle et al. [1998] proposent que le flux d'impacts s'écrive :

$$ \phi(\gamma) = \left( 1 + \frac{v_l}{\sqrt{2v_l^2 + U'^2}} \cos \gamma \right). \qquad (6.23) $$

Dans un deuxième temps, en déterminant empiriquement une solution analytique que reproduit au mieux les résultats de leurs simulations, Zahnle et al. [2001] donnent

$$ \phi(\gamma) = \left( 1 + \frac{v_l}{\sqrt{2v_l^2 + U'^2}} \cos \gamma \right)^2 . \qquad (6.24) $$
Ces deux relations sont basées sur l'hypothèse que la masse du satellite est nulle, et que le bombardement du système Terre-Lune est isotrope. On constate que la deuxième expression reproduit de manière très satisfaisante notre simulation qui fait l'hypothèse de l'isotropie (courbe noire), mais pas notre simulation plus rigoureuse (courbe rouge). Notons enfin que la masse de la Lune ne joue pas un rôle significatif à ces vitesses (la courbe rouge, recalculée pour une masse lunaire nulle, est quasi identique.)

## 6.5 Résultats

<figure><img src="image_51.png" /><figcaption>FIG. 6.10 – Variations latitudinales du flux d'impacts à la surface des planètes telluriques.</figcaption></figure>

La figure 6.10 présente les variations du flux d'impacts en fonction de la latitude, calculées à partir de la distribution des conditions d'approche déterminées au chapitre 5 pour les planètes telluriques et la Lune. Le ratio pôle / équateur est calculé à 1.05, 1.00, 0.96, 0.90 et 1.26 pour Mercure, Venus, la Terre, la Lune et Mars, respectivement. Mercure et Mars, en raison de la grande proportion d'ap-
proches de forte inclinaison, subissent un bombardement plus intense aux pôles. Qui plus est, $\Delta$ est en moyenne suffisamment faible pour ne pas homogénéiser le flux. Venus ne présente aucune variation significative, d'une part parce que la source est presque isotrope, d'autre part parce que $\Delta$ est relativement élevé. La Terre, quant à elle, entre en collision avec des projectiles préférentiellement peu inclinés, ce qui entraîne un très léger excès d'impacts à l'équateur. La distribution des $i_{\infty}$ est quasiment identique pour la Lune, mais les faibles valeurs de $\Delta$ par rapport à la Terre (vitesses d'approche plus élevées, masse beaucoup plus faible), et l'obliquité plus réduite, produisent un plus grand effet. La ratio pôle / équateur calculé pour Mars sur le long terme est de 1.14. Cette réduction est presque entièrement due aux fortes variations d'obliquité (on l'a vu, les conditions d'approche ne diffèrent pas réellement). Pour l'ensemble des autres planètes, la prise en compte des variations séculaires ne fait pas varier le ratio pôle / équateur de plus de 2 %, ce que l'on tiendra pour négligeable. De la même façon, la distance Terre-Lune n'a pas d'incidence sur les variations en latitude du flux d'impacts lunaire, comme le montre la figure 6.11.

<figure><img src="image_52.png" /><figcaption>FIG. 6.11 – Variations latitudinales du flux d'impacts à la surface de la Lune, en fonction de la séparation Terre-Lune.</figcaption></figure>

La figure 6.12 montre la dépendance en latitude de l'angle moyen d'impact. Le comportement de l'angle d'impact est comparable à celui du flux. L'angle est en moyenne plus grand d'environ 2° à l'équateur qu'au pôles terrestres, et l'écart atteint 2.5° pour la Lune. A l'inverse, pour Mercure, Venus et Mars, il est plus grand au pôle d'environ 2°, 0.5° et 9°, respectivement (5° pour Mars sur le long-
terme). La vitesse d'impact est relativement constante et ne varie que de 0.2-0.5 km.s⁻¹ entre le pôle et l'équateur pour l'ensemble des planètes (figure 6.13).

<figure><img src="image_53.png" /><figcaption>FIG. 6.12 – Variations latitudinales de l’angle moyen d’impact à la surface des planètes telluriques.</figcaption></figure>

On note que chaque impact est susceptible de se produire avec une vitesse et un angle d'incidence qui diffèrent sensiblement de ses valeurs moyennes. Pour l'ensemble des planètes, l'écart-type est de l'ordre de 20° pour l'angle d'impact et d'environ 40 % de la valeur moyenne globale à la surface de la planète pour la vitesse d'impact. L'allure de la distribution de probabilité $p_λ(u, θ)$ dépend elle aussi de la latitude. Par exemple, sur la Lune, une plus grande proportion de collisions combinant faible incidence et grande vitesse auront lieu aux pôles. Nous verrons quelle influence cela peut avoir sur la formation des cratères d'impact aux chapitres suivants.

On attend de la Lune qu'elle présente également des variations longitudinales, en raison de la synchronicité de sa rotation. La figure 6.14 représente le flux d'impacts lunaire relatif, pour la distance Terre-Lune actuelle. Outre l'effet de latitude, un effet apex/antapex très significatif est observable. On remarque que, sur l'ensemble des croiseurs, aucun effet face visible / face cachée n'est observable (le ratio hémisphérique est de 1.001). Puisque l'effet apex/antapex se conjugue à l'effet de latitude, le flux est à son minimum entre le pôle et l'antapex, en l'occurrence vers (90°E, ±50°N). Le bombardement maximum s'observe à l'apex. La figure 6.15 représente la vitesse et l'angle d'impact moyens en fonction de la
<figure><img src="image_54.png" /></figure>

FIG. 6.13 – Variations latitudinales de la vitesse d'impact moyenne à la surface des planètes telluriques.

position. Très logiquement, la vitesse d'impact est plus importante à l'apex qu'à l'antapex. C'est aussi le cas de l'angle d'impact.

Il faut noter qu'au moment de la prime jeunesse du système solaire, lors de la phase d'accrétion des planètes, le bombardement était probablement généré en grande majorité par des projectiles en orbite planétocentrique. Nous ne tentons pas dans ce travail de quantifier les asymétries pouvant résulter d'un tel bombardement. Le lecteur pourra par exemple consulter le travail de Pinet [1985] appliqué au cas de la Lune. L'auteur conclut à une concentration du flux autour de l'équateur et de la face visible qui va s'amenuisant à mesure que la Lune s'éloigne de la Terre, jusqu'à présenter un excès polaire.

L'effet apex/antapex sera d'autant plus prononcé que la vitesse orbitale lu-
naire est grande, et donc que la distance Terre-Lune est petite. Sur la figure 6.16,
le flux d'impacts relatif est tracé en fonction de la distance angulaire à l'apex pour
trois distances Terre-Lune : 20, 40 et 60 rayons terrestres. Le ratio apex/antapex
est respectivement de 1.29, 1.35 et 1.54.

Intéressons nous maintenant au bombardement terrestre en fonction des co-
ordonnées orbitales (ce qui revient à rester dans le repère (OXYZ)). On procède
ainsi afin de pouvoir évaluer l'intensité du bombardement en fonction de l'heure
locale, ce qui est d'une grande utilité dans la mise en oeuvre de protocoles d'obser-
vation. Une longitude orbitale de 90°E correspond à 18h, 270°E à 6h. La partie
gauche de la figure 6.17 montre que, sur l'ensemble des projectiles , le bom-
<figure><img src="image_55.png" /><figcaption>FIG. 6.14 – Flux d'impacts à la surface de la Lune. Distance Terre-Lune actuelle.</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>FIG. 6.15 – Vitesses et angles d'impact moyens à la surface de la Lune. Distance Terre-Lune actuelle.</figcaption></figure>

bardement est beaucoup plus intense autour de midi et de minuit. Cependant, les méthodes de détection d'entrées atmosphériques sont en général sensibles à l'énergie cinétique des projectiles. La partie droite de la figure 6.17 pondère le flux d'impacts selon l'énergie cinétique. On constate alors que la majorité des impacts à haute énergie se produisent à l'aube. En restreignant nos résultats de telle façon que la vitesse d'impact est comprise entre 20 et 30 km/s, les observations radars (qui observent typiquement des particules de masses allant du micro au
<figure><img src="image_57.png" /><figcaption>FIG. 6.16 – Effet apex / antapex à la surface de la Lune, en fonction de la séparation Terre-Lune.</figcaption></figure>

milligramme) effectuées par Webster et al. [2004] et Jones et al. [2005] dans cette même gamme de vitesse sont relativement bien reproduites, comme le montre la figure 6.18. Cela laisse supposer que les particules qui bombardent la Terre sont des débris générés par de plus gros astéroïdes, et que leurs orbites sont par conséquent relativement similaires à celles utilisées dans notre modèle [Gallant et al., 2006].

Notons pour conclure que nous sommes capables, en utilisant les mêmes hypothèses, de reproduire à l'identique les résultats de Halliday [1964] et Halliday et Griffin [1982] avec la méthode analytique décrite dans ce chapitre. Dans Halliday [1964], le flux d'impacts relatif est calculé numériquement sur Terre pour des projectiles lancés parallèlement au plan orbital, et pour différentes vitesses d'impact. La partie gauche de la figure 6.19 reproduit exactement ces résultats. On rappelle cependant que l'hypothèse utilisée ($i_∞ = 0$) n'est pas réaliste. Dans Halliday et Griffin [1982], un jeu d'orbites supposé représentatif est extrait des observations
<figure><img src="image_58.png" /><figcaption>FIG. 6.17 – Flux d'impacts sur Terre en fonction de la latitude et de la longitude orbitales. A gauche Flux pour une taille donnée de projectiles. A droite Flux pondéré par l'énergie cinétique d'impact.</figcaption></figure>

d'entrées atmosphériques par le réseau canadien de caméras MORP. Le bombardement terrestre est alors calculé numériquement. La encore, nous retrouvons exactement les mêmes résultats, à savoir un ratio pôle / équateur de 0.85, et une variabilité du flux selon l'heure locale qui présente un maximum en soirée. Les orbites utilisées se basent néanmoins sur un très petit nombre d'observations.

## 6.6 Modéliser la bruit sismique lunaire

Pour conclure ce chapitre, nous nous proposons de comparer notre modèle aux impacts enregistrés durant les sept ans d'activité des stations sismiques lunaires Apollo 12, 14, 15 et 16. Celles-ci ont enregistré environ 2000 impacts pendant cette période. Le modèle nous permet également de modéliser le bruit sismique généré par le continuel bombardement de la surface lunaire par des micro-météorites, bruit qui n'a pas été enregistré par les stations sismiques, de résolutions insuffisantes. Une telle démarche est importante pour la réalisation d'une nouvelle génération de sismomètres lunaires, mais aussi pour la mise en place d'expériences in situ, qui requièrent un sol stable, à l'instar de mesures d'interférométrie ou d'ondes de gravité. Ce travail a été effectué sous la direction de Philippe Lognonné, avec la collaboration de Catherine Johnson.

L'amplitude du signal généré par un impact est considérée comme résultant d'une force ponctuelle (point force source). Le signal s enregistré à un temps t pour un impact ayant lieu à t = 0 s'exprime donc

$$s(t) = A \frac{m u}{L} e^{-\frac{\omega T}{2Q}}, \qquad (6.25)$$
<figure><img src="image_59.png" /><figcaption>FIG. 6.18 – Comparaison des observations radars d'entrées atmosphériques terrestres [Webster et al., 2004] avec notre simulation.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>FIG. 6.19 – Reproduction à l'identique des résultats de Halliday [1964] et Halliday et Griffin [1982].</figcaption></figure>

où *m* et *u* sont la masse et la vitesse de l'impacteur, *L* est la distance épicentrale qui sépare impact et station, *ω* = *π* est la fréquence angulaire correspondant au pic de fréquence de la station Apollo (0.5 Hz), et *T* le temps de propagation de
l'onde P (fonction de D). A et Q sont deux paramètres à déterminer, Q étant le facteur de qualité d'atténuation.

La figure 6.20 (à gauche) présente l'inversion par moindres carrés des paramètres A et Q, à partir des impacts artificiels générés par le programme Apollo (Saturn IV et Lunar Ascent Vehicle). Les amplitudes sont reproduites dans une gamme de ±50%, sur l'ensemble des distances épicentrales (de 67 à 1099 km), des valeurs de la quantité de mouvement $p = mu$ (de $0.4 \times 10^7$ à $3.6 \times 10^7$ kg.m.s⁻¹), et des angles d'impact (4 à 80 degrés). L'axe des ordonnées est donné en unités digitales Apollo (DU), telles que $1 \text{ DU} = 0.5 \times 10^{-10} \text{ m}$ en déplacement du sol et environ $0.5 \times 10^{-9} \text{ m.s}^{-2}$ en accélération. Une unité correspond à la limite de résolution des sismomètres installés actuellement sur la Lune.

<figure><img src="image_61.png" /><figcaption>FIG. 6.20 – A gauche. Ajustement des paramètres A et Q pour reproduire au mieux (courbe verte) les données. A droite. Exemple de signal généré par un impact synthétique (courbe noire), en comparaison des enregistrements Apollo.</figcaption></figure>

Un impact synthétique à une distance épicentrale quelconque est généré par une simple combinaison linéaire des deux signaux produits par les impacts artificiels aux distances épicentrales les plus proches, soit

$$s(t-T) = \frac{pE}{L} \left[ s_1(t-T_1) \frac{L_1}{p_1 E_1} + s_2(t-T_2) \frac{L_2}{p_2 E_2} \right], \quad (6.26)$$

avec $p, p_1$ et $p_2$ les quantités de mouvement de l'impact synthétique et des deux impacts artificiels, $T, T_1$ et $T_2$ les temps de propagation respectifs, et $E = e^{-\omega T/2Q}$, $E_1 = e^{-\omega T_1/2Q}$ et $E_2 = e^{-\omega T_2/2Q}$ les corrections d'amplitude par atténuation. La figure 6.20 (à droite) montre un exemple typique d'impact synthétique (en noir).

Il reste à présent à produire une « population » d'impacts à partir de notre modèle. Les incertitudes inhérentes à notre approche nous conduisent à négliger
ici les variations spatiales du bombardement. Le nombre d'impacts par unité de temps sur l'ensemble de la surface lunaire s'écrit en fonction de la masse du projectile comme

$$F(> m) = F(m > m_0) m^{-f}, \qquad (6.27)$$

où la masse s'exprime en grammes et $m_0 = 1$g. A partir du flux d'impacts lunaire absolu déterminé précédemment en fonction du diamètre ($\Phi_Q (> d)$), et en considérant que la densité des projectiles est de 2700 kg.m⁻³, on a $F(m > m_0) = 2.33 \times 10^6$ S.I. et $f = 0.97$. Cette loi de puissance est considéré comme valide jusqu'à $m \approx 1.5 \times 10^{12}$g ($d > 100$ m, figure 3.10).

Pendant un laps de temps donné $\Delta t$, le nombre total d'impacts pour lesquels $m > m_0$ est

$$N = \text{Int}[F(> m_0)\Delta t]. \qquad (6.28)$$

Sur une période de 7 ans (durée d'enregistrement des stations sismiques), on tire donc $N$ impacts aléatoires, dont la position est donnée par l'équation (6.1), dont la vitesse est échantillonnée selon la distribution des vitesses d'impact calculée pour la Lune (figure 5.8) et dont la masse est choisie comme

$$m = m_0 x^{-1/f}, \qquad (6.29)$$

$x$ étant un nombre uniformément distribué entre 0 et 1.

La figure 6.21 compare le modèle avec les enregistrements de la station Apollo 12. Le nombre absolu d'impacts synthétiques « enregistrés » dépend fortement de la valeur du facteur d'atténuation $Q$. Beaucoup des impacts simulés ont lieu sur la face cachée de la Lune (la moitié), tandis que peu d'événements ont été effectivement observés à des distances épicentrales de plus de 120° (rappelons que les stations sismiques sont toutes localisées sur la face visible). Le modèle reproduit bien le nombre absolu d'observations si l'on considère que $Q$ est divisé par 2 pour $L > 1500$ km. Une atténuation plus forte dans la Lune profonde a été suggérée par de nombreuses études (e.g. Nakamura [2004]). Contraindre le flux absolu de manière plus poussée nécessite de mieux connaître l'atténuation lunaire. Cependant, la pente de la distribution synthétique est en excellent accord avec les observations, et semble valider la distribution en taille de notre modèle. La courbe tiretée est un exemple de réalisation aléatoire sur 7 ans qui reproduit de manière particulièrement bonne les données, tandis que la courbe pleine est la distribution statistique sur une période de temps plus longue.

Le bruit sismique généré par le modèle est présenté en figure 6.22 pour une période de 10 jours. Le bruit de fond fluctue entre 0.02 et 0.1 DU, parsemé
<figure><img src="image_62.png" /></figure>

FIG. 6.21 – Comparaison du modèle avec les impacts enregistrés par la station sismique Apollo 12.

de plusieurs pics journaliers qui peuvent atteindre quelques dizaines d'unités. Le signal produit par les séismes profonds n'a pas été pris en compte, et est également susceptible de générer un bruit de fond significatif. Cependant, puisque ces événements ont une activité cyclique [Bulow et al, 2006], on peut s'attendre à ce que le signal présente une variation périodique dont le minimum correspond au bruit associé aux impacts, estimé ici.
<figure><img src="image_63.png" /><figcaption>FIG. 6.22 – Bruit sismique (modélisé) généré par le bombardement lunaire.</figcaption></figure>
# Chapitre 7

## De l'impact au cratère

Chaque impact produit un cratère, à condition que le projectile heurte une surface solide. Les propriétés de ce cratère, en particulier sa taille, dépendent de quantité de paramètres. Les propriétés du projectile, d'une part, qui regroupent sa taille, sa densité, sa vitesse et son angle d'incidence au moment de l'impact; celles de la cible, d'autre part, qui se composent de la densité et la cohésion du terrain, ainsi que de la gravité de surface. A cela s'ajoute la présence d'une éventuelle atmosphère, dont la traversée est susceptible d'avoir modifié les propriétés du projectile au moment où il atteint le sol (si elle ne lui a pas été fatale). Jusqu'à présent, nous n'avons étudié que le flux d'impacts. Une grandeur autrement intéressante, notamment en terme de chronologie, est, on l'aura compris, le taux de formation des cratères, c'est-à-dire le nombre de cratères plus grands qu'une certaine taille se formant par unité de temps et de surface. Pour relier l'un à l'autre, il nous faut connaître précisément les lois d'échelles empiriques qui reproduisent la formation d'un cratère d'impact.

Ces lois d'échelles ont fait l'objet de nombreuses études, mais certains paramètres (notamment l'influence de l'angle d'incidence) restent mal contraints. La plus complète synthèse sur le sujet reste à ce jour celle de Melosh [1989]. La démarche clef consiste en l'étude adimensionnée d'expériences de tirs supersoniques en laboratoire ou d'explosions nucléaires, propre à isoler les grandeurs indépendantes de l'échelle du phénomène et les coefficients qui y sont associés [Schmidt et Housen, 1987; Holsapple et Schmidt, 1987; Holsapple, 1993; Gault, 1974]. L'applicabilité à grande échelle des relations ainsi obtenues reste délicate, puisque certains processus physiques associés aux grands impacts ne sont pas reproductibles. C'est le cas notamment de la vaporisation produite par des pressions et des températures extrêmes, ou encore de l'influence de la gravité dans les stades terminaux de la formation du cratère. L'usage de simulations numériques [Pierazzo et Melosh, 2000] et de l'observation de phénomènes d'échelle planétaire (par exemple l'impact de la comète Shoemaker-Levy avec Jupiter) procurent des
contraintes supplémentaires. La plus récente mise à jour des lois d'échelle est due
à Holsapple et Housen [2007], où les auteurs se consacrent à l'étude des résultats
de la mission *Deep Impact*.

Nous donnons dans un premier temps les relations qui permettent de calculer
la taille du cratère. Nous appliquons ensuite ces relations à la construction de
distributions taille / fréquence des cratères, à partir de notre modèle de bom-
bardement, en ajustant les paramètres mal contraints à l'aide des distributions
observées, notamment sur la Lune. Nous faisons alors état de nos résultats concer-
nant le taux relatif de formation des cratères entre les planètes. Les variations
spatiales de ce taux à la surface des planètes sont laissées au chapitre suivant.

7.1 Formation du cratère

Comme le mentionnent Holsapple et Housen [2007], le couplage de l'énergie
et de la quantité de mouvement de l'impact avec le matériau de la cible s'étend
sur une zone dont la taille est de l'ordre de celle du projectile. Dans la mesure où
le cratère formé est bien plus grand que le projectile, le processus peut-être traité
comme résultant d'une source ponctuelle. Pour un impact vertical et un type
de matériau cible donné, les processus qui caractérisent la formation du cratère
(vitesse d'excavation notamment) dépendent d'une seule mesure, le paramètre
de couplage : $du^\mu \rho^\nu$, où d est le diamètre du projectile, u sa vitesse d'impact,
ρ sa densité et μ et ν deux exposants dont la valeur dépend de la nature de
la cible. L'exposant μ (1/3 < μ < 2/3) décrit dans quelle mesure c'est, soit
l'énergie cinétique, soit la quantité de mouvement qui gouverne le processus.
Pour μ = 2/3, il n'y a pas de dissipation de l'énergie et c'est cette dernière qui
contrôle l'excavation; à mesure que la porosité de la cible augmente, l'énergie se
dissipe plus facilement, jusqu'à laisser la main-mise à la quantité de mouvement
quand μ = 1/3.

7.1.1 Cratère transitoire

On rappelle que le cratère transitoire est défini comme le cratère qui s'est
formé à la fin du processus d'excavation et de relaxation, avant d'éventuelles
modifications par effondrement gravitaire qui ne sont pas reproductibles à l'échelle
du laboratoire. On note $D_⊥$ le cratère transitoire formé par un impact purement
vertical. Dans le cas le plus général, Holsapple et Housen [2007] donnent la relation

$$ \frac{D_{\perp}}{d} = K \left[ \frac{g}{u^2} \left( \frac{\rho_0}{\rho} \right)^{\frac{2\nu}{\mu}} + \left( \frac{\bar{Y}}{\rho_0 u^2} \right)^{\frac{2+\mu}{2}} \left( \frac{\rho_0}{\rho} \right)^{\frac{\nu(2+\mu)}{\mu}} \right]^{-\frac{\mu}{2+\mu}}, \quad (7.1) $$

où g est la gravité de surface, ρ la masse volumique du projectile, ρ₀ et ¯Y respectivement la masse volumique et la cohésion moyenne du sol. La constante K s'ajoute à ν et μ pour caractériser la réponse de ce sol à l'impact. Notons que ¯Y représente la moyenne de la cohésion en régime de tension et de la cohésion en régime de compression.

On se situe en régime de contrainte lorsque l'échelle du phénomène est suffisamment réduite pour que la gravité n'intervienne pas. On pose alors g = 0 dans l'équation précédente, et c'est la cohésion du matériau cible qui est déterminante. A l'inverse, les grands cratères se forment en régime gravitaire, et l'on pose ¯Y = 0.

En régime gravitaire, on peut donc écrire

$$ \frac{D_{\perp}}{d} = K \left( \frac{u^2}{gd} \right)^{\nu_1} \left( \frac{\rho}{\rho_0} \right)^{\nu_2}. \quad (7.2) $$

Pour des roches non poreuses, les constantes s'expriment K = 1.17, ν₁ = 0.22 et ν₂ = 0.31, tandis que pour un sable sec on prendra K = 1.03, ν₁ = 0.17 et ν₂ = 0.332 [Schmidt et Housen, 1987; Holsapple et Housen, 2007].

On sait que la valeur de l'angle d'impact θ influe sur la taille du cratère, mais la relation fonctionnelle n'est pas clairement établie. Beaucoup d'auteurs considèrent que seule la composante verticale de la vitesse d'impact contribue à former le cratère [e.g. Pierazzo et al., 1997]. En revanche, des mesures de laboratoire effectuées par Gault et Wedekind [1978] suggèrent que le volume de matière éjecté est proportionnel à sin θ. Si le ratio profondeur / diamètre du cratère ne change pas, le diamètre sera alors proportionnel à (sin θ)¹/³ [voir Melosh, 1989]. Beaucoup d'études considèrent simplement que tous les impacts sont verticaux [e.g. Shoemaker, 1983]. On a vu (figure 6.12) que la distribution de l'angle d'impact varie en fonction de la position à la surface des corps bombardés. Le taux de formation des cratères, qui s'exprime pour une taille donnée, dépendra donc de la manière dont l'angle d'incidence est pris en compte. Ecrivons que le diamètre transitoire est donné en fonction de θ comme

$$ D_T = D_{\perp}(\sin \theta)^{\nu_3}, \quad (7.3) $$

et deux possibilités seront considérées. Dans un premier cas, on considère que l'angle d'impact n'intervient pas et l'on pose ν₃ = 0. Comme alternative, on
suppose que seule la composante verticale de la vitesse intervient et on pose $\nu_3 = 2\nu_1$. Chacune des relations fonctionnelles proposées plus haut donneront des résultats encadrés par ces deux cas limites.

### 7.1.2 Cratères simples

Les cratères simples, typiquement de taille kilométrique, résultent de certaines modifications du cratère transitoire et présentent une taille finale plus importante. D'après Melosh [1989], le glissement vers le bas des bords du cratère augmente le diamètre d'un facteur 1.25. Un autre facteur 1.25 tient compte du fait que la bordure du cratère se forme à une élévation plus grande que la surface initiale. Puisque l'on mesure le diamètre d'un cratère bord à bord, celui-ci est défini plus grand que s'il était mesuré au niveau de la topographie alentour.

Le diamètre d'un cratère simple est donc

$$D_s = 1.56D_T. \qquad (7.4)$$

### 7.1.3 Cratères complexes

Au-delà d'une certaine taille, la gravité provoque un effondrement encore plus important des parois du cratère. On note $D_*$ le diamètre de transition au-delà duquel le phénomène se produit. En première (et, jusqu'à présent, meilleure) approximation, $D_*$ est inversement proportionnel à la gravité de surface [Pike, 1980]. Sur Terre, $D_*$ est relativement bien contraint autour de 3 km. Les valeurs correspondantes sont 18 km pour la Lune, 8 km pour Mercure et Mars, et également 3 km pour Venus. La relation qui lie cratères simples et cratères complexes n'est pas non plus très bien définie, mais on utilise communément celle donnée par Croft [1985] :

$$D = \begin{cases} D_s & \text{pour } D_s < D_* \\ D_*^{-\nu_4} D_s^{1+\nu_4} & \text{pour } D_s \ge D_*, \end{cases} \qquad (7.5)$$

avec $\nu_4 = 0.18$.

### 7.1.4 Cratères en régime de contrainte

Dans le cas des petits cratères, c'est, on l'a dit, la cohésion de sol qui gouverne le processus de formation. On constate en posant $g = 0$ dans l'équation (7.1) que le diamètre de transition du cratère (qui est aussi son diamètre final ici)
est directement proportionnel à celui du projectile. Plutôt que de chercher à déterminer une valeur représentative de la cohésion des sols lunaires ou martiens, Ivanov [2001] construit la transition cratères de contrainte / cratères simples en créant une fonction continue et lisse donnée par les équations (7.2) et (7.4) quand $D \gg D^\circ$, et proportionnelle à $d$ quand $D \ll D^\circ$, où $D^\circ$ est le diamètre de transition observé entre régime de contrainte et régime de gravité.

C'est-à-dire,

$$D (D + D^\circ)^{\frac{\nu_1}{1-\nu_1}} = d \left[ 1.56K \left( \frac{(u \sin \theta)^2}{g} \right)^{\nu_1} \left( \frac{\rho}{\rho_0} \right)^{\nu_2} \right]^{\frac{1}{1-\nu_1}} . \quad (7.6)$$

Là encore, on peut considérer que $D^\circ$ est inversement proportionnel à la gravité de surface. Sur Terre, on observe que $D^\circ \simeq 50\text{m}$. On a donc $D^\circ \simeq 300\text{ m}$ sur la Lune et sur Mercure, et $D^\circ \simeq 150\text{m}$ sur Mars.

# 7.2 Distributions taille-fréquence

Nous décrivons ici les outils nécessaires au calcul des distributions taille-fréquence des cratères à partir du flux et des conditions d'impact déterminées au chapitre 6. On ne considère ici que la formation des cratères d'impact à l'échelle de la planète entière. Les variations spatiales sont laissées au chapitre suivant.

## 7.2.1 Formalisme

On rappelle que le flux de projectiles d'un diamètre supérieur à $d$ qui entrent en collision avec la planète à une certaine époque est donné par

$$\Phi(> d, t) = \eta(t) \Phi(> d) ,$$

avec

$$\Phi(> d) = \Phi(d > 1) S(> d) ,$$

en supposant que la distribution orbitale des astéroïdes est indépendante du diamètre considéré et du temps. De la même manière, on écrit le taux de formation des cratères, c'est-à-dire le nombre de cratères d'un diamètre supérieur à $D$ qui se forment par unité de temps et de surface (en moyenne sur toute la surface de la planète), comme

$$\Phi_c(> D, t) = \eta(t) \Phi_c(> D) . \quad (7.7)$$
En utilisant la relation $d(D, u, \theta)$, qui lie le diamètre du projectile au diamètre du cratère et est donnée par les équations (7.2)–(7.5), on écrit

$$ \Phi_c(> D) = \Phi(d > 1) \bar{S}(> d(D)), \qquad (7.8) $$

avec

$$ \bar{S}(> d(D)) = \frac{1}{2\pi^2} \int_{\lambda_0, \varphi_0} \int_{U_\infty} S(> d(D, u(U_\infty), \theta(\lambda_0, \varphi_0, U_\infty))) \quad (7.9) $$

$$ \times \phi(\lambda_0, \varphi_0, U_\infty) p(U_\infty) dU_\infty d\lambda_0 d\varphi_0. \qquad (7.10) $$

Une autre notation, plus claire, pourrait être

$$ \bar{S}(> d(D)) = \iint S(> d(D, u, \theta)) p(u, \theta) du d\theta, \qquad (7.11) $$

en supposant que l'on a construit la distribution de probabilité des conditions d'impact $p(u, \theta)$ telle que $\iint p(u, \theta) du d\theta = 1$. Une approximation qui peut s'avérer utile car facile d'utilisation est

$$ \bar{S}(> d(D)) \approx S(> d(D, \bar{u}, \bar{\theta} = \pi/4)) . \qquad (7.12) $$

La distribution taille / fréquence de la densité moyenne des cratères qui re-couvrent une planète dont la surface possède un âge $A$ est simplement

$$ N(> D, A) = \int_{0}^{t} \Phi_{c}(> D, t) dt = \Phi_{c}(> D) \int_{0}^{A} \eta(t) dt, \qquad (7.13) $$

quantité directement proportionnelle à $A$ si l'on suppose que le flux d'impact est resté constant durant la période considéré ($\eta(t) = 1$). La distribution taille / fréquence incrémentale est donnée par

$$ N(D, A) = \frac{dN(> D, A)}{dD}, \qquad (7.14) $$

et la distribution relative est

$$ r(D, A) = D^3 N(D, A). \qquad (7.15) $$

Le taux de formation des cratères sur une planète donnée comparativement à la Lune s'exprime à travers la grandeur $R_c$ telle que

$$ R_c(> D) = \frac{\Phi_c(> D)}{\Phi_c^c(> D)} = R \frac{\bar{S}(> d(D))}{\bar{S}(> d(D))_c}, \qquad (7.16) $$

où l'on rappelle que $R$ est le flux d'impacts relatif à la Lune.

Même si $\Phi_c(> D)$ s'obtient numériquement sans difficulté pour n'importe quelle distribution supposée de la taille des projectiles (c'est d'ailleurs ce qui est
fait dans la section suivante à partir de la distribution établie au chapitre 3), il est utile de supposer que celle-ci suit une simple loi de puissance afin de dégager des équations précédentes un comportement lisible. Considérons donc un moment que $S(> d) = d^{-n}$, avec $n$ la pente de la distribution, soit

$$ \Phi(> d) = \Phi(d > 1) d^{-n}. \qquad (7.17) $$

On fait également l'hypothèse que tous les cratères sont simples. En insérant les équations (7.4)–(7.2) dans l'équation (7.8), on a

$$ \Phi_c(> D) = \Phi_c(D > 1)D^{-s}. \qquad (7.18) $$

où le taux de formation de cratères d'un diamètre supérieur au kilomètre est

$$ \Phi_c(D > 1) = \Phi(d > 1) \left[ 1.56 K g^{-\nu_1} \left( \frac{\rho}{\rho_0} \right)^{\nu_2} \right]^s \langle \left( u^{2\nu_1} \sin^{\nu_3} \theta \right)^s \rangle, \qquad (7.19) $$

avec $\langle \dots \rangle$ la moyenne sur la distribution des conditions d'impact, et où la pente de la distribution taille / fréquence des cratères est

$$ s = \frac{n}{1 - \nu_1}. \qquad (7.20) $$

Si on suppose que la densité de surface est la même, la production relative de cratères par rapport à la Lune est donc,

$$ R_c = R \left( \frac{g\zeta}{g} \right)^{\nu_1 s} \frac{\langle (u^{2\nu_1} \sin^{\nu_3} \theta)^s \rangle}{\langle (u^{2\nu_1} \sin^{\nu_3} \theta)^s \rangle \zeta}, \qquad (7.21) $$

et l'on note que $R_c$ ne dépend pas de $D$ dans ce cas. Une approximation raisonnable est donnée en prenant directement les valeurs moyennes de la vitesse d'impact et de l'angle d'impact, plutôt que $\langle (u^{2\nu_1} \sin^{\nu_3} \theta)^s \rangle$ (équation (7.12)). L'angle d'impact moyen étant toujours de $\pi/4$, on a

$$ R_c(> D) = R \left[ \frac{g\zeta}{g} \left( \frac{\bar{u}}{\bar{u}\zeta} \right) \right]^{\nu_1 s}. \qquad (7.22) $$

Pour $s \approx 2-3$, la pente généralement observée sur l'ensemble des planètes, on se fait une idée rapide du taux de production relatif des cratères. On constate que le ratio pour un diamètre de cratère donné, peut, selon les différences de gravité de surface et de vitesse d'impact, être accru ou réduit par rapport au ratio pour une taille de projectile donnée : un objet d'une certaine taille ne produira pas le même cratère selon la planète sur laquelle il échoue. Sur Mars, par exemple, il sera plus petit que sur la Lune, car la vitesse d'impact est en moyenne deux fois moindre. Ainsi, même si $R$ est grand, nous verrons que $R_c$ est proche de l'unité pour cette planète.
### 7.2.2 Reconcilier l'observation des cratères et celle des boîtes

Il s'agit maintenant de vérifier que la distribution en taille des astéroïdes, telle qu'elle a été déterminée à partir des observations spatiales, est à même de générer des distributions de cratères qui reproduisent les observables. On choisit comme référence la fonction polynomiale établie par Neukum et al. [2001a] pour la distribution en taille des cratères lunaires (équation (1.5)). Celle-ci, d'après les auteurs, est restée inchangée durant toute l'histoire lunaire, et seul le nombre absolu de cratères varie avec le temps. Qui plus est, ils s'appuient sur les âges des échantillons lunaires pour conclure (en comparant ceux-ci avec les densité de cratères qui recouvrent les unités géologiques associées) que le bombardement est resté relativement constant pendant les trois derniers milliards d'années. On choisit donc d'effectuer la comparaison pour une surface vieille de 3 milliards d'années, en calculant $N_Q(D > 1,3 \text{ Ga})$. Pour cela, on fait usage des équations (7.13), (7.10) et (7.8) appliquées à la Lune, en posant $A = 3 \text{ Ga}$ et $\eta(t) = 1$. On considère également que $d(D, u, \theta)$ est donné pour $\nu_3 = 2\nu_1$ (composante verticale de la vitesse d'impact). Enfin, les masses volumiques moyennes des projectiles et de la cible sont prises comme $\rho = 2050 \text{ kg.m}^{-3}$ [Stuart et Binzel, 2004] et $\rho_0 = 2800 \text{ kg.m}^{-3}$.

La figure 7.1 compare la relation établie par Neukum et al. [2001a] avec notre synthétique. On constate qu'il a été nécessaire de faire intervenir deux régimes distincts de formation des cratères : le régime poreux pour $D_T < 5 \text{ km}$, le régime non poreux pour $D_T > 15 \text{ km}$ (équation (7.2) avec les valeurs de $K$, $\nu_1$ et $\nu_2$ appropriées). Dans l'intervalle, une simple transition linéaire a été utilisée. Sans cette transition, superposer les deux distributions n'est pas possible (courbes tiretées). Cela mis à part, l'accord est excellent. En général, la majorité des auteurs considèrent que les cratères se forment en régime non poreux sur les planètes telluriques. Il semble, comme Ivanov [2006] l'a mis en évidence, que le seul moyen de réconcilier l'observation des astéroïdes et l'observation des cratères soit de considérer que le processus de formation des cratères dans la proche subsurface, c'est-à-dire dans le régolithe, est bien représenté par un régime poreux.

Cela peut se comprendre : le régolithe, agrégat de morceaux de roches de tailles variées, peut présenter un comportement similaire au sable des laboratoires à l'échelle planétaire. Les travaux de Lognonné et al. [2003], par exemple, mettent en évidence par inversion des données sismiques lunaires que le premier kilomètre
<figure><img src="image_64.png" /><figcaption>FIG. 7.1 – SFD synthétique en comparaison de la courbe de Neukum, en incluant une transition linéaire du régime d'impact.</figcaption></figure>

de la subsurface est marqué par une zone de faible vitesse sismique (~1 km.s⁻¹ pour les ondes P). Puisque la profondeur d'un cratère est de l'ordre de D/5 [Melosh, 1989], les cratères se forment entièrement dans cette zone jusque $D_T \simeq 5$km, et l'on retrouve notre transition régime poreux / régime non poreux. Reste à savoir si l'on peut raisonnablement supposer que ce comportement est le même sur les autres planètes. La question reste ouverte, et peut remettre en cause l'exportation de la courbe de Neukum sur les autres planètes. En l'absence de réponse, on considérera que la transition est identique sur les autres objets de cette étude.

On remarque que l'hypothèse d'un flux de projectiles constant sur les trois derniers milliards d'années, associée à une loi d'échelle qui se base sur la composante verticale de la vitesse d'impact, reproduit de manière excellente le nombre absolu de cratères formés sur cette période de temps. Bien sûr, il faut garder à l'esprit que de grandes incertitudes existent sur l'ensemble des étapes qui ont
mené à l'obtention de cette courbe, et que celles-ci sont susceptible de se contre-balancer favorablement pour aboutir à ce résultat (étonnamment ?) en accord. On note en particulier que notre modèle n'inclut pas la contribution des cratères secondaires, potentiellement importants pour $D \sim 1$ km.

La relation $d(D, u, \theta)$ qui a permis produire la figure 7.1 s'écrit donc

$$ 
\left\{ 
\begin{array}{ll} 
D_s = D & \text{si } D < D_*, \\ 
D_s = D_*^{0.15} D_*^{0.85} & \text{si } D \ge D_*, \\ 
D_T = D_s/1.56, \\ 
d_1 = \left(0.97 D_1 (u \sin \theta)^{-0.34} g^{0.17} \left(\frac{\rho_0}{\rho}\right)^{0.33}\right)^{1.2}, & \\ 
d_2 = \left(0.85 D_2 (u \sin \theta)^{-0.44} g^{0.22} \left(\frac{\rho_0}{\rho}\right)^{0.31}\right)^{1.28}, & \\ 
d = \left(0.97 D_T (u \sin \theta)^{-0.34} g^{0.17} \left(\frac{\rho_0}{\rho}\right)^{0.33}\right)^{1.2} & \text{si } D_T < D_1, \\ 
d = d_1 + \frac{d_2 - d_1}{D_2 - D_1} (D_T - D_1) & \text{si } D_1 \ge D_T \le D_2, \\ 
d = \left(0.85 D_T (u \sin \theta)^{-0.44} g^{0.22} \left(\frac{\rho_0}{\rho}\right)^{0.31}\right)^{1.28} & \text{si } D_T > D_2. 
\end{array} 
\right. 
\quad (7.23) 
$$

avec $D_1 = 5$ km et $D_2 = 15$ km.

Dans ces conditions, la grandeur $\bar{S}(> d(D))$ se calcule pour la Lune comme

$$ 
\begin{aligned} 
\bar{S}(> d(D))_{\mathcal{C}} &= 442.5 D^{-3.5} && \text{pour } 0.1 \le D < 2 \text{ km,} \\ 
&= 156.3 D^{-2} && \text{pour } 2 \le D < 8 \text{ km,} \\ 
&= 19.5 D^{-1} && \text{pour } 8 \le D < 26 \text{ km,} \\ 
&= 974.4 D^{-2.2} && \text{pour } 26 \le D < 58 \text{ km,} \\ 
&= 11138.1 D^{-2.8} && \text{pour } 58 \le D < 300 \text{ km,} \\ 
&= 3.8 D^{-1.4} && \text{pour } 300 \le D < 1000 \text{ km.} 
\end{aligned} 
\quad (7.24) 
$$

et le taux de caractérisation lunaire actuel (par unité de temps et de surface) est donné par l'équation (7.8). Puisque le flux de projectiles a été calculé comme $\Phi_c^c(d > 1) = 2.1 \times 10^{-15} \text{km}^{-2}\text{an}^{-1}$, on a, par exemple pour $D > 1$ km, $\Phi_c^c(D > 1) = 9.3 \times 10^{-13} \text{km}^{-2}\text{an}^{-1}$.

### 7.2.3 Comparaison interplanétaire

En considérant que les mêmes hypothèses s'appliquent sur l'ensemble des planètes telluriques (transition régime poreux / non poreux entre 5 et 15 km; densité du sol de 2800 kg.m$^{-3}$; même distribution taille / fréquence des astéroides), nous calculons ici le taux de formation des cratères relatif à la Lune, $R_c(> D)$, qui dépend, outre du bombardement relatif $R$, des conditions d'impact différentes et de la gravité de surface. Les valeurs obtenues à partir de l'équation (7.16) sont
<table><caption>TAB. 7.1 – Taux de cratérisation relatif à la Lune, pour $D > 1$km et $D > 30$km, comparés aux bombardement relatifs.</caption><thead><tr><th>Planète</th><th>Mercure</th><th>Vénus</th><th>Terre</th><th>Lune</th><th>Mars<br>(moyenne séculaire)</th></tr></thead><tbody><tr><td>$R$</td><td>1.81</td><td>1.73</td><td>1.57</td><td>1</td><td>2.74</td></tr><tr><td>$R_c$ ($D > 1$ km)</td><td>2.8</td><td>0.86</td><td>0.59</td><td>1</td><td>0.82</td></tr><tr><td>$R_c$ ($D > 30$ km)</td><td>3.5</td><td>1.7</td><td>1.2</td><td>1</td><td>1.2</td></tr></tbody></table>

regroupées en Table (7.1), pour $D$ supérieur à 1 et 30 km. La fonction $R_c(> D)$ est tracée sur la figure 7.2. Il faut noter que la valeur de $R_c$ est assez fortement sensible aux différents paramètres de cette étude, notamment la pente de la distribution d'astéroïdes et la loi d'échelle utilisée. Les changements de l'un et de l'autre en fonction de la taille produisent des fluctuations conséquentes (de l'ordre du facteur 2). Les valeurs de $R_c$ données dans la publication présentée en annexe D se basent sur une reproduction plus grossière de la distribution taille-fréquence de Neukum et al. [2001a], et les valeurs listées pour $D > 1$ km correspondent ici à $D > 2$ km. On constate sur la figure 7.2 une transition marquée entre 1 et 2 km (entre autres) qui provient du changement de pente de la distribution d'astéroïdes à $d = 100$ m ; dans la mesure où $R_c(D > 1$km) est un outil fondamental de l'exportation extra-lunaire de la méthode de datation par comptage de cratères, de plus amples travaux sont requis pour s'assurer du bon emplacement de cette transition.

L'atmosphère n'a pas été prise en compte dans ces estimations. Cela n'a évidemment aucune incidence pour la Lune et Mercure. L'estimation de l'influence de l'atmosphère passée de Mars (très probablement beaucoup plus épaisse et dense qu'aujourd'hui) sur la formation des cratères prévoit que seuls les cratères de taille inférieure à 250 m voient leur nombre significativement réduit par le passage atmosphérique [Chappelow et Sharpton, 2005]. L'effet est beaucoup plus important pour Venus. Les travaux de McKinnon et al. [1997] et de Korycansky et Zahnle [2005] laissent entendre que l'influence atmosphérique est sensible jusqu'à des diamètres de plusieurs dizaines de kilomètres.

On constate que $R_c$ est en général réduit par rapport à $R$. Seule Mercure, à la surface de laquelle la gravité est très faible et les vitesses d'impact très grandes, montre un taux de formation des cratères nettement plus élevé que les autres
<figure><img src="image_65.png" /></figure>

FIG. 7.2 – Taux de formation des cratères relatif à la Lune, pour des cratères plus grands que $D$, en fonction de $D$.

planètes. En d'autre termes, on peut dire qu'à âge égal, une surface mercurienne aura accumulé beaucoup plus de cratères d'une certaine taille que n'importe quelle autre surface planétaire. La figure 7.3 trace les distributions taille / fréquence calculées pour une période de 3 milliards d'années sur la Lune, Mercure, Venus et Mars (moyenne séculaire). On rappelle qu'il s'agit de distributions moyennes calculées sur l'ensemble de la planète. Les variations spatiales sont abordées au chapitre suivant.
<figure><img src="image_66.png" /><figcaption>FIG. 7.3 – Comparaison des SFD planétaires globales.</figcaption></figure>
# Chapitre 8

## Variations spatiales du taux de cratérisation

Le chapitre 6 a vu la détermination des conditions d'impact en fonction de la position, tandis que le chapitre 7 a donné les outils nécessaires au passage de l'impact au cratère. Ici, nous en déduisons les variations spatiales du taux de formation des cratères à la surface de la Lune, Mercure, Vénus, la Terre et Mars. En particulier, nous nous intéressons au taux de formation des cratères d'un diamètre supérieur à 1 kilomètre, puisque c'est la mesure principalement utilisée dans la détermination d'un âge par comptage des cratères d'impact.

### 8.1 Formalisme

Rappelons pour commencer que la variation relative du flux d'impacts en un point donné, par rapport à la moyenne sur l'ensemble de la planète, a été baptisée $\phi(\lambda, \varphi)$. De manière similaire, on notera $\phi_c(> D, \lambda, \varphi)$ le taux relatif de formation des cratères. Si les variations spatiales du flux d'impacts ne dépendent pas de la taille $d$ des objets considérés, le taux de cratérisation dépend a priori de la valeur de $D$ considérée, puisque les lois qui relient $d$ à $D$ dépendent des conditions d'impact de manière non linéaire. Il s'exprime de manière similaire à l'équation (7.16) comme

$$ \phi_c(> D, \lambda_0, \varphi_0) = \phi(\lambda_0, \varphi_0) \frac{\bar{S}(> d(D))_{\lambda_0, \varphi_0}}{\bar{S}(> d(D))} \quad (8.1) $$

avec

$$ \bar{S}(> d(D))_{\lambda_0, \varphi_0} = \int_{U_{\infty}} S(> d(D, u(U_{\infty}), \theta(\lambda_0, \varphi_0, U_{\infty}))) \quad (8.2) $$

$$ \times \phi(\lambda_0, \varphi_0, U_{\infty}) p(U_{\infty}) dU_{\infty} \quad (8.3) $$

ou encore

$$ \bar{S}(> d(D))_{\lambda_0, \varphi_0} = \iint S(> d(D; u, \theta)) p_{\lambda_0, \varphi_0}(u, \theta) du d\theta, \quad (8.4) $$
avec $p_{\lambda_0, \varphi_0}(u, \theta)$ la distribution de probabilité des conditions d'impact à la position $(\lambda_0, \varphi_0)$. $\bar{S}(> d(D))$, la moyenne planétaire globale, est donnée par l'équation (7.10) (ou 7.11). Si la planète possède une obliquité significative, le passage du repère orbital aux coordonnées géographiques $(\lambda, \varphi)$ se fait par rotation de l'axe de spin, comme détaillé au chapitre 6 (équation (6.20), en substituant $\Phi_c$ à $\Phi$).

La distribution taille / fréquence des cratères qui recouvrent une surface d'âge $A$ et de position $(\lambda, \varphi)$ est simplement donnée par

$$N(> D, A, \lambda, \varphi) = \phi_c(> D, \lambda, \varphi) N(> D, A), \quad (8.5)$$

où $N(> D, A)$, la distribution taille / fréquence moyenne sur la totalité de la planète, est donnée par l'équation (7.13).

Pour extraire un sens physique clair des équations précédentes, supposons un instant, comme au chapitre précédent, que la distribution de taille des projectiles suit une simple loi de puissance telle que $S(> d) = d^{-n}$. Dans ce cas, si on considère que tous les cratères formés sont simples, et qu'il ne se forment que dans un seul type de régime (poreux ou non poreux) on a

$$\phi_c(\lambda, \varphi) = \phi(\lambda, \varphi) \frac{\langle u^{2\nu_1 s} \sin^{\nu_3 s} \theta \rangle_{\lambda, \varphi}}{\langle u^{2\nu_1 s} \sin^{\nu_3 s} \theta \rangle}, \quad (8.6)$$

indépendant de $D$ dans ce cas, mais pas dans le cas général. Notons qu'une approximation raisonnable de $\phi_c$ peut être obtenue en considérant que $s \simeq 2-3$, ce qui reproduit approximativement la pente des distributions observables sur l'ensemble des planètes telluriques, et en insérant directement dans l'équation précédente les valeurs moyennes de la vitesse et de l'angle d'impact à une position donnée, telles qu'elles ont été déterminées au chapitre 6, plutôt que de calculer $\langle u^{2\nu_1 s} \sin^{\nu_3 s} \theta \rangle_{\lambda, \varphi}$ sur l'ensemble de la distribution des conditions d'impact à cette position.

Sur la Lune, par exemple, nous avons calculé que non seulement le flux d'impacts, mais aussi la vitesse et l'angle d'impact, étaient en moyenne plus faibles au pôle qu'à l'équateur. Un projectile d'une taille donnée produira donc, en moyenne, un cratère plus petit s'il tombe au pôle plutôt qu'à l'équateur. Cela accentue les variations du taux de formation des cratères par rapport à celles du flux d'impacts, comme nous allons le voir.
## 8.2 Résultats

La figure 8.1 représente la variation latitudinale du taux de formation des cratères de diamètre supérieur à 1 kilomètre ($\phi_c(D > 1, \lambda)$) sur les planètes telluriques et la Lune. On rappelle que l'on considère que la vitesse d'impact contribue à la taille du cratère uniquement par l'intermédiaire de sa composante verticale ($\nu_3 = 2\nu_1$). Puisque l'angle d'impact est fonction de la latitude, et se comporte de manière similaire au flux d'impacts (un flux plus élevé est associé à un angle d'impact moyen plus grand), les variations en latitude sont accrues lorsque l'on parle de cratères d'une certaine taille plutôt que d'impacts (un fort angle d'impact crée un cratère plus grand). Les ratios pôle / équateur sont résumés en Table 8.1, et comparés au cas où la relation $d(D)$ ne prend pas en compte l'angle d'impact ($\nu_3 = 0$). Dans cette hypothèse, seules les variations latitudinales de la norme du vecteur vitesse d'impact sont susceptibles d'accroître ou de réduire les variations par rapport à celles du flux d'impact. On constate que cet effet est plutôt modéré.

<figure><img src="image_67.png" /><figcaption>FIG. 8.1 – Taux de formation des cratères de diamètre supérieur à 1 km à la surface des planètes telluriques, en fonction de la latitude.</figcaption></figure>
<table><caption>TAB. 8.1 – Quotients pôle / équateur du flux d'impacts et du taux de formation des cratères pour $D > 1$ km. (1) Loi d'échelle basée sur la composante verticale de la vitesse d'impact. (2) Loi d'échelle pour des impacts verticaux.</caption><thead><tr><th>Ratio</th><th>Mercure</th><th>Venus</th><th>Terre</th><th>Lune</th><th>Mars</th><th>Mars<br>(long terme)</th></tr></thead><tbody><tr><td>Flux d'impacts</td><td>1.05</td><td>1.00</td><td>0.96</td><td>0.90</td><td>1.26</td><td>1.14</td></tr><tr><td>Production de cratères (1)</td><td>1.08</td><td>1.00</td><td>0.90</td><td>0.83</td><td>1.56</td><td>1.27</td></tr><tr><td>Production de cratères (2)</td><td>1.04</td><td>0.99</td><td>0.95</td><td>0.87</td><td>1.30</td><td>1.16</td></tr></tbody></table>

Même si $\phi_c(> D, \lambda)$ est en toute rigueur fonction du diamètre $D$ considéré, nous calculons qu’un écart de moins de 5% existe entre les ratios présentés en Table 8.1 pour $D > 1$ km et les ratios correspondant à n’importe quel autre diamètre. Ceci est illustré en figure 8.3, qui montre les distributions taille / fréquence des cratères martiens et lunaires, au pôle et à l'équateur. On constate qu'en dépit des distributions de probabilité différentes en terme de conditions d'impact, seul un décalage vertical est évident entre les courbes produites sur un même corps. La figure 8.2 représente les variations en fonction de $D$ du ratio pôle/équateur de cratérisation.

<figure><img src="image_68.png" /><figcaption>FIG. 8.2 – Variations du ratio pôle/équateur de cratérisation selon $D$.</figcaption></figure>

Il faut cependant noter que le fait que les conditions d'impact ne présentent pas la même distribution de probabilité selon la latitude influe sur la proportion de cratères elliptiques. Ceux-ci, caractérisés par une forme oblongue comme leur nom l'indique, sont supposés se former pour un angle d'impact inférieur à 12°
<figure><img src="image_69.png" /></figure>

FIG. 8.3 – Variations latitudinales maximales des distributions taille / fréquence générées sur la Lune et Mars (long terme).

[Bottke et al., 2000b]. Sur la Lune, la probabilité d'un faible angle d'impact est plus élevée au pôle, et l'on calcule que l'on peut attendre un excèdent de cratères elliptiques de 30 % par rapport à l'équateur (par unité de surface). Sur Mars, et sur la moyenne des variations séculaires, on détermine un effet exactement opposé.

Sur la Lune, les variations du taux de production des cratères sont aussi fonction de la longitude, comme le montre la figure 8.4. A la différence des variations en latitude, celles-ci sont sensibles à la distance Terre-Lune. Nous considérons que, au cours des 3 derniers milliards d'années, l'évolution de cette distance s'est faite linéairement entre 40 et 60 rayons terrestres (chapitre 4). L'amplitude de l'effet est légèrement accrue par rapport à la distance Terre-Lune actuelle. Tout comme pour le flux d'impacts, le maximum sur l'ensemble de la surface lunaire est localisé à l'apex, tandis que le minimum, qui résulte de la compétition de l'effet apex / antapex et de l'effet de latitude, se situe aux environs de (±60°N, 90°E). Les variations sont accrues par rapport au flux d'impacts. Le rapport du taux maximal sur le taux minimal est de 1.47.
<figure><img src="image_70.png" /><figcaption>FIG. 8.4 – Variations du taux de formation des cratères de diamètre supérieur à 1 km à la surface de la Lune, estimé pour les trois derniers milliards d'années.</figcaption></figure>

## 8.3 Confrontation aux observations

Nous tentons ici de comparer nos estimations avec les observations. Une comparaison approfondie des distributions taille / fréquence réelles des cratères planétaires recouvrant des unités géologiques d'âges similaires et de positions différentes est en dehors du cadre de ce travail. Nous nous contentons de comparer nos résultats avec les bases de données globales disponibles. Nous ne tentons pas d'utiliser les observations martiennes pour vérifier notre estimation des variations en latitude, car les pôles martiens et l'hémisphère nord, géologiquement jeunes, viendraient fausser la comparaison. La tectonique des plaques et les océans interdisent toute exploitation des données terrestres. La seule comparaison possible est celle du bombardement atmosphérique, et celle-ci a été réalisée (avec succès) au chapitres 5 et 6. Nous n'avons pas connaissance d'une base de données exploitable dans le cas de Mercure. Restent Venus et la Lune.

La population des cratères vénusiens a été observée à partir de l'imagerie SAR de la mission Magellan [Phillips et al., 1992]. On pense de Vénus que sa surface entière s'est reformée simultanément suite à un épisode d'activité volcanique intense, et des tests statistiques ont montré que la distribution des cratères d'impact n'est pas distinguable d'une distribution uniforme [Strom et al., 1994]. Cela est en parfait accord avec nos estimations, qui prédissent une dépendance en latitude quasi-nulle.
<figure><img src="image_71.png" /></figure>

FIG. 8.5 – Comparaison entre les variations en latitude prédites et les observations lunaires.

Dans le cas de la Lune, nous utilisons la base de données<fnref n="1" /> des cratères baptisés, compilée par J. McDowell à partir du travail de Andersson et Whitaker [1982]. Dans la mesure où tous les cratères lunaires n'ont évidemment pas été nommés, nous nous restreignons aux cratères de diamètre supérieur à 25 ou 50 km, en supposant que ceux-ci sont suffisamment visibles pour tous avoir été reconnus et nommés. La résolution des images lunaires se dégradant quelque peu vers les pôles (Tableau 2 de Wilhelms et al. [1987]), la différence entre les deux jeux de données (> 25 et > 50 km) est supposée donner une idée raisonnable de l'influence de la précision variable des observations. L'histoire géologique de la face visible est rendue complexe par la présence des mers basaltiques, plus jeunes que les haut-plateaux environnants. Nous nous restreignons donc aux observations effectuées sur la face cachée, dont on peut raisonnablement supposer qu'elle a été affectée uniformément par les mêmes processus géologiques tout au cours de son histoire. Précisément, la zone exploitée est délimitée par les longitudes 100 E et 250 E, afin d'éviter les bassins Orientale et Australe. Tout ceci laisse un total de 1349 cratères plus grands que 25 km et 539 plus grands que 50 km. Dans le but de

<footnote n="1">¹http://www.planet4589.org/astro/lunar</footnote>
réduire les incertitudes statistiques autant que possible, la moyenne est faite sur les hémisphères nord et sud, ainsi que sur des bandes de 30 degrés en latitude. La distribution résultante est tracée en figure 8.5. On y a également représenté l'ensemble des cratères Coperniciens et Eratosthèniens plus grands que 25 km (127 cratères sur l'ensemble de la Lune), postérieurs aux mers basaltiques et à toute activité géologique majeure. Ceux-ci nous ont été fournis par N. Petro, à partir des données publiées par Wilhelms et al. [1987], Grier et al. [2001], McEwen et al. [1993] et McEwen et Robinson [1997]. On peut voir que la densité des cratères montre une tendance à être plus réduite au pôle qu'à l'équateur, comme nos résultats le prédissent. Néanmoins, les grandes incertitudes statistiques associées au faible nombre de données, ainsi que la rigueur discutable des bases de données utilisées, ne permettent pas d'être affirmatif.

Nous comparons également nos résultats aux données récoltées par Morota et Furumoto [2003] et Morota et al. [2005], qui ont mesuré la densité des cratères étoilés (donc récents et aisément identifiables) en fonction de la distance angulaire à l'apex, $\gamma$. Ces cratères, au regard de la durée de vie des éjectas en étoile, constituent le bombardement lunaire de ce dernier milliard d'années. La comparaison, favorable, est montrée en figure 8.6. On a également représenté l'expression analytique calculée par Morota et al. [2005], selon la forme fonctionnelle
$$ \phi_c(\gamma) = (1 + a \cos \gamma)^b $$
prédite par Zahnle et al. [2001], en ajustant les valeurs de $a$ et de $b$. Ils déduisent de la meilleure valeur de $a$ une vitesse moyenne d'impact d'environ $15 \text{ km.s}^{-1}$. Pour une vitesse moyenne d'impact de près de $20 \text{ km.s}^{-1}$, notre modèle produit une courbe aussi convaincante au regard des données, même si l'effet mesuré semble supérieur à l'effet calculé. De nombreux facteurs peuvent être responsables de cette différence, notamment l'évolution de la distance Terre-Lune, mal contrainte. On rappelle que la forme fonctionnelle établie par Zahnle et al. [2001] est prédite lorsque le bombardement est supposé isotrope. La courbe produite par notre modèle semble mieux reproduire le quasi plateau observable pour $\gamma > \pi/2$, caractéristique d'un effet de latitude.

## 8.4 Troyens autour de la Terre ?

Nous avons pu constater qu'aucune asymétrie face visible / face cachée n'est observable dans nos résultats pour la Lune. Pour qu'un tel effet existe, il faut que l'inclinaison du vecteur vitesse d'approche avec le système Terre-Lune, $i'$, soit quasi nulle, permettant ainsi à la Terre de concentrer les trajectoires sur une
<figure><img src="image_72.png" /></figure>

FIG. 8.6 – Taux de formation des cratères lunaires en fonction de la distance angulaire à l’apex. Comparaison avec les données récoltées par Morota et Furumoto [2003].

ligne qui coupe l'orbite du satellite. La distribution des i' générée par l'ensemble
de la population d'astéroïdes ne permet pas à cette asymétrie de se faire sentir.

Cependant, une population d'astéroïdes telle que d'hypothétiques troyens ter-
restres, piégés aux points de Lagrange puis éjectés, serait susceptible de réunir les
conditions nécessaires à l'existence d'un tel effet. Les travaux de Mikkola et Inna-
nen [1995] et Wiegert et al. [2000] prédissent que ces objets, que l'on n'observe pas
aujourd'hui, pénétreraient préférentiellement dans le système Terre-Lune avec un
vitesse autour de 1-2 kilomètres par seconde et une inclinaison comprise entre 0
et 6 degrés.

La figure 8.7 représente les variations spatiales du taux de formation des
cratères lunaires produites par une population d'objets qui vérifie $U' = 2\text{km.s}^{-1}$
et $i' \in [0°, 6°]$. On a choisi de considérer que $i'$ est uniformément distribué dans cet
intervalle. La distance Terre-Lune est posée à 40 rayons terrestres. On constate
alors que le maximum est situé entre l'apex et le centre de la face visible. Le
ratio hémisphère visible / hémisphère caché est de 1.4. Ce ratio peut atteindre
1.5 pour une distance de 30 rayons terrestres.
Si une telle population a bien existé, à l'instar des troyens joviens, et si elle a massivement contribué, à une époque, au bombardement lunaire, alors l'hémi-sphère visible aura été significativement plus exposé.

<figure><img src="image_73.png" /><figcaption>FIG. 8.7 – Taux de formation des cratères lunaires générés par une population hypothétique de troyens. Distance Terre-Lune de 40 rayons terrestres.</figcaption></figure>
# Chapitre 9

## Révision de la datation par comptage de cratères

Ayant calculé les variations spatiales du taux de cratérisation sur la Lune, nous sommes en mesure de corriger la relation âge / densité de cratères. Cette dernière, bâtie à partir des âges mesurés sur les sites d'échantillonnage des missions lunaires, ne tient en effet pas compte de la position de ces sites. Une fois ceci fait, nous formalisons la méthodologie qui permet, à partir de cette relation, d'attribuer un âge à une zone quelconque de la Lune (en tenant compte encore une fois des variations spatiales), mais également à une zone de n'importe quelle autre planète tellurique, puisque nous avons estimé les taux de production de cratères relatifs à la Lune. Nous donnons pour conclure de nouvelles estimations de certains âges sur la Lune, Mercure, Venus, Mars et la Terre.

### 9.1 Recalibration de la chronologie

Sur la Lune, on rappelle que la relation empirique

$$N_{\mathcal{Q}} (D > 1\text{km}, A) = 5.44 \times 10^{-14} (e^{6.93A} - 1) + 8.38 \times 10^{-4} A, \quad (9.1)$$

établie par Neukum et al. [2001a], lie l'âge A d'une unité géologique (en Ga) à la densité de cratères de diamètre supérieur à 1 km. Les points de calibration de cette courbe proviennent essentiellement des âges absolus mesurés autour des sites d'alunissage des missions Apollo et Luna, représentés en figure 9.1.

A partir des variations spatiales données par la figure 8.4, nous tenons compte de la position à laquelle les densités de cratères ont été mesurées pour convertir celles-ci en densités moyennes sur l'ensemble de la surface lunaire, que nous appelons « équivalent global moyen ». L'équivalent global moyen de la densité de cratères mesurée à une position $(\lambda, \varphi)$ donnée est

$$N_{\mathcal{Q}} (D > 1) = N_{\mathcal{Q}} (D > 1, \lambda, \varphi) / \phi_c^{\mathcal{Q}} (D > 1, \lambda, \varphi), \quad (9.2)$$
<figure><img src="image_74.png" /><figcaption>FIG. 9.1 – Sites d'alunissage des missions Apollo et Luna.</figcaption></figure>

avec $\phi_c^c (D > 1, \lambda, \varphi)$ donné par la figure 8.4.

Les points de calibration corrigés sont indiqués sur la figure 9.3, et listés en Table 9.1. La position des unités géologiques correspondantes est indiqué en figure 9.2. A défaut de meilleure estimation, la correction a été appliquée au delà de 3.5 milliards d'années, même si la population d'astéroïdes est soupçonnée d'avoir été très différente. On constate que sur l'ensemble des points, la différence n'est pas flagrante, la plupart des unités géologiques concernées étant éloignées des zones d'extrema du taux de cratérisation, ou suffisamment étendues pour qu'en moyenne, les variations spatiales à leur surface soient proches de la moyenne planétaire. La forme analytique déterminée par Neukum et al. [2001a] n'est pas remise en cause. Au contraire, les points corrigés, en tout cas pour la période où ils peuvent être considérés comme relativement fiables (0–~3.5 Gyr), semblent, assez nettement, en meilleur accord avec l'expression analytique que les points préexistant. On continuera donc à utiliser l'expression donnée par l'équation (1.1) comme base de la méthode de datation. Notons que le point 4 sur la figure 9.3 correspond à l'équivalent lunaire de l'accumulation des cratons terrestres durant les quatre derniers millions d'années. Leur cas est traité plus loin, avec pour résultat un parfait accord avec le modèle, contrairement aux résultats de Neukum et Ivanov [1994].
<figure><img src="image_75.png" /><figcaption>FIG. 9.2 – Mers de lave lunaires.</figcaption></figure>

Les variations maximales sont également représentées sur la figure 9.3, pour des sites situés respectivement aux coordonnées (0°N, 270°E) et (±60°N, 90°E), ainsi qu'un exemple de l'erreur que cela peut engendrer lors de l'estimation d'un âge (rectangle grisé). Puisque la relation est linéaire pour les trois derniers milliards d'années (si le bombardement a bien été constant), une variation spatiale donnée se reporte directement sur un biais en âge dans la même proportion. Supposons que deux unités géologiques, situées aux coordonnées correspondant aux variations maximales, se soient vues attribuer un âge identique, par exemple 2 milliards d'années, sans tenir compte de leur position. Les variations spatiales étant au maximum de 40 %, l'écart en terme d'âge sera lui aussi de 40% de la valeur estimée : 800 millions d'années séparent alors les deux surfaces.

## 9.2 Méthodologie révisée

Puisque la relation entre âge et densité de cratères lunaires reste identique à la relation établie par Neukum et al. [2001a], notre fonction $\eta(t)$ (qui traduit l'évolution temporelle du flux de projectiles et est posée égale à 1 de nos jours)
<figure><img src="image_76.png" /><figcaption>FIG. 9.3 – Correction des points de calibration de la relation âge / densité de cratères, établie sur la Lune.</figcaption></figure>

s'écrit

$$ \eta(t) = \text{Const.} \times \frac{dN_{\mathcal{Q}} (> 1, t)}{dt} = \frac{a \ b}{c} e^{b \ t} + 1, \quad (9.3) $$

avec $a = 5.44 \times 10^{-14}$, $b = 6.93$, $c = 8.38 \times 10^{-4}$.

La distribution taille / fréquence des cratères qui recouvrent une surface d'âge A et de position $(\lambda, \varphi)$ est donc

$$ \begin{aligned} N(> D, A, \lambda, \varphi) &= \phi_c(> D, \lambda, \varphi) N(> D, A), \\ &= \phi_c(> D, \lambda, \varphi) R_c(> D) N_{\mathcal{Q}}(> D, A). \end{aligned} \quad (9.4) $$

où $\phi_c(> D, \lambda, \varphi)$ est la variation relative du taux de cratérisation et où l'on s'est ramené à l'équivalent lunaire global par le biais du taux de cratérisation relatif, $R_c$ (figure 7.2). Il arrive souvent que, par commodité, seule la densité de cratères de diamètres supérieurs à 1 km, $N(D > 1, A, \lambda, \varphi)$, soit déterminée par l'observation. L'estimation de l'âge s'appuie alors sur cette seule valeur, telle que

$$ N_{\mathcal{Q}}(D > 1, A) = \frac{N(D > 1, A, \lambda, \varphi)}{\phi_c(D > 1, \lambda, \varphi) R_c(D > 1)}, \quad (9.5) $$

avec $\phi_c(D > 1, \lambda, \varphi)$ donné par les figures 8.1 et 8.4 (pour les planètes, $\phi_c$ ne dépend pas de $\varphi$, mais cette variable est conservée ici pour que les relations
<table><caption>TAB. 9.1 – Conversion des densité de cratères lunaires données par [Stöffler et Ryder, 2001] en équivalent global moyen, par correction spatiale.</caption><thead><tr><th rowspan="2">Unité</th><th rowspan="2">Localisation</th><th>Diamètre</th><th>Age</th><th>Densité de cratères<br>[Stöffler et Ryder, 2001]</th><th>Equivalent global moyen<br>(cette étude)</th></tr><tr><th>(km)</th><th>(Ga)</th><th>N(&gt; 1, λ, φ)<br>(10<sup>-4</sup> km<sup>-2</sup>)</th><th>N(&gt; 1)<br>(10<sup>-4</sup> km<sup>-2</sup>)</th></tr></thead><tbody><tr><td>1 Cratère Cone</td><td>-3.7°N, -17.5°E</td><td>0.34</td><td>0.025</td><td>0.21</td><td>0.188</td></tr><tr><td>2 Cratère North Ray</td><td>-9°N, -15.6°E</td><td>1</td><td>0.053</td><td>0.44</td><td>0.436</td></tr><tr><td>3 Cratère Tycho</td><td>-43.4°N, -11.1°E</td><td>102</td><td>0.109</td><td>0.90</td><td>0.868</td></tr><tr><td>5 Copernicus</td><td>9.7°N, -20.1°E</td><td>93</td><td>0.80</td><td>13</td><td>12.1</td></tr><tr><td>6 O. Procellarum (A12)</td><td>18.4°N, -57.4°E</td><td>2568</td><td>3.15</td><td>36</td><td>30.0</td></tr><tr><td>7 M. Crisium (L24)</td><td>17.0°N, 59.1°E</td><td>418</td><td>3.22</td><td>30</td><td>33.8</td></tr><tr><td>8 M. Imbrium (A15)</td><td>32.8°N, -15.6°E</td><td>1123</td><td>3.30</td><td>32</td><td>31.3</td></tr><tr><td>9 M. Fecunditatis (L16)</td><td>-7.8°N, 51.3°E</td><td>909</td><td>3.41</td><td>33</td><td>39.0</td></tr><tr><td>10 M. Tranquilitatis (jeune) (A11)</td><td>8.5°N, 31.4°E</td><td>873</td><td>3.58</td><td>64</td><td>67.0</td></tr><tr><td>11 Orientale Basin / ejecta</td><td>-19.4°N, -92.8°E</td><td>930</td><td>3.73</td><td>110</td><td>93.3</td></tr><tr><td>12 M. Serenitatis (A17)</td><td>28.0°N, 17.5°E</td><td>707</td><td>3.75</td><td>100</td><td>101.8</td></tr><tr><td>13 M. Tranquilitatis (ancienne) (A11)</td><td>8.5°N, 31.4°E</td><td>873</td><td>3.80</td><td>90</td><td>94.3</td></tr><tr><td>14 Fra Mauro (A14)</td><td>-3.7°N, -17.5°E</td><td></td><td>3.77</td><td>370</td><td>331.4</td></tr><tr><td>15 Descartes (A16)</td><td>-9.0°N, 15.6°E</td><td></td><td>3.92</td><td>340</td><td>336.7</td></tr><tr><td>16 Crisium Basin</td><td>17.0°N, 59.1°E</td><td>740</td><td>3.84</td><td>570</td><td>643.8</td></tr><tr><td>17 Nectaris Basin</td><td>-15.2°N, 35.5°E</td><td>333</td><td>3.92</td><td>1200</td><td>1265.3</td></tr></tbody></table>

soient valides dans le cas lunaire – en posant $R_c = 1$). L'âge est alors simplement donné à partir de $N_c(D > 1, A)$ par la relation (9.1).

Si cette manière de procéder est très commode, elle reste beaucoup moins robuste que la recherche du meilleur accord entre une distribution taille-fréquence prédite et la distribution observée sur un grand intervalle de taille. Cette dernière approche est plus fastidieuse, car elle nécessite de trier avec précision les cratères en fonction de leur diamètre $D$. Qui plus est, les grandeurs $\phi_c$ et $R_c$, fonctions de $D$, doivent être connues à chacun de ces diamètres. Cependant, les observations, localisées, reproduisent rarement la distribution idéalement prédite, pour différentes raisons. Notamment, des variations purement statistiques existent. Certains processus géologiques sont également susceptibles d'avoir modifié (par érosion, recouvrement, etc.) le nombre de cratères visibles dans une gamme de taille préférentielle. Se limiter à $D > 1$ km revient à ne pas s'assurer que la distribution de cratère se comporte « normalement », et induit le risque de s'appuyer sur une mesure biaisée. Sachant que $N_c(D > D, A) = \Phi_c^c(D > D) \int_0^A \eta(t) dt$, on écrit à partir de l'équation (9.4) :

$$N(D > D, A, \lambda, \varphi) = \phi_c(D > D, \lambda, \varphi) R_c(D > D) \Phi_c^c(D > D) \left( \frac{a}{c} e^{bA} + A \right), \quad (9.6)$$

où $\Phi_c^c(D > D)$ est le taux de cratérisation actuel et absolu, en moyenne sur l'ensemble de la surface lunaire, tel que $\Phi_c^c(D > D) = \Phi^c(D > 1) \bar{S}(D > d(D))$. Avec un flux de projectile calculé comme $\Phi^c(D > 1) = 2.1 \times 10^{-15} \text{km}^{-2}\text{an}^{-1}$ et la
grandeur $\bar{S}(> d(D))_c$ donnée par l'équation (7.24), on dispose de toutes les informations nécessaires au tracé de la distribution taille / fréquence, c'est-à-dire $N(> D, A, \lambda, \varphi)$ en fonction de $D$. Faire varier l'âge $A$ revient alors à translater verticalement la courbe sur le graphe. On cherche donc $A$ tel que la distribution prédite reproduise au mieux les données.

## 9.3 Attribution de nouveaux âges

Pour finir, nous appliquons la méthode décrite à l'instant pour estimer l'âge de diverses unités géologiques lunaires et planétaires. L'ensemble des âges recalculés est listé en table 9.2, et comparé avec diverses valeurs publiées au préalable.

<table><caption>TAB. 9.2 – Ages estimés. Valeurs comparées à : <sup>a</sup>Neukum et al. [2001b]; <sup>b</sup>Strom et Neukum [1988]; <sup>c</sup>McKinnon et al. [1997]; <sup>d</sup>Neukum et Ivanov [1994]; <sup>e</sup>Strom et al. [1994]; <sup>f</sup>Grieve et Dence [1979]; <sup>g</sup>Hartmann et Neukum [2001].</caption><thead><tr><th rowspan="2">Planète</th><th rowspan="2">Unité</th><th>Age</th><th>Age</th></tr><tr><th>(littérature)</th><th>(cette étude)</th></tr></thead><tbody><tr><td rowspan="2">Mercure</td><td>Mare Caloris</td><td>3.77<sup>a</sup>, 3.85<sup>b</sup></td><td>3.70 Ga</td></tr><tr><td>Highlands</td><td>4.18<sup>a</sup>, 4.07<sup>b</sup></td><td>3.89 Ga</td></tr><tr><td>Vénus</td><td>surface globale</td><td>750<sup>c</sup>, 650<sup>d</sup>, 290<sup>e</sup></td><td>230 Ma</td></tr><tr><td>Terre</td><td>Cratons phanérozoïques</td><td>700<sup>d</sup>, 375<sup>f</sup></td><td>380 Ma</td></tr><tr><td rowspan="2">Mars</td><td>Amazonis Planitia D &lt; 1km</td><td>100 – 200<sup>g</sup></td><td>100 – 200 Ma</td></tr><tr><td>Amazonis Planitia D &gt; 1km</td><td>0.5 – 2<sup>g</sup></td><td>~ 1 Ga</td></tr></tbody></table>

Sur Terre, le temps nécessaire à l'accumulation des cratons Nord américains et européens a été estimé géologiquement à environ 375 millions d'années par Grieve et Dence [1979]. A partir de la distribution taille / fréquence établie par ces auteurs, Neukum, G. [1983], qui établit pour la première fois une relation âge / densité relativement robuste pour la Lune, exporte la méthode de datation sur Terre et en déduit un âge relativement proche, qui vient conforter l'hypothèse d'un bombardement constant dans les dernières centaines de millions d'années (point 4 de la figure 9.3). Cependant, Neukum et Ivanov [1994] mettent en doute la validité de ce résultat. En effet, le passage des conditions lunaires aux conditions terrestres est fait de manière très approximative : il ne tient pas compte de la gravité de la Terre. En appliquant un traitement plus rigoureux, Neukum et Ivanov [1994]
<figure><img src="image_77.png" /><figcaption>FIG. 9.4 – Distribution des cratons terrestres et âge donné par notre modèle.</figcaption></figure>

constatent que l'âge produit par la méthode de datation est de ~ 700 Ma, en désaccord d'un facteur 2 avec l'âge (robuste) donné par Grieve et Dence [1979]. Ils concluent en évoquant la possibilité d'un accroissement du bombardement bref et intense dans les temps géologiques récents. Pour autant, les auteurs conservent le point de calibration établi par Neukum, G. [1983]. La figure 9.4 représente la distribution taille / fréquence des cratons, et la distribution synthétique, calculée par nos soins pour un âge qui reproduit au mieux les données. Le modèle donne un âge de 380 Ma, en parfait accord avec la valeur de Grieve et Dence [1979]. Ce résultat est un bon indicateur de confiance pour notre modèle. Notons que le déficit de cratères de diamètre inférieur à 10 km est attribué à l'érosion.

Il est généralement admis que l'ensemble de la surface de Venus s'est renou-
velée à peu près simultanément suite à un épisode de volcanisme intense. La
figure 9.5 représente la distribution de l'ensemble des cratères vénusiens obser-
vés, accompagnée de notre distribution modèle. Nous estimons que la surface de
Vénus est vieille de 230 Ma. De nombreuses tentatives ont été faites dans la lit-
térature pour déterminer cet âge, avec des estimations allant de 290 à 800 Ma,
pour des hypothèses très variées. Une revue assez complète peut être trouvée
dans McKinnon et al. [1997]. Ces derniers incluent un modèle atmosphérique qui
reproduit très fidèlement la distribution de cratères observée, même aux petites
tailles, et estiment l'âge à 700 Ma. Cependant, ils utilisent une estimation de la
<figure><img src="image_78.png" /><figcaption>FIG. 9.5 – Distribution des cratères vénusiens et âge donné par notre modèle.</figcaption></figure>

population d'astéroïdes et de comètes faite par Shoemaker et al. [1990], qui se base sur des observations bien moins complètes qu'aujourd'hui. Notre modèle, basé sur une distribution orbitale et une distribution en taille des projectiles plus récente, reproduit très bien la distribution des cratères pour lesquels $D >\sim 10$ km, semblant indiquer que l'atmosphère ne freine pas le processus de cratération au-delà de cette taille. Cependant, si l'hypothèse que les petits cratères se forment en régime poreux ne s'applique pas à Vénus, le seuil se déplace à $\sim 30$ km. Dans leurs simulations, McKinnon et al. [1997] prédisent que l'effet de l'atmosphère se fait notablement sentir jusqu'à $D \simeq 80$ km. Cependant, les parameètres qui contrôlent la traversée atmosphérique (coefficients d'ablation, seuil de fragmentation, etc.), mal contraints, ont été ajustés de manière à ce que les données soient reproduites par la population de projectile choisie par les auteurs, très différente de la nôtre.

La figure 9.6, qui concerne Mercure, donne l'estimation de l'âge des Highlands et de Caloris Mare. Les données sont représentées sous forme de distribution relative, qui permet de mieux observer les variations selon $D$. L'âge des highlands donné par le modèle est 3.89 Ga, contre 4.07 pour Neukum et al. [2001b] et 4.18 pour Strom et Neukum [1988]. L'âge de Caloris Mare est calculé à 3.70 Ga, contre 3.77 pour Neukum et al. [2001b] et 3.85 pour Strom et Neukum [1988]. Puisque notre distribution en taille de projectiles suit une succession de lois de puissances,
<figure><img src="image_79.png" /></figure>

FIG. 9.6 – Distribution des cratères de Caloris Mare et des Highlands de Mercure, et âges donnés par notre modèle.

la distribution cumulée des cratères fait de même. Cela produit d'abruptes variations de la distribution relative. Pour autant, la distribution synthétique reproduit bien les données, peut-être même mieux que l'expression polynomiale de Neukum et al. [2001b] (voir la figure 1.5 (a)).

Enfin, la figure 9.7 correspond à Amazonis Planitia, sur Mars. Il semble que les grands cratères soient plus anciens, puisque la distribution synthétique ne peut reproduire convenablement les données sur l'ensemble de la gamme de diamètres. C'est en tout cas ce que proposent Hartmann et Neukum [2001], en suggérant que l'épanchement de lave qui s'étend sur Amazonis Planitia, mince et récent, est recouvert par les petits cratères (le jeune âge de la coulée laisse peu de chances à de grands cratères, moins fréquents, de se former), tout en laissant apparaître les grands cratères préexistants recouverts. Bien sûr, il est possible que le modèle ne soit pas adapté dans ce cas, pour une raison ou pour une autre (autre population de projectiles, autre régime de formation des cratères, etc.). Quoiqu'il en soit, l'âge est estimé à environ 1 Ga pour $D \gg 1$km, et 100 – 200 Ma pour les cratères plus petits, c'est-à-dire pour la coulée de lave. Hartmann et Neukum [2001] donnent une estimation équivalente.
<figure><img src="image_80.png" /><figcaption>FIG. 9.7 – Distribution des cratères de Amazonis Planitia, sur Mars, et âges estimés par notre modèle.</figcaption></figure>
Conclusions

En nous appuyant sur les derniers modèles et les plus complètes observa-
tions de la population d'astéroïdes, nous avons, au cours de ce travail, modélisé
le bombardement des planètes telluriques et de la Lune. Le chemin qui sépare
l'astéroïde du cratère d'impact est long et parsemé d'incertitudes. Pourtant, un
certain nombre de résultats présentés sont confirmés par des observations variées.
Les progrès constants dans le domaine de l'observation spatiale, de la dynamique
du système solaire et de la mécanique de formation des cratères permettront à
court terme d'affiner les estimations établies dans ce manuscrit.

La distribution orbitale et la distribution en taille de la population d'asté-
roïdes et de comètes ont été bâties en faisant la synthèse de nombreux travaux
indépendants [Bottke et al., 2002; Morbidelli et al., 2002; Stuart et Binzel, 2004;
Halliday et al., 1996; ReVelle, 2001; Brown et al., 2002; Harris, A. W., 2002; Rabi-
nowitz et al., 2000]. Les fomules d'Öpik [Öpik, 1951; Wetherill, 1967; Greenberg,
1982; Shoemaker et Wolfe, 1982] ont été utilisées pour déterminer les probabilités
d'impact avec les planètes, et adaptées au cas de la Lune. Les derniers résultats
concernant l'évolution temporelle des orbites planétaires ont été inclus dans les
estimations Laskar [2008]; Laskar et al. [2004]; Eriksson et Simpson [2000]; Sonett
et Chan [1998]; Touma et Wisdom [1994]. Une méthode analytique a été déve-
loppée pour calculer le bombardement comme fonction de la position à la surface
des planètes, à partir de la distribution des conditions d'approche déterminées
par le calcul des probabilités. La conversion du flux d'impacts en taux de formac-
tion de cratères a été rendue possible par l'usage des lois empiriques établies en
laboratoire [Melosh, 1989; Schmidt et Housen, 1987; Holsapple et Housen, 2007].

L'évolution journalière du bombardement terrestre est en remarquable ac-
cord avec les données radar d'entrées micro-météoritiques dans l'atmopshère
[Campbell-Brown, 2008; Jones et al., 2005; Webster et al., 2004]. Cela suggère
fortement que les météoroïdes reproduisent le comportement orbital de leurs as-
téroïdes parents, lesquels, plus facilement observables, ont servi à bâtir le modèle.
Sur l'ensemble de la population de croiseurs, on observe deux maxima du flux d'impacts terrestre : l'un à midi, l'autre à minuit. A mesure que l'on considère des projectiles de plus en plus rapides, le flux se concentre autour de 6 heures du matin.

Le pente de la distribution amplitude / fréquence des enregistrements sismiques des stations lunaires Appolo (corrigée des séismes profonds) est très bien reproduite par le modèle, suggérant qu'il est apte à prédire le bombardement micro-météoritique du satellite. Nous déterminons un bruit de fond météoritique qui varie entre deux centièmes et un dixième du pouvoir de résolution des sismomètres Appolo.

Sur la Lune, nous calculons qu'un cratère de diamètre supérieur au kilomètre se forme tous les 3500 ans. Cela est en accord avec le taux de cratérisation établi par Neukum et al. [2001a] pour les trois derniers milliards d'années. La distribution taille-fréquence des cratères lunaires, modélisée à partir de notre distribution d'astéroïdes, reproduit de manière excellente celle de Neukum et al. [2001a], donnée valide pour les 4 derniers milliards d'années. Il est nécessaire pour cela de considérer que les petits cratères (diamètre transitoire inférieur à 5 km) se forment en régime poreux, comme suggéré par Ivanov [2006]. Cela revient à considérer que le premier kilomètre de la croûte possède une porosité importante, ce qui semble en accord avec les faibles vitesses sismiques déterminées dans cette zone [Lognonné et al., 2003]. La réconciliation entre les distributions observées d'astéroïdes et de cratères est possible sans l'inclusion dans le modèle des cratères secondaires de tailles subkilométriques. Cela demande une analyse plus poussée [McEwen et Bierhaus, 2006].

Nous calculons le taux de cratérisation relatif entre les planètes telluriques et la Lune, en moyenne sur les derniers milliards d'années. Celui-ci est, par unité de surface et pour des cratères de diamètre supérieur au kilomètre, et sans tenir compte d'une éventuelle atmosphère, de 2.8 pour Mercure, 0.9 pour Vénus, 0.6 pour la Terre et 0.8 pour Mars. Pour des cratères de taille supérieure à 30 kilomètres, il est de 3.5, 1.7, 1.2 et 1.2, respectivement.

Nous prédisons des variations spatiales significatives du bombardement à la surface de la Lune et de Mars. Les pôles martiens, une fois les variations d'obliquité de la planète prises en compte, accumulent les cratères 30 % plus rapidement que l'équateur (par unité de surface). Sur la Lune, des variations en longitude s'ajoutent, engendrées par sa rotation synchrone. Le minimum est localisé en (±60°N, 90°E), et le maximum, près de 50 % plus élevé, se situe en (0°N, 90°W).
Les variations longitudinales lunaires sont validées par le comptage des cratères étoilées [Morota et Furumoto, 2003; Morota et al., 2005]. Les variations en latitude nécessitent d'être confrontées à de plus amples observations, mais une tendance semble être visible dans les jeux de données utilisés, même si les incertitudes statistiques sont grandes. Ces résultats ont d'importantes implications dans la recherche des meilleurs sites d'implantation de futures bases lunaires ou martiennes. Ils ont aussi une conséquence directe sur les âges prédits par la méthode de datation par comptage de cratères [Neukum et Ivanov, 1994; Neukum et al., 2001a; Hartmann et Neukum, 2001], qui considère que le bombardement est isotrope. Si les incertitudes inhérentes à la méthode sont en général plus grandes que les erreurs induites par les variations spatiales, le biais associé peut être de 800 millions d'années dans le pire des cas.

Nous donnons de nouvelles estimations de certains âges géologiques. En particulier, le temps d'accumulation des cratons terrestres est estimé à 380 millions d'années, un âge extrêmement proche de celui donné géologiquement par Grieve et Dence [1979], à savoir 375 Ma. Cet âge ne pouvait jusqu'ici être retrouvé par la méthode de datation (Neukum et Ivanov [1994] donnent 700 Ma). Ce résultat procure une certaine confiance dans le modèle. La surface de Vénus est estimée vieille de 230 Ma, un âge plus jeune que l'ensemble des estimations données jusqu'à présent (McKinnon et al. [1997] donnent 700 Ma). Les hautes terres martiennes sont d'après nos résultats plus jeunes de 200 millions d'années que ne l'estime Neukum et al. [2001b]. Les résultats obtenus pour Mars sont sensible-ment identiques à ceux proposés par Hartmann, W. K. [2005], tant que les zones datées ne sont pas proches du pôle.

Les méthodes développées ou reprises dans ce manuscrit permettent de modéliser rapidement le bombardement d'une planète quelconque, et procurent un outil facile d'utilisation. La prolongement de ce travail verra la confrontation plus poussée des estimations obtenues à l'aide de ces méthodes avec les données géologiques et spatiales.
# Annexe A

## Trajectoires hyperboliques générées par une phase d'approche

Au chapitre 5, nous déterminons la probabilité d'impact avec une cible donnée (planète, lune, voire système planète-lune) en fonction des conditions d'approche $U_\infty, i_\infty, j_\infty$. Chaque approche décrite par une triade ($U_\infty, i_\infty, j_\infty$) génère une nuée de projectiles uniformément distribués sur un disque situé à l'infini de la cible (figures 5.2 et 5.6). La position sur ce disque est donnée par le paramètre d'impact $b$ et l'angle $\delta$ (figure 5.3). L'objet de cette annexe est d'exprimer les éléments orbitaux d'un projectile au cours de la phase d'approche, dans le repère orbital (OXYZ) centré sur la cible, en fonction des conditions d'approche et de la position initiale sur le disque. Les relations dérivées ici sont appliquées à deux cas distincts abordés dans cette étude : (1) la description de la trajectoire du projectile à l'approche de la cible (planète ou satellite) afin de déterminer les coordonnées de l'impact ; (2) la description des orbites qui pénètrent dans un système planète-satellite, en vue de déterminer la probabilité d'impact avec le satellite (chapitre 6, section 2.2).

On a défini le repère ($Ox_1x_2x_3$) comme indiqué par la figure 5.3. Quand les angles d'approche sont tous deux nuls, le repère ($Ox_1x_2x_3$) se confond avec (OXYZ). Dans ($Ox_1x_2x_3$), le vecteur position sur le disque s'exprime

$$ \mathbf{r}_{\infty}(Ox_1x_2x_3) = \begin{pmatrix} \infty \\ b \cos \delta \\ b \sin \delta \end{pmatrix}, \qquad (A.1) $$

tandis que le vecteur vitesse initial est donné par

$$ \mathbf{U}_{\infty}(Ox_1x_2x_3) = \begin{pmatrix} -U_{\infty} \\ 0 \\ 0 \end{pmatrix}. \qquad (A.2) $$
Le passage au repère orbital se fait par l'intermédiaire de la matrice P telle
que

$$
P_{(Ox_1x_2x_3) \to (OXYZ)} =
\begin{bmatrix}
\cos i_{\infty} & 0 & -\sin i_{\infty} \\
0 & 1 & 0 \\
\sin i_{\infty} & 0 & \cos i_{\infty}
\end{bmatrix}
\times
\begin{bmatrix}
\cos j_{\infty} & -\sin j_{\infty} & 0 \\
\sin j_{\infty} & \cos j_{\infty} & 0 \\
0 & 0 & 1
\end{bmatrix}
. \tag{A.3}
$$

Il faut noter que l’uniformité sur le disque se traduit par le fait que $\delta$ est uniformément distribué entre 0 et $2\pi$, mais non $b$. C’est $b^2$ qui est uniformément distribué entre 0 et $b_{max}^2$. Dans l’hypothèse où l’on souhaite générer aléatoirement une nuée de projectile, ce qui est un des moyens de simuler le bombardement (une alternative plus coûteuse en temps de calcul que la méthode décrite au chapitre 7), il faut échantillonner $\delta$ et $b$ comme

$$
\begin{equation}
\begin{cases}
\delta = x_1 \ 2\pi, \\
b^2 = x_2 b_{\text{max}}^2,
\end{cases}
\tag{A.4}
\end{equation}
$$

où $(x_1, x_2) \in ]0, 1[$ sont des nombres aléatoires uniformément distribués $^1$.

A.0.1 Hyperboles d'impact

On rappelle que la distance à la cible est donnée par

$$
r = \frac{p}{1 + e \cos f}. \tag{A.5}
$$

La constante des aires s'exprime comme (Éq. (2.4),

$$
h = |\mathbf{r}_{\infty} \times \mathbf{U}_{\infty}| = b U_{\infty} . \tag{A.6}
$$

Le paramètre de la conique est donc (Éq. (2.6),

$$
p = \frac{b^2 U_{\infty}^2}{G M_0}, \tag{A.7}
$$

où $M_0$ est la masse de la cible. D'après l'équation (2.22), qui relie l'énergie à
l'excentricité, cette dernière s'exprime (avec l'énergie donnée par l'équation (5.9))
comme

$$
e = \sqrt{1 + \frac{b^2 U_{\infty}^4}{(GM_0)^2}} . \qquad (A.8)
$$

<footnote n="1">¹En effet, échantillonner une grandeur revient à identifier sa distribution cumulée de probabilité avec un nombre aléatoire. En coordonnées cartésiennes, l'uniformité sur le disque se traduit par une probabilité constante, $p(x, y) = A$. Grâce à la matrice jacobienne qui relie le repère cartésien $(x, y)$ au repère polaire $(b, \delta)$, on déduit $p(b, \delta) = b p(x, y) = b A$. La probabiltié (normalisée) d'avoir un angle plus petit que $\delta$ pour une valeur quelconque de $b$ est $p(< \delta) = \int_0^\delta p(b, \delta) \, d\delta / \int_0^{2\pi} p(b, \delta) \, d\delta = \delta / 2\pi$. De la même façon, $p(< b) = b^2 / b_{\text{max}}^2$.</footnote>
Au point d'impact, on a $r = R_0$. Le cosinus de l'anomalie vraie correspondant à l'impact, $f_i$, est donc donné par

$$C_i = \frac{p/R_0 - 1}{e}, \qquad (A.9)$$

où $R_0$ est le rayon de la cible. Initialement, sur le disque, la distance qui sépare le projectile et la planète est infinie. D'après l'équation (A.5), le cosinus de l'anomalie vraie initiale, $f_\infty$, vérifie

$$C_\infty = -1/e. \qquad (A.10)$$

L'anomalie vraie à l'infini est négative (vérifiable par le fait que $\mathbf{r}_\infty \cdot \mathbf{U}_\infty < 0$, cf. chapitre 2). Puisque la collision a nécessairement lieu avant ou au péricentre ($f=0$), c'est aussi le cas pour l'anomalie vraie à l'impact. On a donc

$$S_i = -\sqrt{1 - C_i^2} \quad \text{et} \quad S_\infty = -\sqrt{1 - C_\infty^2}. \qquad (A.11)$$

Il y a impact quand le péricentre de l'hyperbole est inférieur au rayon de la cible, soit $q \le R_0$. On vérifie aisément en exprimant $q = p/(1+e)$ (éq. (2.16)) en fonction de $b$ à l'aide des équations (A.7) et (A.8), que cela correspond à $b \le \tau$. En d'autres termes, le rayon du disque ($\sigma$), $b_{max}$, est le rayon de la section gravitationnelle de la cible, $\tau$ (Eq. (5.10)). Quand les projectiles ont une vitesse relative suffisamment élevée au regard de la masse de la cible, leur trajectoire s'assimile à une ligne droite et $\tau$ tend vers $R_0$. Dans le cas extrême opposé, la masse de la cible est très élevée par rapport à l'énergie cinétique d'approche, et $\tau$ tend vers l'infini, c'est-à-dire que l'hyperbole de la trajectoire peut s'incurver sans limite pour permettre aux projectiles les plus éloignés d'entrer en collision avec la cible.

Considérons à ce stade un nouveau repère, $(On_1n_2n_3)$, tel que $\mathbf{n}_2$ est colinéaire à $\mathbf{b}$, et qui est confondu à $(Ox_1x_2x_3)$ quand $\delta = 0$. La trajectoire du projectile appartient donc au plan $(On_1n_2)$. La matrice de passage est

$$P'_{(On_1n_2n_3) \to (Ox_1x_2x_3)} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \delta & -\sin \theta \\ 0 & \sin \delta & \cos \theta \end{bmatrix}. \qquad (A.12)$$

Les coordonnées du point d'impact dans $(On_1n_2n_3)$ sont

$$\mathbf{r}_{i(On_1n_2n_3)} = R_0 \begin{pmatrix} -\sin(f_i - f_\infty) \\ -\cos(f_i - f_\infty) \\ 0 \end{pmatrix} = R_0 \begin{pmatrix} -S_i C_\infty + S_i C_\infty \\ -(C_i C_\infty + S_i S_\infty) \\ 0 \end{pmatrix}, \qquad (A.13)$$
En conséquence, les coordonnées de l'impact dans le repère orbital (OXYZ) sont

$$ \mathbf{r}_{i(\text{OXYZ})} = P \times P' \times \mathbf{r}_{i(\text{On}_1\text{n}_2\text{n}_3)} \quad (\text{A.14}) $$

et ne dépendent, par l'utilisation successive des équations (A.7)–(A.11), (A.13), (A.12) et (A.3), que de $\mathbf{U}_\infty$, $b$ et $\delta$. Si nécessaire, on tient compte de l'obliquité de la cible et de sa rotation en faisant intervenir un nouveau repère mobile pour exprimer les coordonnées de l'impact selon la latitude et la longitude géographique. Il est alors possible de simuler le bombardement en générant un grand nombre de projectiles par tirage aléatoire de $b$ et $\delta$ selon l'équation (A.4) pour chaque $\mathbf{U}_\infty$ (en pondérant la contribution de chaque nuée de projectiles approchant selon $\mathbf{U}_\infty$ par sa probabilité, telle qu'elle à été déterminée au chapitre 5). Cette approche a été utilisée dans la publication figurant en annexe D. Une approche moins coûteuse en temps de calcul, qui exprime le flux d'impact de façon analytique en utilisant les équations (A.7)–(A.11), est présentée au chapitre 6.

La vitesse d'impact $u$ est donnée par la conservation de l'énergie. On a

$$ \frac{U_{\infty}^{2}}{2} = \frac{u^{2}}{2} - \frac{GM_{0}}{R_{0}}, \qquad (\text{A.15}) $$

soit

$$ u = \sqrt{U_{\infty}^{2} + v_{e}^{2}}. \qquad (\text{A.16}) $$

Enfin, il est utile de déterminer l'angle d'incidence de l'impact par rapport à l'horizontale, $\theta$. On déduit de la conservation du moment angulaire spécifique $\mathbf{h} = \mathbf{r} \times \mathbf{v}$ que

$$ h = \|\mathbf{r}_i \times \mathbf{u}\| = R_0 u \sin(\theta + \pi/2), \quad (\text{A.17}) $$

soit, d'après l'équation (A.6),

$$ \cos \theta = b/\tau, \qquad (\text{A.18}) $$

où $\theta \in [0, \pi/2]$.

### A.0.2 Hyperboles dans le système planète-satellite

Nous nous plaçons ici dans le cas où la « cible » est un système planète-satellite, typiquement le système Terre-Lune. Une fois que le projectile entre dans ce système, on considère qu'il n'est soumis qu'à l'influence gravitationnelle de la planète, et que par conséquent il décrit une hyperbole susceptible de le mener à percuter la planète ou sa lune. Les conditions d'approche avec le système global
sont notées $U'$, $i'$ et $j'$, pour les distinguer de $U_\infty$, $i_\infty$ et $j_\infty$, qui s'appliquent
aux planètes. Comme dans le cas précédent, chaque approche décrite par un $\mathbf{U}'$
donné génère une nuée de projectiles uniformément distribués sur un disque ($\sigma$).
Dans le cas précédent, $b_{max}$ correspond à la limite à laquelle le projectile entre
en collision avec la cible, au-delà il se contente de la frôler, en-deça il ne peut
que la percuter. Ici, $b_{max}$ correspond à la plus grande valeur de $b$ permettant un
impact avec la lune, mais pour $b < b_{max}$, l'impact n'est pas garanti : il possède une
certaine probabilité, et n'est possible que pour certaines valeurs de $b \in [0, b_{max}]$
quand $\mathbf{U}'$ et $\delta$ sont fixés. C'est précisément en vue de déterminer cette probabilité
(chapitre 5, section 3) que nous établissons un certain nombre de relations ici. Une
fois les probabilités d'impact avec la lune obtenues, on se ramène au cas précédent,
où de nouvelles trajectoires hyperboliques sont générées à l'approche de la lune,
qui est alors considérée comme le seul corps massique. Ici, nous décrivons en
fonction de $\mathbf{U}'$, $b$ et $\delta$ les trajectoires hyperboliques générées pendant la phase
intermédiaire où la planète domine, après la phase où les orbites sont déterminées
par le Soleil, et avant la phase d'approche avec la lune. Commençons par établir
$b_{max}$.

Exprimons pour cela l'excentricité de l'hyperbole en faisant intervenir la vi-
tesse orbitale de la lune $v_l$. On considère que $v_l = \sqrt{\frac{GM_p}{a_l}}$, où $M_p$ est la masse de la
planète et $a_l$ le rayon de l'orbite lunaire (on néglige donc la masse de la lune ici,
et on suppose que son orbite est circulaire). Dans ce cas, l'excentricité s'exprime
(Eq. (A.8) avec $M_0 = M_p$),

$$
e = \sqrt{1 + \left(\frac{b}{a_l}\right)^2 \left(\frac{U'}{v_l}\right)^4}, \quad (A.19)
$$

et le paramètre de la conique est (Eq. (A.7) avec $M_0 = M_p$),

$$
p = \left( \frac{b U'}{a_l v_l} \right)^2, \qquad (A.20)
$$

Le péricentre $q$ est alors donné par (Eq (2.16))

$$
q = a_l (e - 1) / \left( \frac{U'}{v_l} \right)^2 \quad (A.21)
$$

Un impact avec la Lune requiert $q \le a_l$, ce qui entraîne

$$
b_{max} = a_l \sqrt{1 + 2 \left( \frac{v_l}{U'} \right)^2}. \qquad (A.22)
$$

La rotation de la lune engendre une symétrie qui nous permet de faire abstrac-
tion de l'angle d'approche $j'$. On pose $j' = 0$. Les axes $(OX)$ et $(Ox_1)$ sont alors
confondus. Dans $(Ox_1x_2x_3)$, les vecteurs position et vitesse initiaux sont donnés par les équations (A.1) et (A.2). Le vecteur **h** s'exprime d'après l'équation (2.4) comme

$$ \mathbf{h}_{(Ox_1x_2x_3)} = \begin{pmatrix} 0 \\ -U'b \sin \delta \\ U'b \cos \delta \end{pmatrix}, \qquad (A.23) $$

et le vecteur excentricité est donné par l'équation (2.28) comme

$$ \mathbf{e}_{(Ox_1x_2x_3)} = \begin{pmatrix} -1 \\ \frac{b}{a_l} \left(\frac{U'}{v_l}\right)^2 \cos \delta \\ \frac{b}{a_l} \left(\frac{U'}{v_l}\right)^2 \sin \delta \end{pmatrix}, \qquad (A.24) $$

On exprime **h** et **e** en faisant usage de la matrice de passage P :

$$ \mathbf{h}_{(OXYZ)} = \begin{pmatrix} -U'b \cos \delta \sin i' \\ -U'b \sin \delta \\ U'b \cos \delta \cos i' \end{pmatrix}, \qquad (A.25) $$

et

$$ \mathbf{e}_{(OXYZ)} = \begin{pmatrix} -\cos i' - \sin i' \frac{b}{a_l} \left(\frac{U'}{v_l}\right)^2 \sin \delta \\ \frac{b}{a_l} \left(\frac{U'}{v_l}\right)^2 \cos \delta \\ -\sin i' + \cos i' \frac{b}{a_l} \left(\frac{U'}{v_l}\right)^2 \sin \delta \end{pmatrix}. \qquad (A.26) $$

A ce stade, toutes les grandeurs seront exprimées dans (OXYZ). On calcule le vecteur nodal avec l'équation (2.29) comme

$$ \mathbf{n} = \begin{pmatrix} \sin \delta \\ -\cos \delta \sin i' \\ 0 \end{pmatrix}. \qquad (A.27) $$

L'inclinaison par rapport au plan de référence (plan orbital de la planète) est donnée par (Eq. (2.30)),

$$ \cos i = \cos \delta \cos i', \qquad (A.28) $$

et

$$ \sin i = \sqrt{1 - \cos^2 \delta \cos^2 i'}. \qquad (A.29) $$

La longitude du noeud ascendant, d'après les équations (2.31) et (2.32), s'exprime

$$ \cos \Omega = \sin \delta / \sin i, \qquad (A.30) $$
et

$$ \sin \Omega = - \frac{\cos \delta \sin i' / \sin i}{\sin \delta \cos i' / \sin i}. \quad (A.31) $$

L'argument du péricentre se trouve à partir de l'équation (2.33), tel que

$$ \cos \omega = - \frac{\sqrt{e^2 - 1} \sin i' + \sin \delta \cos i'}{e \sin i}, \quad (A.32) $$

où

$$ \omega \in [0, \pi] \quad \text{si} \quad \sqrt{e^2 - 1} \sin \delta \cos i' \geq \sin i' \quad (A.33) $$

$$ \omega \in ]\pi, 2\pi[ \quad \text{si} \quad \sqrt{e^2 - 1} \sin \delta \cos i' < \sin i'. \quad (A.34) $$
# Annexe B

## Etablissement des probabilités d'Öpik

Cette annexe détaille l'obtention de la probabilité d'impact présentée en section 5.2.1 dans le cas où la cible est en orbite circulaire.

<figure><img src="image_81.png" /><figcaption>FIG. B.1 – Géométrie du mouvement autour du noeud commun.</figcaption></figure>

Considérons la géométrie du mouvement au voisinage du noeud commun illustrée par la figure B.1. Les axes (OX) et (OY) du repère orbital (OXYZ) défini par la figure 5.1 sont ici inversés par souci de cohérence avec les articles traitant du sujet. Le projectile et la cible possèdent respectivement une vitesse **v** et **v**₀.
On suppose que, dans le voisinage du noeud commun, le mouvement des deux objets s'assimile à une ligne droite, et que la vitesse du projectile est constante.
Puisque que la collision se fera dans ce voisinage, la vitesse relative d'approche
du projectile vers la cible sera calculée directement à partir des vitesses de l'un et de l'autre évaluées au noeud commun. L'angle que fait **v** avec (-OX) est noté α, tandis que $\mathbf{v}_0$ est selon (-OY) puisque l'orbite est circulaire. On se trouve sur la figure au niveau du noeud ascendant, et la composante de **v** sur (OZ) est positive. La distance $\ell$ est donnée par $\ell = r_{\text{noeud}} - a_0$, et $\ell = 0$ correspond à l'une ou l'autre des quatre valeurs $\omega_n$. A un temps donné, projectile et cible se trouvent à une position arbitraire sur leur orbite, ce que l'on traduit par l'usage de la grandeur $\ell'$ : on décide qu'à $t = 0$, le projectile coupe le plan (OXY), tandis que la cible est éloignée d'une distance $\ell'$ de l'origine.

Les relations géométriques sur l'ellipse permettent d'exprimer $\alpha$ comme

$$ \cos \alpha = \frac{eS}{\sqrt{1 + 2eC + e^2}} \quad \text{et} \quad \sin \alpha = \frac{1 + eC}{\sqrt{1 + 2eC + e^2}}. \qquad (B.1) $$

La vitesse du projectile au noeud s'exprime dans le repère (OXYZ) comme

$$ \left\{ \begin{array}{l} v_X = -v \cos \alpha, \\ v_Y = -v \sin \alpha \cos i, \\ v_Z = v \sin \alpha \sin i. \end{array} \right. \qquad (B.2) $$

avec $v = \sqrt{GM_s/a}\sqrt{\frac{1+e^2+2eC}{1-e^2}}$ d'après les relations sur une orbite keplérienne. La vitesse de la cible, selon -(OY), est donnée par $v_0 = \sqrt{GM_s/a_0}$. Avec C donné par l'équation (5.13) et $S^2 = 1 - C^2$, on exprime la vitesse relative d'approche $\mathbf{U}_\infty = \mathbf{v} - \mathbf{v}_0$ en fonction des éléments orbitaux selon l'équation (5.17).

Sous l'approximation que les mouvements sont rectilignes, la distance minmale $d_{min}$ qui sépare les deux trajectoires s'écrit

$$ d_{min} = |\ell| \frac{\sin i}{\sqrt{\cot^2 \alpha + \sin^2 i}}. \qquad (B.3) $$

On en déduit que $d_{min} \le \tau$ correspond à une valeur maximale de $l$ donnée par

$$ \ell_{max} = \tau \frac{\sqrt{\cot^2 \alpha + \sin^2 i}}{\sin i} = \tau \sqrt{1 + U_X^2/U_Z^2}, \qquad (B.4) $$

où la dernière égalité fait intervenir les relations données par l'équation (B.2).

Exprimons $\frac{\Delta\omega}{2\pi}$. En dérivant $r$ par rapport à $\omega$, on trouve (Eq. (2.13)),

$$ \left(\frac{dr}{d\omega}\right)_{\text{noeud}} = a(1 - e^2) \frac{eS}{(1 + eC)^2} = a_0 \cot \alpha = a_0 \frac{U_X}{U_Z} \frac{1}{\sin i}. \qquad (B.5) $$

La collision est possible autour de $r = a_0$ pour $d_{min} < \tau$, soit $\Delta\ell = \Delta r = 2\ell_{max}$. Avec $\Delta r/\Delta\omega \simeq |(dr/d\omega)_{noeud}|$, on trouve

$$ \frac{\Delta\omega}{2\pi} = \frac{\tau}{a_0 \pi} \frac{1}{\sin i} \frac{\sqrt{U_{\infty}^2 - U_Y^2}}{|U_X|}. \qquad (B.6) $$
Il reste à exprimer $\frac{\Delta t}{T_0}$. Appelons $\Delta_0 t$ l'intervalle de temps pendant lequel la collision est possible quand $\omega = \omega_n$, c'est-à dire quand $\ell = 0$ et que le projectile coupe l'origine du repère (OXYZ) (Fig. B.1). Alors, en moyenne sur toutes les valeurs de $\ell \in [-\ell_{max}, \ell_{max}]$, $\Delta t = \frac{\pi}{4}\Delta_0 t$. Considérons donc le cas particulier $\ell = 0$. À $t = 0$, on décide que le projectile et l'origine du repère coïncident ($X = 0, Y = 0, Z = 0$, soit $\ell = 0$). La cible, elle, est à une distance $\ell'$ de l'origine ($X_0 = 0, Y_0 = -\ell', Z_0 = 0$). Au temps $t = \delta t$, la distance $L$ qui sépare les deux objets s'exprime donc

$$L^2 = (X_0 - U_X \delta t)^2 + (Y_0 - U_Y \delta t)^2 + (Z_0 - U_Z \delta t)^2 \quad (B.7)$$

$$= (U_X^2 + U_Z^2)\delta t^2 + (\ell' - U_Y \delta t)^2 . \quad (B.8)$$

La valeur minimale de $L$ est obtenue pour le temps qui annule la dérivée temporelle de l'équation (B.8), telle que

$$L_{min}^2 = \ell'^2 \left( 1 - \frac{U_y^2}{U_\infty^2} \right). \quad (B.9)$$

La valeur maximale de $\ell'$ qui permet $L_{min} < \tau$ est donc

$$\ell'_{max} = \tau \frac{U_{\infty}}{\sqrt{U_{\infty}^2 - U_y^2}}. \quad (B.10)$$

L'intervalle de temps durant lequel la collision est possible ($|\ell'| < \ell'_{max}$) est alors donné par

$$\Delta_0 t = \frac{2\ell'_{max}}{v_0}. \quad (B.11)$$

Avec $\Delta t = \frac{\pi}{4}\Delta_0 t$ et $T_0 = 2\pi a_0^{3/2}/\sqrt{GM_s}$, on obtient

$$\frac{\Delta t}{T_0} = \frac{\tau}{4a_0} \frac{U_\infty}{\sqrt{U_\infty^2 - U_y^2}}.$$
# Annexe C

## Simulation numérique dans le système Terre-Lune

Aux premiers stades de cette étude, le bombardement au sein du système Terre-Lune était simulé de manière purement numérique. La méthode semi-analytique, développée aux chapitres 5 et 6 notamment, permet d'obtenir les mêmes résultats en un temps extrêmement plus réduit en se basant sur certaines approximations : une dizaine de minutes contre une ou deux semaines sur un ordinateur récent. Nous présentons néanmoins la méthode numérique ici, essentiellement comme outil de vérification.

### C.1 Principe

Les conditions d'approche avec l'ensemble du système Terre-Lune sont calculées à partir des probabilités d'Öpik décrites au chapitre 5. La cible est ici la sphère de Hill terrestre, de rayon $R_H = 1.5 \times 10^6$ km et de masse nulle (par définition, la masse terrestre n'influe sur le projectile que lorsqu'il a pénétré dans cette sphère). On définit le repère (OXYZ), centré sur la Terre, tel que (OXY) définisse le plan de l'écliptique. La rotation de la Lune autour de la Terre rend le problème insensible à la direction d'approche dans le plan de l'écliptique. On choisit (OX) comme direction de référence. L'un après l'autre, des projectiles sont choisis au hasard et leur trajectoire intégrée numériquement jusqu'à l'impact ou l'éjection hors du système. Pour ce faire, on tire au hasard un jeu de conditions d'approche ($U_\infty, i_\infty$) selon leur distribution de probabilité. On rappelle que $U_\infty$ est la vitesse relative d'approche avec le système, et $i_\infty$ l'inclinaison de $U_\infty$ par rapport à l'écliptique. La vitesse initiale est alors donnée par $\dot{x}_0 = 0$, $\dot{y}_0 = U_\infty \cos i_\infty$ et $\dot{z}_0 = -U_\infty \sin i_\infty$. Le projectile est susceptible de pénétrer dans le système en n'importe quel point de la sphère de Hill : on choisit $x_0, y_0$ et $z_0$ au hasard, tels que l'ensemble des projectiles soient initialement distribués
uniformément sur l'hémisphère de Hill définit par la direction (OX). La position initiale de la Lune sur son orbite est elle aussi choisie au hasard à chaque tir.

La trajectoire du projectile, soumise à la fois à l'attraction terrestre et lunaire, vérifie

$$ \ddot{x} = -GM_T \cdot \frac{x}{(x^2 + y^2 + z^2)^{\frac{3}{2}}} - GM_l \cdot \frac{x - x_l}{[(x - x_l)^2 + (y - y_l)^2 + (z - z_l)^2]^{\frac{3}{2}}}, $$

$$ \ddot{y} = -GM_T \cdot \frac{y}{(x^2 + y^2 + z^2)^{\frac{3}{2}}} - GM_l \cdot \frac{y - y_l}{[(x - x_l)^2 + (y - y_l)^2 + (z - z_l)^2]^{\frac{3}{2}}}, $$

$$ \ddot{z} = -GM_T \cdot \frac{z}{(x^2 + y^2 + z^2)^{\frac{3}{2}}} - GM_l \cdot \frac{z - z_l}{[(x - x_l)^2 + (y - y_l)^2 + (z - z_l)^2]^{\frac{3}{2}}}, $$

où les indices T et l font respectivement référence à la Terre et à la Lune, et M désigne la masse. Les coordonnées de la Lune sont données par une rotation uniforme (orbite circulaire) dans son plan orbital, de vitesse

$$ v_l = \sqrt{\frac{G(M_T + M_l)}{a_l}} . $$

Notons que si on prend en compte l'inclinaison de l'orbite lunaire, il faut choisir à chaque tir une valeur de $\Omega_l$, l'argument du noeud ascendant, uniformément compris entre 0 et $2\pi$.

La trajectoire d'un projectile est calculée en résolvant numériquement la relation fondamentale de la dynamique à l'aide de l'algorithme de Bulirsch-Stoer [Press et al., 1992]. Cet algorithme est un schéma de différences finies qui adapte le pas de calcul (ici le pas de temps) selon les besoins. Ainsi, plus d'itérations sont effectuées dans les zones à fort gradient de vitesse, c'est-à-dire à proximité de la Terre ou de la Lune. Ce type de procédé permet d'économiser un précieux temps de calcul sans pour autant sacrifier à la précision. L'idée clef est la suivante : le calcul effectué dans l'intervalle $\Delta t$ est considéré comme une fonction polynomiale de variable $\Delta t^2$. Un pas de temps initial relativement grand est subdivisé de plus en plus finement, et la convergence des calculs est extrapolée pour donner le résultat à $t + \Delta t$ comme s'il avait été calculé par une infinité de sous-pas. La figure C.1 illustre l'adaptation du pas de temps selon les besoins. La distance entre les points qui définissent la position du projectile correspond au pas de temps choisi par l'algorithme. On peut voir que plus de points sont évalués à proximité de la Lune. On constate également que, dans certains cas extrêmes, la trajectoire du projectile diverge complètement de l'hyperbole du problème à deux corps.
<figure><img src="image_82.png" /><figcaption>FIG. C.1 – Trajectoire numérique du projectile.</figcaption></figure>

## C.2 Détermination de l'impact.

Puisque l'algorithme adapte lui même le pas de temps, et que la Lune est en mouvement dans l'intervalle, il est nécessaire de prendre certaines précautions afin de ne manquer aucun impact. La méthode mise en oeuvre est la suivante :
Considérons $\Delta r(t)$ tel que

$$ \Delta r(t) = \| \mathbf{r}(t) - \mathbf{r}_l(t) \| - R_l, $$

où $\mathbf{r}$ est la position du projectile, $\mathbf{r}_l$ la position de la Lune, $R_l$ son rayon.

1. On tire un projectile au hasard. Sa trajectoire est calculée jusqu'à ce que :

(i) $\Delta r$ augmente après avoir décrût,
soit $\Delta r(t_i) > \Delta r(t_{i-1})$ et $\Delta r(t_{i-1}) < \Delta r(t_{i-2})$.
Aller à 2.
ou

(ii) $\Delta r(t_i) < 0$.

Aller à 3.

2. Dans le cas (i) seulement, les trois valeurs de $\Delta r$ sont interpolées par un polynôme de degré 3, et le minimum $\Delta r_{min}$ est obtenu en annulant la dérivée $\frac{\partial \Delta r}{\partial t}$. Si $\Delta r_{min} > D_l$, avec $D_l = 4R_l$, alors le calcul de la trajectoire continue normalement jusqu'à ce que les cas (i) ou (ii) se reproduisent, ou bien jusqu'à ce que le projectile quitte définitivement le système Terre-Lune.
Sinon, aller à 3.

3. Le calcul reprend à $t_{i-2}$ dans le cas (i), ou à $t_{i-1}$ dans le cas (ii). Cette fois, $D_l = D_l/2$ et l'on impose à l'algorithme un pas constant $\Delta t_{imp} = \min(t_i - t_{i-1}, t_{i-1} - t_{i-2})/2$ (ou $\Delta t_{imp} = \Delta t_{imp}/2$ si la boucle a déjà été effectuée).
Aller à 1.

La boucle s'effectue tant que $\Delta r_{min} < D_l$, pour $D_l$ et $\Delta t$ de plus en plus petits, comme illustré par la figure C.2. Si cette condition vient à ne plus être respectée, le calcul de la trajectoire se poursuit et l'algorithme est de nouveau libre d'adapter le pas de temps à sa guise. L'approximation polynomiale se justifie par le fait qu'à proximité de la planète ou du satellite, le mouvement s'apparente à la conique d'un problème à deux corps. On peut voir sur la figure C.2 la validité de ce choix.

On considère qu'il y a impact lorsque $\Delta r(t_i) < 0$ avec $||\dot{\mathbf{r}}||\Delta t_{imp} < 30$ km, distance qui correspond à 1 degré sur la surface de la Lune. En cas de collision, le moment $t_{imp}$ de l'impact est interpolé pour correspondre à $\Delta r(t) = 0$, en supposant une variation linéaire entre $t$ et $t - \Delta t_{imp}$. Toujours de manière linéaire, les vecteurs position et vitesse du projectile et de la Lune sont déterminés à $t_{imp}$ dans le référentiel (OXYZ). On se place alors dans le référentiel géographique lunaire (en tenant compte d'une éventuelle obliquité) pour déduire les coordonnées, la vitesse et l'angle d'impact. La même méthode est appliquée à la Terre, à ceci près qu'elle est fixée dans le référentiel choisi. Le choix de la valeur de $D_l$ est fait empiriquement : pour une valeur supérieure à $4R_l$, on ne recense pas plus d'impacts. On lance un nouveau projectile lorsque le précédent à quitté le système, soit $r > R_H$.
<figure><img src="image_83.png" /><figcaption>FIG. C.2 – Détermination d'un impact dans les simulations numériques. Illustration de l'algorithme décrit dans le texte.</figcaption></figure>

## C.3 Précision et optimisation du code

La précision du code est ajustée de la manière suivante. On calcule avec l'algorithme la position du projectile à la sortie du système Terre-Lune, pour une Lune de masse nulle. Puis on compare, au même temps, cette position avec celle déterminée analytiquement dans le cadre d'un problème à deux corps (Terre et projectile).

Dans la mesure où l'équation de la conique ne fait pas intervenir le temps expliquement, il est nécessaire d'avoir recours à une méthode itérative pour connaître la position du projectile sur l'hyperbole du problème à deux corps à un temps précis. L'objectif est de trouver la valeur $u_0$ qui annule la fonction

$$F(u) = \sigma_0 u^2 c(\alpha u^2) + (1 - \alpha r_0) u^3 s(\alpha u^2) + r_0 u - \sqrt{GM_T}(t - t_0),$$

où

$$\sigma_0 = \frac{\mathbf{r}_0 \cdot \dot{\mathbf{r}}_0}{\sqrt{GM_T}},$$

$$\alpha = \frac{2}{r_0} - \frac{\dot{\mathbf{r}}_0^2}{\sqrt{GM_T}},$$
et

$$c(u') = \frac{\cosh(\sqrt{-u'}) - 1}{-u'}$$ $$s(u') = \frac{\sinh(\sqrt{-u'}) - \sqrt{-u'}}{\sqrt{-u'^3}}$$

Ayant déterminé $u_0$ par itération, par exemple à l'aide de l'algorithme de Newton [Press et al., 1992], la position à l'instant $t$ est donnée (voir par exemple Murray et Dermott [2000]) par

$$\mathbf{r}(t) = \left[ 1 - \frac{u_0^2}{r_0} c(\alpha u_0^2) \right] \mathbf{r}_0 + \left[ (t - t_0) - \frac{u_0^3}{\sqrt{GM_T}} s(\alpha u_0^2) \right] \dot{\mathbf{r}}_0 .$$

Le paramètre de précision de l'algorithme est choisi tel que l'écart entre la trajectoire numérique (en l'absence de Lune) et la trajectoire analytique n'excède pas la dizaine de mètres dans les cas les plus sensibles (qui correspondent à des trajectoires qui frôlent la Terre, où le projectile subit un fort gradient de vitesse). Il est vérifié qu'une fois la Lune incluse dans les simulations, une précision accrue n'induit pas une différence significative dans la localisation des impacts. Au final, sachant que la position de l'impact est interpolée sur un intervalle de temps au cours duquel la distance parcourue par le projectile n'excède pas 30 km, on peut considérer que l'erreur sur la surface lunaire est inférieure au degré.

## C.4 Comparaison avec la méthode semi-analytique.

A la différence de la méthode semi-analytique construite aux chapitres 5 et 6, l'intégration purement numérique des trajectoires dans le système Terre-Lune possède l'intérêt d'éviter un grand nombre d'approximations. On ne calcule plus ici les conditions d'approche lunaires « à l'infini », mais directement les conditions de l'impact, sans qu'il soit nécessaire de diviser artificiellement l'espace au sein du système Terre-Lune en une zone d'influence purement terrestre et une zone d'influence purement lunaire.

La figure C.3 compare les deux méthodes. Elle représente, en chaque point de la surface lunaire, la différence entre le flux d'impacts calculé au chapitre 6 (figure 6.14) et les résultats obtenus numériquement. La variation des résidus, visiblement décorrélée de la position géographique, n'est due qu'au bruit statis- tique de la simulation numérique (qui se compose de deux millions d'impacts, et a nécessité une semaine de calculs). On note que la simulation numérique prend en compte l'inclinaison du plan orbital lunaire par rapport à l'écliptique (5°), à la différence de la méthode semi-analytique, sans que cela n'ait d'incidence sur
la validité des résultats de cette dernière. Des résultats similaires sont obtenus
pour tous les tests que nous avons pu effectuer. Notamment pour n'importe quel
type de distance Terre-Lune comprise entre 10 et 60 rayons terrestres. Afin de
s'assurer que d'éventuelles différences systématiques n'étaient pas masquées par
la diversité des conditions d'approche susceptibles de se contrebalancer, le bom-
bardement généré par un seul couple ($U_\infty, j_\infty$) a également été testé avec succès,
pour des couples très variés.

<figure><img src="image_84.png" /><figcaption>FIG. C.3 – Résidus entre la méthode analytique et la simulation numérique, pour les résultats présentés en figure 6.14.</figcaption></figure>

La comparaison très favorable permet d'utiliser les méthodes semi-analytiques
écrites pour le système planète-lune, extrêmement plus rapides, avec confiance.
Cette confiance peut également s'étendre aux conditions d'approche calculées
pour les planètes par les probabilités d'Öpik, puisque la géométrie du problème
est la même, et que seule l'échelle change (le satellite est remplacé par la planète,
la planète par le Soleil). Cela justifie le fait de ne pas avoir intégré numériquement
la trajectoire des projectiles dans le système Soleil-planète.
# Annexe D

**Article :** Nonuniform cratering of the terrestrial planets.
<figure><img src="image_85.png" /></figure>

# Nonuniform cratering of the terrestrial planets

Mathieu Le Feuvre *, Mark A. Wieczorek

Institut de Physique du Globe de Paris, 4, avenue de Neptune, 94107 Saint-Maur, France

## ARTICLE INFO

**Article history:**

Received 26 July 2007

Revised 28 March 2008

**Keywords:**

Cratering

Terrestrial planets

## ABSTRACT

We estimate the impact flux and cratering rate as a function of latitude on the terrestrial planets using a model distribution of planet crossing asteroids and comets [Bottke, W.F., Morbidelli, A., Jedicke, R., Petit, J.-M., Levison, H.F., Michel, P., Metcalfe, T.S., 2002. Icarus 156, 399–433]. After determining the planetary impact probabilities as a function of the relative encounter velocity and encounter inclination, the impact positions are calculated analytically, assuming the projectiles follow hyperbolic paths during the encounter phase. As the source of projectiles is not isotropic, latitudinal variations of the impact flux are predicted: the calculated ratio between the pole and equator is 1.05 for Mercury, 1.00 for Venus, 0.96 for the Earth, 0.90 for the Moon, and 1.14 for Mars over its long-term obliquity variation history. By taking into account the latitudinal dependence of the impact velocity and impact angle, and by using a crater scaling law that depends on the vertical component of the impact velocity, the latitudinal variations of the cratering rate (the number of craters with a given size formed per unit time and unit area) is in general enhanced. With respect to the equator, the polar cratering rate is about 30% larger on Mars and 10% on Mercury, whereas it is 10% less on the Earth and 20% less on the Moon. The cratering rate is found to be uniform on Venus. The relative global impact fluxes on Mercury, Venus, the Earth and Mars are calculated with respect to the Moon, and we find values of 1.9, 1.8, 1.6, and 2.8, respectively. Our results show that the relative shape of the crater size-frequency distribution does not noticeably depend upon latitude for any of the terrestrial bodies in this study. Nevertheless, by neglecting the expected latitudinal variations of the cratering rate, systematic errors of 20–30% in the age of planetary surfaces could exist between equatorial and polar regions when using the crater chronology method.

© 2008 Elsevier Inc. All rights reserved.

## 1. Introduction

Impact cratering is not only a major geologic process that shapes planetary surfaces (see Melosh, 1989), but it is also a method by which surfaces can be dated. In essence, a fresh surface will accumulate craters at a given rate, and the number of craters per square kilometer can be converted to an approximate age (e.g., Neukum et al., 2001; Stöffler and Ryder, 2001; Hartmann and Neukum, 2001). While certain regions of the Moon have been dated directly by the analysis of radiogenic isotopes in samples with known geologic context, for other planets, such as Mercury, Venus, and Mars, “crater counting” is one of the only reliable means for placing relative geologic stratigraphy into an approximate absolute framework.

When dealing with planetary cratering rates, and in particular when using this crater chronology method, it is commonly assumed that a given bolide is equally likely to impact at any location on the planet's surface. Furthermore, most studies also assume that the impact velocity and impact angle distributions

are independent of position (e.g., Pierazzo and Melosh, 2000; Neukum et al., 2001; Strom et al., 2005). Given that the impact flux is assumed to be everywhere the same, and that the impact velocity and angle are also assumed not to depend upon position, the resulting cratering rate would also be independent of position.

The above assumptions would be expected to be true if the velocity distribution of candidate impactors were isotropic (e.g., Pierazzo and Melosh, 2000), or if the planet were sufficiently massive so that the trajectories of objects with high encounter velocities could be gravitationally deflected. However, the assumption of isotropy is questionable. Furthermore, many impactors have high enough approach velocities such that their trajectories are not significantly deviated by the planet. Indeed, in the case of zero inclination with respect to the equator of the planet, simple geometric considerations predict that objects with infinite approach velocities will produce an impact flux that varies with the cosine of the geocentric latitude. Conversely, if the average relative inclination of projectiles striking an object were extremely high, one would expect more impacts to occur (per unit area) at the poles.

With few exceptions (e.g., Morota and Furumoto, 2003; Morota et al., 2005; Zahnle et al., 2001) spatial variations of the cratering rate have not been considered when applying the crater chronology method. Assuming that the average cratering rate of a planet

<footnote>* Corresponding author. Fax: +33 (0) 1 45 11 42 57.
E-mail address: lefeuvre@ipgp.jussieu.fr (M. Le Feuvre).</footnote>