<figure><img src="image_1.png" /></figure>

HAL
open science

# Predicting Deep Moonquake Source Regions Using Their Temporal and Spatial Patterns and Machine Learning

Josipa Majstorović, Philippe Lognonné, Taichi Kawamura, Mark Panning

<figure><img src="image_2.png" /></figure>

## To cite this version:

Josipa Majstorović, Philippe Lognonné, Taichi Kawamura, Mark Panning. Predicting Deep Moonquake Source Regions Using Their Temporal and Spatial Patterns and Machine Learning. Journal of Geophysical Research. Machine Learning and Computation, 2024, 1 (2), 10.1029/2023JH000117. hal-05070626

HAL Id: hal-05070626

https://hal.science/hal-05070626v1

Submitted on 16 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2023JH000117

**Key Points:**

* We study the seismic data acquired during the Apollo missions for the future single-station deployment on the far side of the Moon

* We propose the classification model of deep moonquake source regions using only their occurrence times

* The models perform with accuracy greater than 70% when trained to classify combinations of four or fewer source regions

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

J. Majstorović,
josipa.majstorovic@protonmail.com

**Citation:**

Majstorović, J., Lognonné, P., Kawamura, T., & Panning, M. P. (2024). Predicting deep moonquake source regions using their temporal and spatial patterns and machine learning. *Journal of Geophysical Research: Machine Learning and Computation, 1*, e2023JH000117. https://doi.org/10.1029/2023JH000117

Received 10 JAN 2024
Accepted 5 JUN 2024

© 2024 The Author(s). Journal of Geophysical Research: Machine Learning and Computation published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_4.png" /></figure>

# Predicting Deep Moonquake Source Regions Using Their Temporal and Spatial Patterns and Machine Learning

Josipa Majstorović¹ ID, Philippe Lognonné¹ ID, Taichi Kawamura¹ ID, and Mark P. Panning² ID

¹Institut de physique du globe de Paris, Université Paris Cité, CNRS, Paris, France, ²Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA

**Abstract** In the near future, lunar exploration should be enhanced by deploying a new seismic station on the farside of the Moon, the Farside Seismic Suite (FSS), a package of two seismometers recently selected by NASA to fly on a commercial lander. The new data should provide us with new insights into Moon's seismic activity, its interior, and composition. To fully benefit from the new data, we need to take advantage of the data acquired during the Apollo missions. The problem of relating new and old data is complex due to the single-station nature of the future deployment. In this study, we tackle this issue by developing a machine learning model in the context of the deep moonquake (DMQ) classification problem. The DMQs form the largest group of detected events from Apollo data, and their source regions have been located and are known to exhibit temporal and spatial patterns connected with the monthly lunar tidal periods. Therefore, we propose to utilize a machine learning (ML) algorithm named random forest to identify DMQs source regions without using waveform information and only using the lunar orbital parameters related to DMQs time occurrences. We show that ML models perform well (with an accuracy >70%) when they are trained to classify 4 or fewer different source regions. This approach gives us a good first location approximation of the DMQs source regions and opens up a new approach to their location estimate when captured by the future FSS single-station seismometers.

**Plain Language Summary** Future space missions will provide us with new in situ ground vibration measurements of the Moon. From the Apollo missions in 1970, we know that the Moon can produce various seismic events. The most numerous ones are deep moonquakes (DMQs), which are events associated with the displacement deep within the lunar interior. These events occur in specific source regions, and their occurrence is related to the monthly motion of the Moon around the Earth. To further study lunar interior and DMQs with future missions, we need to define a way to relate the new data with the existing Apollo data. Due to the single-station nature of the future deployment, the mapping between these two data sets is challenging with classical techniques. Therefore, in this study, we tackle this issue using machine learning (ML). We propose a model that is able to classify DMQs source regions using only the time of the occurrence of quakes. We show that models perform well (with an accuracy >70%) when they are trained to classify 4 or fewer source regions. This gives us a good preliminary location of the DMQs source regions and an opportunity to detect them with future lunar missions.

## 1. Introduction

We are at the beginning of a golden age of lunar exploration as many nations, together with private companies, are establishing numerous efforts to obtain new scientific measurement of the Moon (Kawamura et al., 2022; Lognonné, Schmerr, et al., 2020; Pickrell, 2022; Wang et al., 2023; Weber et al., 2021). In light of this, NASA established the Artemis program, which should land a crewed mission at the lunar south pole (Witze, 2022) and deploy the lunar environment monitoring station (Benna et al., 2020). This would be the first attempt of a crewed landing after the successful Apollo missions in 1970's. Before Artemis missions land on the Moon, NASA has also established the Commercial Lunar Payload Services (CLPS) program to land scientific payloads on the Moon using commercial landers. The Farside Seismic Suite (FSS) is one of the selected payloads, and it will deliver two co-located seismometers to Schrödinger basin on the farside of the Moon (Cutler et al., 2023; Panning et al., 2021; Standley et al., 2023): one vertical very broadband seismometer and a three-axes short period sensor, both spared or derived from the SEIS experiment sensors (Lognonné, Banerdt, et al., 2020; Lognonné et al., 2019) from the InSight mission to Mars (Banerdt et al., 2020). Due to the unique opportunity and diverse application of the future ground motion sensors' deployment, the onboard instruments are capable of measuring wide range of frequencies,
from 0.01 to 50 Hz, although mission requirements focus between 0.1 and 1 Hz, where we have good knowledge of expected signals from the Apollo seismometers.

Apollo missions showed the importance of deploying sensors on the surface of the Moon, since a great amount of our knowledge about the Moon comes from the analysis of data acquired during the Apollo missions (Lognonné & Johnson, 2015). Thus, analyzing ground motion measurements provided the community with the first constraints on the lunar interior and the activity at the surface (Gagnepain-Beyneix et al., 2006; Garcia et al., 2011, 2019; Kawamura et al., 2017; Khan & Mosegaard, 2002; Khan et al., 2000, 2014; Lognonné et al., 2003; Nakamura et al., 1982; Nunn et al., 2020a; Weber et al., 2010). It has also revealed that the Moon can host events of various origins, such as shallow and deep moonquakes, meteoroids, and artificial impacts (Dainty et al., 1975; Lammlein, 1977; Nakamura, 1983, 2003, 2005; Toksöz et al., 1974). Today, we have more than twelve thousand events, out of which the deep moonquakes (DMQs) form the most numerous group with more than 7,000 events reported (Nakamura, 2005; Nakamura et al., 1981).

DMQs are a distinctive group of seismic events that originate from depths between 700 and 1,200 km at high pressure and temperature conditions, where little brittle deformation is expected. Due to a very high waveform similarity between quakes for a fixed pair of station and source region, the DMQs have been clustered into about 300 source regions or nests (Frohlich & Nakamura, 2009; Nakamura, 2003). This has interpreted to be a consequence of DMQs occurring repeatedly at the fixed nests, which are located mostly on the near side of the Moon. It has been shown that time occurrence of the DMQs is correlated with the monthly motion of the Moon around the Earth. Thus, DMQ occurrences exhibit tidal periodicities, and furthermore, the associated high strain rates might explain brittle processes (Kawamura et al., 2017). However, the real causes of their origins are yet to be discovered. There are still significant puzzles about their origin, as cyclic tidal stresses, caused by the monthly motion of Moon around Earth, are far less than the estimated confining pressures where DMQs occur (Cheng & Toksöz, 1978; Minshull & Goulty, 1988) suggesting extremely low cohesion and friction values if a Mohr-Coulomb failure model is used (Pou et al., 2024). This opens the question of whether we need both tectonic and ambient tidal stresses to explain their mechanical origin (Frohlich & Nakamura, 2009). Of course, there is no guarantee that the DMQ nests observed over several years during the Apollo mission are still active nearly 50 years later, but considering that tidal attraction between the Earth and the Moon is still present and relatively consistent over time and that the Apollo data showed nests activated consistently over tens of tidal cycles, it is reasonable to suggest that some may still be active today.

To better constrain the lunar interior and unravel the cause of DMQs, it is important to locate new events and associate them with the Apollo known nests. We are confident that next lunar missions shall capture these new events. However, due to the complexity of space missions, it is highly likely that at the beginning of the lunar exploration era we might have only one ground sensor at a single station such as FSS. Therefore, in this paper, we study the events' classification from particular Apollo nests designed for single-station configuration. We tackle this problem without using the events waveforms recorded during the Apollo mission. This is because the location of the new ground sensor will not match the existing Apollo-era stations, and therefore, the waveform templates from the known nests locations shall be different due to different epicentral distances and propagation paths as well as differences in instrument response for modern instruments compared with the Apollo instruments. We propose a machine learning (ML) model that is trained to classify nests within the set of nests of similar differential travel times. The main features used for the model training are related to the fact that different nests respond differently to lunar cycle. Each new single station on the Moon will indeed add a new differential S-P wave travel time and possibly, due to the modern broadband capability, other phases. Then, these additional data can be used in addition to the previously determined differential travel times to update spherically symmetric models of the Moon (with updated location of the DMQs). The differential travel times of the events occurring at epicentral distances much larger than those sampled during the Apollo era can also be used, if their position is somehow constrained even in probabilistic way. The differential travel times from those events should provide key constraints on seismic velocities in the very deep lunar interior and/or in areas far from the Apollo network.

Very early studies have shown correlation between lunar transient events and position of the Moon relative to the Earth (Cameron & Gilheany, 1967; Middlehurst, 1967; Moore, 1968). This motivated the hypothesis that some moonquakes occur with periods that reflect Earth-Moon-Sun relationship (Ewing et al., 1971). Later, it was shown that the occurrence of DMQs is related to tidal stress cycles, and correlations between DMQs occurrence times and lunar monthly tidal cycles have been indicated (Cheng & Toksöz, 1978; Lammlein, 1977; Lammlein
et al., 1974; Minshull & Goulty, 1988; Toksöz et al., 1977). The lunar cycle can be explained with three lunar months: synodic, draconic, and anomalistic. Synodic month is the period of lunar phases such as New Moon, First Quarter, Full Moon, and Last Quarter. Draconic month is the period between two nodes, ascending or descending, where the nodes are points at which the Moon's orbit plane crosses the ecliptic plane toward which it is inclined of about 5.14°. Thus, the descending node is where Moon transfers from the northern hemisphere to the southern one, and ascending node is the opposite motion. Anomalistic month is the period between two extreme points, perigee and apogee, since the Moon's orbit approximates an ellipse rather than a circle. Earlier studies counted the number of events per day as a function of time and found 0.5 and 1 month signals in the occurrence times related to anomalistic and draconic period of 27 days (Lammlein, 1977; Lammlein et al., 1974). The same studies also indicated 206-day and 6-year periods related to the Sun's perturbation on the lunar orbit and the relative pre- cession of the perigee of the Moon's orbit. Subsequently, many recent papers studied and confirmed the tidal periodicities of DMQs and more (Bills et al., 2008; Bulow et al., 2005, 2007; Frohlich & Nakamura, 2009; Turner et al., 2022; Weber et al., 2009, 2010). Over the years, most studies focused on the most numerous nests to report their occurrence periods. Therefore, it has been reported that the largest nest, A1, exhibits short-term periodicity associated with draconic, anomalistic, and synodic months but also long-term periodicity of 6 years related to beat between the draconic and anomalistic months (Bulow et al., 2007; Frohlich & Nakamura, 2009). Bulow et al. (2007) reported that the nests A8 and A33 exhibit only anomalistic periodicity, while A9 and A20 nests exhibit only draconic one. In the same study, the nests A6, A10 (with new label A24), A14, and A18 show the impact of both draconic and anomalistic months, with a dominant peak associated with the anomalistic month. Next, Weber et al. (2010) applied the delay-coordinate DC map technique to show how most nests (92% out of total analyzed), that have more than 35 events, exhibit some type of periodic behavior. According to this tech- nique, the DMQs exhibit five different DC map patterns: triangle, four-corner, one-corner, four-corner (special case), and no pattern. The nests A8 and A33, that Bulow et al. (2007) previously labeled to exhibit only anomalistic period, have been sorted in two separate groups. The A8 nest manifests a triangular DC map pattern, which implies that this nest's activity may be explained with two periods. Other nests that belong to this group alongside the A8 nest are A1, A5, A6, A7, A24 (previously A10), A13, A14, A18, A25, A44, A51, A73, A97, A218, and A224. Next, the nest A33, also previously labeled to exhibit only anomalistic period, manifests one- corner DC map. The one-corner DC map might indicate existence of a single group with a single monthly periodicity. Others nests that belong to this group are A3, A19, A32, A33, A34, A37, A40, A41, A42, A96, and A223. Likewise, Bulow et al. (2007) has associated nests A9 and A20 only with the draconic period, while according to the DC map technique the periodicity of these nests can be associated with four-corners and four- corners special case maps, respectively. The four-corners map indicates a single periodicity, meaning that all events have the same period. Alongside the nest A9, there are few other nests, such as A15, A16, A17, A21, A22, A28, A30, A36, A38, A65, A86, A201, A202, A204, A234, A238, and A257, that exhibit four-corner DC map patterns. The four-corners special case map, associated with the nest A20 and also the nest A39, is related to the fact that events from these nests have different polarities, and isolating events of certain polarity exhibit four- corner DC map patterns. The mentioned association between nests does not imply that nests belonging to the same group have absolutely similar temporal patterns in terms of frequency and released amplitude.

Based on the previous papers, it is clear that DMQ nests exhibit temporal patterns in their occurrences and that these are correlated with the Moon-Earth system. Therefore, the open question is whether we can design features, which would reflect these temporal patterns and to further use those features to study the nest classification with one lunar station. In this paper, we tackle the question of defining optimal features and the classification model. We show that the designed features, when used as an input to the models trained to classify nests, might reflect more profound physical meaning behind the nests' activities, such as their fault plane orientation and geographical location. The paper is organized as follows: first, we discuss the existing catalog of DMQ events. Second, we discuss the feature design. Third, we introduce a classification model. Forth, we discuss successes and pitfalls of this model for nest classification applied to different combination of nests. We conclude how this study can offer some first estimates of the nest locations in the future lunar missions.

## 2. Defining Database for Nests Classification Using Lunar Catalog and Orbital Features

We start with the existing catalog of lunar events (Nakamura et al., 1981), which was updated in 2008 and modified in 2018 (Nunn et al., 2020a). The catalog contains a list of events (shallow and deep moonquakes,
meteoroids, and artificial impacts) with attributes such as date and time of the event occurrence, signal envelope amplitude as measured in mm on a standard plot, data availability per station, and the nest (source) classification for DMQs. It is important to note that the source classification is not an exact location defined by latitude and longitude rather a result from the waveform cross-correlation and single-link cluster analysis (Nakamura, 2003). This analysis positively clustered approximately 7,000 DMQs into about 300 nests (Frohlich & Nakamura, 2009), where the largest nest is associated with label A1, as shown in Figure S1 of the Supporting Information S1. This nest contains 441 quakes, and it is placed on the near side of the Moon.

## 2.1. Catalog Processing: Nest Sets Based on Travel Time Information

Earlier studies published lunar interior models and location of DMQ nests in terms of latitude, longitude, and depth by picking P and S travel times on the quake waveforms (see Garcia et al. (2019) for a review of these picks and Nunn et al. (2020a) for data and instrument overview). Using these lunar interior models and calculated nest locations, we can define nests that are close by in distance if we consider only the $t_{sp} = t_s - t_p$ travel time measurements, where $t_p$ and $t_s$ stand for travel time of P and S waves, respectively. To do so, we assume that (a) our single lunar station is located on the farside in Schrödinger basin with FSS landing site at 71.378°S, 138.248°E (Panning et al., 2021), (b) nests' latitudes and longitudes are fixed and are defined by Lognonné et al. (2003), and (c) calculated P- and S- wave travel times ($t_p$ and $t_s$, respectively) are obtained using a lunar velocity model between the landing site and the locations of DMQ nests. By fixing the nest locations and landing site, we are examining the typical error due to lunar models. We further consider only 23 nests, as analyzed by Lognonné et al. (2003). By having $t_s$ and $t_p$, we can calculate $t_{sp}$ for all nests and existing lunar models to assess the variability of $t_{sp}$ travel times. The considered lunar models are from papers Garcia et al. (2019) (three models labeled as ISSI 1, ISSI 2, and ISSI 3), Garcia et al. (2011), Khan et al. (2014), Matsumoto et al. (2015), and Weber et al. (2011). While their overall comparison is presented by Garcia et al. (2019). The largest discrepancies between these models exist approximately below 1,200 km due to the lack of data to constrain the deepest layer of the Moon. It is worth repeating that all models used some subset of the existing lunar data to invert for the velocity structure, and differences may arise from different travel time readings, inversion techniques, and parametrization of the models. In our approach, we first calculate epicentral distances between nests and station and then travel time of P- and S- seismic waves for all seven lunar models using Python (Obspy package with "taup" module) and Java packages, both algorithms based on paper by Crotwell et al. (1999). Since, we find small difference in the travel time estimates of two algorithms, we further use the average values of $t_p$ and $t_s$ times (see Text S1 and Figure S2 in Supporting Information S1) that are eventually used to calculate the final $t_{sp}$ values for all nests and all seven models (see Figure S3 in Supporting Information S1). Some combination of nests and lunar models don't have estimates of $t_{sp}$ values (model from Weber et al. (2011)), or these estimates are overestimated or underestimated compared to the average values (for ISSI 1 model from Garcia et al. (2019) and model from Matsumoto et al. (2015)). The reason for this is unresolvable $t_p$ and/or $t_s$ values for either Python or Java package (see Figure S4 in Supporting Information S1), which relates to the nests with higher values of epicentral values (see Figure S5 in Supporting Information S1). This might reflect the complexity these three models reflect in the deeper velocity structure, since P and S waves sample a deeper layer with larger epicentral distances. Finally, we decided to work with only four lunar models, which are Garcia et al. (2011) labeled as G11, two models from Garcia et al. (2019) labeled as issi2 and issi3, and model from Khan et al. (2014) labeled as K14.

Next, we classify nests based on these travel time measurement by having $t_{sp}$ as a distance measurement. During our classification, to take into account error that might rise from the picking process itself, uncertainty due to conflicting lunar models, uncertainty on the nests' locations, and lateral variations, we assume a threshold error on $t_{sp}$ measurement of 5 s. This assumption is well-grounded due to several reasons. First, the average picking error in Lognonné et al. (2003) is 5 s. Second, if from the same study we consider picking errors for A1 nest on all four Apollo stations, we obtain the average value of 1.4 s. Third, if for A1 nest we consider different location estimates from studies Lognonné et al. (2003), Nakamura (2005), and Gagnepain-Beyneix et al. (2006) and four lunar models from previous paragraph (G11, issi2, issi3, K14), and for each pair we calculate $t_{sp}$ values, eventually we obtain the standard deviation of 1.04 s. Fourth, if we consider the possible existence of lateral variations at the depth of DMQs, based on the study by Zhao et al. (2008), we can assume lateral variation of 2%, which translates to 3.2–4.8 s travel time error for $t_{sp}$ range from 160 to 240 s. Therefore adding all these errors together into one error budget, we obtain $\sqrt{1.4^2 + 1.04^2 + 4.8^2} = 5.1$ seconds. The future very broadband seismometer should have in general better sensitivity than the Apollo sesimometers (see Figure 5 in Lognonné et al. (2019)), therefore
having an improvement by, at least, a factor of 10 at 1 Hz (Lognonné, Banerdt, et al., 2020; Lognonné et al., 2019). If with the current data set acquired with the Apollo seismometers, we can arrive at the lower limit of the threshold error of 5 s, it seems like a valid choice to keep this threshold for data that shall be acquired with the future improved seismometers. The classification is achieved by counting how many nests have similar travel times not exceeding the 5-s error rule. Therefore, if each nest *i* is associated with travel time measurement *t*<sub>sp,i</sub>, we can form intervals $\delta t_i = [t_{sp,i} - 5, t_{sp,i} + 5]$, $i = 1, ..., 23$. Following, for each nest *i* we count nests *j* that satisfy $t_{sp,i} \neq j \in \delta t_i$, where $i, j = 1, ..., 23$. This counting forms *Si* sets, but since we select only the unique sets, we obtain *Sk*, where $k < i$. Counting is performed for each of four lunar models, and final result is a union. Therefore, if $Sk_{G11}$, $Sk_{issi2}$, $Sk_{iss3}$, and $Sk_{K14}$ are sets obtained for four lunar models, we proceed as $Sk = Sk_{G11} \cup Sk_{issi2} \cup Sk_{iss3} \cup Sk_{K14}$, for each *k*. In Figure 1, we show distance measurements based on differential travel times $t_{sp}$ in time (Figure 1a) and space (Figure 1b) domain for each nest and four lunar models. The formed sets are listed in Table 1. Since, we consider union of sets from different lunar models, nests such as A8 and A14 belong to a set S12 but also to sets S9, S10, S11, and S12 due to the estimates of one lunar model $t_{sp}^{issi3}$ shown in Figure 1a.

Finally, these sets signify that if we are able to infer $t_s$ and $t_p$ of the new lunar event with an accuracy within 5 s, we are not able to uniquely differentiate between nests that belong to one particular set Sk. It is important to emphasize that the number and content of sets strongly depend on the picking error threshold. If this threshold is smaller than 5 s, we shall obtain some sets that contain one nest, meaning that those nests can solely be constrained with travel time measurement. On the contrary, having larger threshold, we shall obtain sets that contain larger number of nets, making the classification task undoubtedly more difficult. In the literature, it has been shown that the error picks for travel time arrivals vary between different studies and for different events. The overview of the existing travel times error picks has been recapped in Figure 11 by Nunn et al. (2020a), and generally studies would usually favor less rigorous picking error threshold to consider broad variations in the arrival times of S waves, such as was studied by Hempel et al. (2012), where a 20 s error picking threshold was chosen. Nevertheless, we have reason to believe that 5 s error is a good approximation for the future lunar mission in situ instruments especially due to their broadband capability and better performances (Yamada et al., 2011, 2013). Therefore, it is still useful to approximate distance location of nests with this approach. To further tackle the nest classification into more details, we proceed to associate each event with a combination of lunar orbital measurements.

## 2.2. Feature Selection Based on the Lunar Orbital Information

It has been shown that the DMQ temporal patterns in time occurrences are related to different lunar cycles and that these patterns differ from nest to nest. These temporal patterns reflect physics governed mainly by the three-body problem of the Earth, Moon, and Sun. Therefore, it is our intention to exploit the physics behind these temporal patterns by designing features that directly describe the Moon's motion. Naturally, to do so we chose features that are independent from each other and have physical meaning. These are synodic, draconic, and anomalistic, and they all have similar periods but are marked by different motions, either as the motion between two Full Moons phases, or two nodes (descending or ascending as explained above), or two apsis, respectively. One can list all the events when Moon is in the Full Moon (New Moon) phase, passing through ascending (descending) node or perigee (apogee) by simply looking at the Moon's ephemeris (Meeus, 1991). To make sure that we take into account the temporal patterns, we design the main three features as a time difference between the time of the quake in the nest and the time of the Moon's Full Moon, ascending node and perigee, denoting it as $\Delta t_{FM}$, $\Delta t_{AN}$, $\Delta t_{PG}$, respectively. We can achieve the same effect by taking the other three time axis as referent one (New Moon, descending node, and apogee). The next feature is related to the Moon's position within its orbit as in Frohlich and Nakamura (2009). The angle between the direction of perigee and the current position of the body, as seen from the main focus of the ellipse, is called the true anomaly denoted further as $\gamma$. Furthermore, as one of the feature, we also use the interval time between two quakes in the nest, noted as $t_{i+1} - t_i$ in Weber et al. (2010). The last two features are related to the position of the Moon with respect to the Earth, and these are the distance, *d*, and the rate of the distance change, $\dot{d}$, as given in Bills et al. (2008).

The selected features all have different ranges and we refer to them as raw data. To train a model that is able to generalize well for a given problem sometimes, it is necessary to transform raw data to a form that is more suitable for training (Langer et al., 2019). By applying transformation on the raw data, we may obtain a mapping, which better reveals patterns in our data. Therefore, we chose to apply trigonometric transformation of the true anomaly angle $\gamma$ to properly address the jump discontinuities in the feature when angle goes from $2\pi$ to 0 due to its cyclic
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Location study of the DMQ nests from the perspective of $t_{sp} = t_s - t_p$ travel time measurements if we place station in Schrödinger basin and consider four different lunar models. (a) $t_{sp}$ travel time measurements for four lunar models from Garcia et al. (2011) (G11), Garcia et al. (2019) (issi2, issi3), and Khan et al. (2014) (K14) with nest labels alongside with $t_{sp}^{G11}$ travel time values plotted as $t_{sp}^{G11}$ versus. $t_{sp}^{issi1}$, $t_{sp}^{G11}$ versus. $t_{sp}^{issi2}$, $t_{sp}^{G11}$ versus. $t_{sp}^{K14}$. (b) Lunar map with the nests locations from Lognonné et al. (2003) where the color indicates the median $t_{sp}$ for the four lunar model.</figcaption></figure>

nature. This is addressed by transforming true anomaly angle $\gamma$ to pair of [cos $\gamma$, sin $\gamma$]. Finally, we can characterize each DMQ with eight features, which are $\Delta t_{AN}$, $\Delta t_{FM}$, $\Delta t_{PG}$, cos $\gamma$, sin $\gamma$, $t_{i+1} - t_i$, $d$, $\dot{d}$ summarized in Table 2. To proceed with the training of a machine learning model for the nest classification, we need to calculate these eight features for all quakes that belong to all nests that share similar travel time measurements shown in Table 1. An example of all eight features are shown in Figure 2 for all DMQs from the largest nest A1.
<table><caption>**Table 1**
List of Sets Sk Obtained by Classifying Nests' $t_{sp}$ Travel Time Measurements</caption><thead><tr><th>Set</th><th>Nests</th></tr></thead><tbody><tr><td>S1</td><td>A1, A11, A16, A97</td></tr><tr><td>S2</td><td>A1, A11, A97</td></tr><tr><td>S3</td><td>A6, A17</td></tr><tr><td>S4</td><td>A6, A17, A44</td></tr><tr><td>S5</td><td>A7, A8, A14, A18, A20, A25, A30, A41, A42, A50</td></tr><tr><td>S6</td><td>A7, A18, A20, A27, A30, A34, A41, A51</td></tr><tr><td>S7</td><td>A7, A18, A20, A27, A34, A51</td></tr><tr><td>S8</td><td>A7, A18, A27, A34, A40, A51</td></tr><tr><td>S9</td><td>A8, A14, A18, A20, A25, A30, A41, A42, A50</td></tr><tr><td>S10</td><td>A8, A14, A20, A25, A30, A33, A41, A42, A50</td></tr><tr><td>S11</td><td>A8, A14, A20, A25, A30, A41, A42, A50</td></tr><tr><td>S12</td><td>A8, A14, A33</td></tr><tr><td>S13</td><td>A9, A24</td></tr><tr><td>S14</td><td>A11, A16</td></tr><tr><td>S15</td><td>A17, A44</td></tr><tr><td>S16</td><td>A27, A34, A40, A51</td></tr></tbody></table>

Note. Each set consists of different combination of nests A1, and we are not able to differentiate between those nests that belong to a particular set.

<table><caption>**Table 2**
List of the Physics Inspired Features That Reflect the Motion of the Moon Around the Earth Calculated for the Reference Time, Which Is the Time of the Quake $t_{quake}$. Within the Nest, Used as an Input Data to Train a Machine Learning Model for the Nests Classification</caption><thead><tr><th>Feature</th><th>Description</th></tr></thead><tbody><tr><td>Δt<sub>AN</sub></td><td>Time difference with the ascending node<br>|t<sub>Ascending Node</sub> - t<sub>quake</sub>|</td></tr><tr><td>Δt<sub>FM</sub></td><td>Time difference with the Full Moon<br>|t<sub>Full Moon</sub> - t<sub>quake</sub>|</td></tr><tr><td>Δt<sub>PG</sub></td><td>Time difference with the perigee<br>|t<sub>Perigee</sub> - t<sub>quake</sub>|</td></tr><tr><td>cos γ</td><td>Cosine of the true anomaly angle γ</td></tr><tr><td>sin γ</td><td>Sine of the true anomaly angle γ</td></tr><tr><td>t<sub>i+1</sub> - t<sub>i</sub></td><td>Time difference between two successive quakes</td></tr><tr><td>d</td><td>Distance between the Moon and the Earth</td></tr><tr><td>dot</td><td>Rate of the distance change between the Moon and the Earth</td></tr></tbody></table>

Note. γ is the true anomaly angle defined as the angle between the direction of the perigee and the position of the Moon at time of $t_{quake}$.

### 3. Methodology

When new lunar data arrive, the first objective is to classify events in groups based on the waveform similarity measurements. The next objective is to infer their P and S travel times and thus form sets of nests, as in Section 2.1. The final step would be to associate these new events from a single station, with the existing Apollo nests if possible. This can be achieved using supervised classification. The model is trained in a predictive way by taking into account nest locations as labels and nest lunar orbital parameters as input data. Since, we want to predict a class (nest), but we do not have statistically large data set (as previously mentioned the largest nest A1 has 443 quakes), we choose to train a random forest (RF) classifier, since RF can perform well with any size of data sets and tend not to overfit (Breiman, 2001; Ho, 1995). The method also preserves non-linear relationship between features; thus, we are not excluding the possible non-linear link between features and nests.

Random forest (RF) is a machine learning technique that is based on decision trees (Breiman et al., 1984; Quinlan, 1986) and bootstrap aggregating (Breiman, 1996), where the main output is reached by majority votes among an ensemble of randomized decision trees. A main building unit, a decision tree, is a tree-like learning algorithm that consists of decision (internal) nodes, inequality branches, and leaf nodes that hold the final prediction of the individual trees shown in Figure 3. Thus, within each tree, the beginning is at root node that doesn't have incoming branches. Next in line are internal nodes where based on the available features and inequality operations, the decision whether the feature is smaller or larger than some threshold is made. These translate to leaf nodes, which represent all possible outcomes. Usually, during the training phase thresholds, order and number of inequality operations within internal nodes are learned. The hyperparameters that define an RF structure and need to be defined before a training process are: the number of decision trees, the maximum depth of trees, the measure that maximizes diversity between classes, the minimum samples in the internal node, the minimum number of samples in the leaf node for it to be considered, the maximum number of features when looking for the best split, the maximum number of leaf nodes, and the maximum samples to be drawn from the main training data set when training each decision tree. In our analysis, we define a base RF model as the one trained with 1,500 decision trees, we do not set the maximum depth, but we let that the nodes expand until all leaves include only 1 data point, the quality of the split is measured by Gini impurity (Breiman et al., 1984), the minimum samples in the internal node is set to 2, the minimum number of samples in leaf node is set to 1, the maximum number of features when looking for the best split is set to $\sqrt{8}$ due to eight features we use, we do not set the maximum number of leaf nodes, and the maximum samples to training each decision tree depend on the size of the input data. Please refer to Text S3 in Supporting Information S1 for more details on the testing of our base model.

RF also provides an assessment of the feature or input variable importance, which might give us an insight of how the model reached its prediction. To assess the feature importance, the RF removes one of the features while it keeps the rest constant, and it measures, among others, the accuracy decrease (Breiman, 2001). RF models are able to solve regression and classification problems as well as two- and multi-class problems. It has been shown that RF can perform with high accuracy with only a few parameters to tune (Probst et al., 2019).

In our case, during the training phase, the RF model has access to the extracted features of the individual quakes within the nest and their labels. In
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Time evolution of features during the Apollo mission (left column) and their histograms (right column) for A1 nest. First three features are time difference between quakes and (a) the instance when Moon was passing through ascending node noted as $\Delta t_{AN}$, (b) Full Moon phase noted as $\Delta t_{FM}$, (c) instance when Moon was in perigee noted as $\Delta t_{PG}$, next (d) cos $\gamma$, (e) sin $\gamma$, where $\gamma$ is the true anomaly angle, indicating the position of the Moon in the orbit, (f) time difference between two quakes noted as $t_{i+1} - t_i$, (g) distance between the Moon and the Earth at the quake date noted as *d*, and (h) the rate of distance change from between the Moon and the Earth noted as $\dot{d}$.</figcaption></figure>

the context of classification, we are solving a multi-class problem since some sets have more than two nests to predict (see Table 1). The training is performed on a subset of the data, while the model performance is evaluated on the test subset, which the model has never seen to avoid overfitting of the data. Evaluation is accomplished by comparing the model's predicted class (nest) with the ground truth one. The statistical performance of the model is presented with the elements of the confusion matrix, receiver operating characteristic (ROC) curve, and a set of measurements such as accuracy, precision, recall, and f1-score (see more details and statistical definitions in Text S3 of the Supporting Information S1). In our analysis, we use 75% of the original database for training and the remaining 25% for testing.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Schematic representation of a random forest algorithm. In this example, the model is trained with 3 features, it consists of 3 decision trees with a maximum tree depth equal to 3. The tree consists of decision nodes (circles), followed by inequality branches (dashed lines), and leaf nodes (rectangles). The prediction is taking place in each tree by yes or no questions, while the final prediction is made upon majority voting considering individual tree predictions.</figcaption></figure>

We expect that in the case of the ideal RF classifier the elements of the confusion matrix that signify belonging to the correct class are maximum, and those elements that signify belonging to an incorrect class are zero. Similarly, an ROC curve is passing through the left upper corner for an ideal classifier, and all other statistical measurements should be maximized. Accuracy is the performance matrix that denotes how well the model predicted the correct class with respect to the total number of classes. Precision indicates how well model retrieved the relevant instances over the retrieved instances. Recall shows what is the fraction of the relevant instances that were retrieved (Powers, 2020). F1-score is the harmonic mean of the precision and recall. Since model can have different values of precision and recall, high value of f1-score would indicate model that perform well in both metrics. Accuracy is usually used to report how the overall model performs, while precision, recall and f1-score are used to report how the individual classes score.

# 4. Random Forest Application for DMQs Nests Classification and Its Implication for DMQs Nests Localization

## 4.1. Training and Testing on Two Largest Nests

We first test the hypothesis whether it is possible to differentiate two DMQ nests using the lunar orbital parameters (features). For this, we select the two largest nests, A1 and A8, with a total size of 768 events and ratio A1:A8 = 57%: 43% (see feature distributions in Figure S6 of the Supporting Information S1).

The first training and testing of our base RF model as defined in previous Section 3 is performed using the raw features as defined in Section 2.2 and using the normalized features as an input data. In our case, normalization is performed by scaling all features within the range of 0 and 1. The base model trained with the normalized data performed better, as shown in Figure 4 (compared with the performance of the RF model trained with the raw data
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Performance of the base RF model trained to classify between nests A1 and A8 using the raw normalized feature data: (a) A1–A1 and A8–A8 elements of the confusion matrix that signifies how much A1 and A8 quakes are correctly predicted as A1 and A8 nests, respectfully (these positive classifications are indicated with plus marker), A1–A8 and A8–A1 elements of the confusion matrix that shows how much A1 quakes are predicted as A8 nest, or A8 quakes are predicted as A1 nest, respectfully (these negative classifications are indicated with minus marker), precision, recall, f1-score per nests, and accuracy of the model (b) ROC curve, (c) feature importance attributes used by the model to make a prediction.</figcaption></figure>

shown in Figure S7 of the Supporting Information S1). Therefore, for the rest of our study, we use the normalized features as an input data. The base model trained with the normalized input data performed with an accuracy of 89%, as shown in Figure 4a. We can notice that model is predicting the A1 nest better than A8 (see Figure 4a). This is indicated by having A1–A1 element of the confusion matrix that signifies how much A1 quakes are correctly predicted as A1 nest higher than A8–A8 element of the confusion matrix (how much A8 quakes are correctly predicted as A8 nest). This suggests that A1–A8 element, quakes A1 predicted as A8 nest, is smaller than A8–A1 element, where A8 quakes are predicted as A1 nest (for more details see Text S3 in Supporting Information S1). Precision, recall, and f1-score for the A1 nest are 88%, 94%, and 91%, respectively, and for A8 are 90%, 80%, and 85%. Hence, the A1 nest is better predicted by having a higher f1-score, yet with the A8 nest prediction we can be more confident due to a higher precision value. The ROC curve is above the random classifier curve, meaning that the base model is not randomly classifying A1 and A8 nests (see Figure 4b). Out of eight features, the first five most important are cos(γ), Δ*t<sub>PG</sub>*d, Δ*t<sub>AN</sub>*d, and *t*<sub>*i*+1</sub> − *t*<sub>*i*</sub> (see Figure 4c). We notice that cos(γ) is the feature with the most important contribution to the model learning. This might be because A1 and A8 have reversed distributions of this feature, thus the RF model might associates all quakes closer to one side of this distribution with A1 class while the other side with A8 class (see Figure S6d in Supporting Information S1).

We proceed into testing learning robustness of our base model in a series of experiments. These tests involve (a) changing the randomness of the bootstrapping initialization of the samples that are used when building decision trees, the randomness of the feature sampling when considering for the best split at each internal node as well as the randomness of the training and test data set split (see Figure S8 in Supporting Information S1); (b) changing
the optimal number of decision trees (see Figure S9 in Supporting Information S1); (c) equalizing the size of nests by randomly downsampling the larger A1 nest to be the same size as A8 nest (see Figure S10 in Supporting Information S1); (d) reducing the number of input feature data to the five most important ones of the base RF model, which are cos(γ), Δ*t<sub>PG</sub>*d, Δ*t<sub>AN</sub>*t<sub>i+1</sub> - t<sub>i</sub> (see Figure S11 in Supporting Information S1). In all tests, except those in experiment (a), we keep the random state fixed. Finally, the results do not vary between different tests, indicating that in all above configurations models are able to learn how to classify two nests with the similar performances. This indicates that the model is stable and further implies that the base model generalizes well and does not overfit the results. Furthermore, we notice that in experiment (c) the prediction of A8 nest improved compared to the prediction obtained with the base RF model, indicating that having balanced classes, while training a ML model, is important (see Figure S10 in Supporting Information S1). Moreover, we observe that the model trained with fewer features performs worse than the base RF model (see Figure S11 in Supporting Information S1). This might be because all eight features are uncorrelated, implying that there is no redundant information within the input data. We also use a cross-validation technique on our base model. With this technique, different portions of the training data are used for training and validation during several iterations. Usually, the training data set is divided into *k* equally sized folds, and then *k* iterations are performed. In each iteration, (*k* − 1) folds are used for training, and one fold is used for validation. During the cross-validation, the test set is kept aside. Eventually, the full data set is divided into three sets: training (55%), validation (20%), and test (25%), where training and validation set change in each iteration. We calculate cross-validation with *k* = 5 while keeping the same hyperparameters as the base RF model. The choice of *k* = 5 has been proven to be a good practice (Witten & James, 2013). This test produces 5 models with equal performance, as indicated with ROC curves (see Figure S12 in Supporting Information S1). Additionally, we also perform a grid search analysis, where we search for the best model by sampling a small space of hyperparameters (for the list of parameters see details in Text S3 of the Supporting Information S1). Yet, the base model still performed better than the grid search favorite (see Figure S13 in Supporting Information S1). Additionally, it is worth emphasizing that with each class prediction the base RF model also associates the class probability. This value ranges between 0 and 1, where 1 indicates that a model is absolutely certain that a given quake belongs to a predicted class (nest). If we consider the distribution of these probability values, we notice that for our base RF model more than 81% are higher than 0.80 for the test data set (see Figure S14a in Supporting Information S1). The mislabeled prediction values show a uniform distribution between 0.5 and 1 (see Figure S14b in Supporting Information S1). Furthermore, we might try to unravel why the base model sometimes misclassifies nests by predicting A1 instead of A8 and vice versa. To do so, it seems reasonable to look at the features with the highest importance attributes (cos(γ), Δ*t<sub>PG</sub>*d shown in Figures S14c–S14e of the Supporting Information S1, respectively). Yet, we could argue that the answer is not straightforward from the perspective of single feature, since these misclassified events are the consequence of non-linear patterns presented probably in all three features.

## 4.2. Training and Testing on Three and More Nests

In this section, we study how the performance of our base RF model from Section 4.1 changes by adding more nests. This is motivated by the fact that the classification problem we are solving is quite unbalanced, since nests' have rather different sizes (see Figure S1 in Supporting Information S1). Therefore, it should be insightful how the performance of the classifying models change by adding more nests into original configuration from Section 4.1, where the sizes of added nests are gradually decreasing. We carry out three tests for the next nests combinations and their ratios: A1–A8–A24 (44%–33%–23%), A1–A8–A24–A18 (36%–27%–19%–18%), and A1–A8–A24–A18–A6 (31%–24%–17%–15%–13%), where the three added nests are three largest nests after A1 and A8 (see Figure S1 in Supporting Information S1).

Analysis shows that by adding more nests, the performance of our base model deteriorates since the accuracy drops from 88% to 77%, as shown in Figure 5a. By adding a 3rd nest, we notice that A1 recall improves and A8 recall deteriorates slightly (see Figure 5a), and 50% of A24 events are classified either as A1 or A8 with the recall of 53% (see Figure 5b). The precision of A24 is the highest. Yet, the overall f1-score is the lowest for A24, and according to this measure nest, A1 has the most balanced precision and recall values (see Figure 5a). Five features, Δ*AN*, Δ*PG*, cos(γ), tᵢ₊₁ - tᵢ, and *d*, have the same feature importance (see Figure 5c). Thus, it seems that by adding a third nest, the importance of all features becomes more equalized. This might reflect the complexity of classifying three instead of two classes, since now the RF model cannot rely on limited features to make a decision, since three source mechanisms we are classifying might be responsive to the same orbital parameters.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Performance of the base RF model trained to classify between 4 groups of nests: A1–A8 (2 nests), A1–A8-A24 (3 nests), A1–A8–A24–A18 (4 nests), and A1–A8–A24–A18–A6 (5 nests): (a) accuracy of each model and precision, recall, f1-score of individual classes (nests), (b) elements of the confusion matrix (with minus marker for positively classified quakes and plus marker for negatively classified quakes), and (c) feature importance attributes used by the models to make a prediction.</figcaption></figure>

By adding a 4th nest, the recall of A1 nest improves, recall of A8 nest deteriorates even more than before recall of A24 nest deteriorates, and the new added nest A18 has a recall of 52% (see Figure 5a). The precision values of three nests do not exceed 70%, and f1-score of A1 nest is still dominating over other nests. Most of the events of
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Accuracy performance of the base RF model trained to classify between nests belonging to different sets from Table 1. Individual nests can belong to multiple sets as presented in Table 1. Sets with similar number of nests are represented with the same color, with the accuracy value above the bars.</figcaption></figure>

A18 nest are misclassified as events from A1 nest (see Figure 5b), with less than 1% of events being classified as events from A8 to A24 nests. This might imply that A18 nests have a source mechanism that is the most similar to the source mechanism of A1 nests. Less notably than before, the importance of all features is still equalized (see Figure 5c).

Lastly, by adding a 5th nest, the recall of A1 nest is the highest, while A8, A24, and A18 nests perform with recall less than 51% and the newly added A6 nest has recall of 63% (see Figure 5a). The precision of the newly added nest is the highest, while precision of A1 nest still decreases (see Figure 5a). We notice that A1 and A6 nests have the highest f1-score (around 70%), while other nests are performing worse by more than 10%. This might imply that A6 nest has a very unique pattern that helps this nest to compete with the most dominant nest in this configuration in terms of the number of quakes, which is A1 nest. The most quakes of A6 nests are mislabeled as quakes of A1 nest, while not a single event is mislabeled as A24 or A18 nests (see Figure 5b). We also notice that the importance of feature $t_{i+1} - t_i$ significantly changes (see Figure 5c).

If we remove A1 nest and train a model to classify between A8–A24–A18–A6, we notice that overall recall of all nests, besides A24 nest, improves (see Figure S15 in Supporting Information S1). Also, A16 still performs the foremost based on f1-score, and the feature importance is dominated by $\Delta t_{AN}$ (see Figure S15 in Supporting Information S1). We notice that interplay between features changes by removing A1 nest and by allowing A6 nest to dominate the prediction, which is shown by A6 nest having the highest f1-score.

These results might imply that by adding more nests, we add more complexity into the problem we are solving, since we might be adding nests that have similar source mechanisms. Also, having similar source mechanisms means that sources are triggered by tides in the same way, so their lunar orbital features have similar characteristics, and we cannot differentiate between nests without having more data. This directly affects the statistic measures of individual nests, and it is not obvious which nest should be having the foremost performance before training the classification model. Generally, having high values of precision over recall for most nests in all four groups shows that we are able to classify instances of particular nest correctly, while this fraction is usually 50%. In case where these results might be used to improve the lunar models, precision over recall is more important, since we want to be sure that we are predicting a correct location of the new quakes.

### 4.3. Training and Testing on Nest Sets

We have shown in Section 2.1 that we cannot differentiate between nests that belong to the same sets listed in Table 1 using only S-P travel times. Therefore, here, we use the same base RF model as in Section 4.1 to train and test how well we can differentiate nests that belong to the same sets using lunar orbital parameters. For ease in visualization, we present results in three groups, where we group sets by the number of nests they contain but also by the frequency of the repeating nest within the set. Therefore, we form group A with S1, S2, S3, S4, S12, S13, S14, S15, and S16 set that contain 4 or less nests, group B with S5, S6, S7, and S8, and group C with S9, S10, and S11 that share many common nests.
In Figure 6, we observe that all models reach accuracy values above 50%. Those sets from group A that contain 4 or less nests (blue bars in Figure 6) achieve accuracy values higher than 70%. If we look how the individual nests within different sets perform, we notice that most nests have high precision values as shown in Figure 7 as well as high recall shown in Figure 8 and f1-score (see Figure S16 in Supporting Information S1) values. Mean values of precision, recall, and f1-score for all nests are 82%, 59%, and 63%, respectively. In general, sets that have ≤4 nests perform better than those with more nests, as shown in Figure 7a for group A. Almost all nests from group A score with precision values higher than the average but A51, A34, and A44. Also, from the same group only A97 and A11 nests do not reach recall values higher than the average (see Figure 8a). Therefore, based on this result, we can determine that RF models are well trained to classify nests within the sets from group A. This further implies that locations of the nests from these sets are well constrained in a probabilistic sense. Since in a group A there are 9 sets, it means that for more than 10 individual nests we can predict a preliminary probabilistic location for the future lunar missions. As seen in Section 4.2 having more than 4 nests, we add more complexity into our classification problem, and this is reflected on the performance for group B (see Figures 7b and 8b) and group C (see Figures 7c and 8c). We still notice that many nests score with precision values higher than the average (condition is not met for A8, A7, and A18 nests shown in Figure 7b and A8, A18 nests shown in Figure 7c). Most nests from group B also have recall values higher than the average recall value (condition not met for A14, A51, A41, and A30 nests in Figure 8b), while for group C it is the opposite, most nests do not meet this condition (as A18, A20, A41, A30, and A42 nests shown in Figure 8c). When the recall values are very low or zero (as for A11 nest in Figure 8a, A51, and A30 nests in Figure 8b, A41, A30, A42, and A50 nests in Figure 8c), it signifies a nest with very few events (see ratio of nests in all sets in Figure S17 in Supporting Information S1). Nevertheless, if we take an example of A20 nest, we notice that it has constant recall in many sets (see Figures 8b and 8c) even though it is not the biggest nest in the set (see Figure S17 in Supporting Information S1). Thus, not only the size but also the uniqueness of the features determines the success of identifying the nest. As mentioned earlier, to be sure that we are correctly locating different nests within sets, we prefer that the model performance scores better in precision over recall metric.

Feature importance attributes for each model trained to classify nests within sets from Table 1 are shown in Figure 9 again in three groups. We notice that, on one hand, removing just one nest could change the feature importance, as in the case of S2 (where we remove A16) versus S1. On the other hand, we notice that the feature importance does not drastically change when comparing results for sets S3 and S4, where we add nest A44, even though the nest itself is large in size (see Figure S17 in Supporting Information S1). For the sets in the group B (S5, S6, S7, S8), the feature importance is stable with respect to adding or removing nests. This is very true for the group C (S9, S10, S11). We notice that sets which contain ≤4 nests (as in group A), there is usually one or two important features, while for sets with more than 4 nests there is homogenization of the feature importance (as in groups B and C). Having fewer nests might imply that a single lunar orbital parameter is enough to explain the occurrence of the nests, which are unique in nature. Moreover, the feature importance changes for different sets that have unique combinations of nests, which may hint that these nests have different source mechanisms. Mixing more nests suggests that we might be mixing nests with similar temporal patterns, thus learning how to classify them is more challenging.

The feature importance attributes for each model trained to classify nests within sets from Table 1 is shown in Figure 9a again in three groups. To better understand how different features are from one set to another, we calculate Pearson correlation coefficient for each pair of set following this expression

$$C_{ij} = \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\left(\sum_{i=1}^{n} (x_i - \bar{x})^2\right) \left(\sum_{i=1}^{n} (y_i - \bar{y})^2\right)}} \quad (1)$$

where n = 8 stands for eight features shown in Table 2, xᵢ stands for features associated with model trained to classify nests in set Sᵢ, and yᵢ for features associated with model trained to classify nests in set Sⱼ with i ≠ j, x, y stand for the mean values. The correlation coefficients varies between -1 and 1, where it represents how similar are two sequences with a value of 1 signifying the perfect match between the two, -1 representing perfect anti-correlation, and 0 no correlation at all. These coefficients calculated between all pair of sets are presented in Figure 9b, where the diagonal is equal to 1 since it is calculated for the pair of the two same nests. It is interesting to notice that removing just one nest could change the feature importance as in the case of S2 set, where we
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Performance of RF models designed to classify nests within different sets in terms of precision metric that each nest scored within different sets from Table 1. Individual nests can belong to multiple sets as presented in Table 1. Sets are indicated with the same color, nests are indicated alongside with the travel time distance based on Garcia et al. (2011) (G11) lunar model ($t_{sp}^{G11}$), the precision value is indicated for each nest, and those nests that have precision value equal to 0 are indicated with cross marks. Results are presented in three groups: (a) S1, S2, S3, S4, S12, S13, S14, S15, S16; (b) S5, S6, S7, S8; (c) S9, S10, S11.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Performance of RF models designed to classify nests within different sets in terms of recall metric that each nest scored within different sets from Table 1. Individual nests can belong to multiple sets as presented in Table 1. Sets are indicated with the same color, nests are indicated alongside with the travel time distance based on Garcia et al. (2011) (G11) lunar model ($t_{sp}^{G11}$), the recall value is indicated for each nest, and those nests that have recall value equal to 0 are indicated with cross marks. Results are presented in three groups: (a) S1, S2, S3, S4, S12, S13, S14, S15, S16; (b) S5, S6, S7, S8; (c) S9, S10, S11.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 9.** (a) Feature importance for RF models trained to classify nests within sets from Table 1. (b) Correlation coefficients calculated between feature vectors for all pair of sets from (a). Individual nests belong to multiple sets as presented in Table 1.</figcaption></figure>

removed nest A16 compared to S1 set. Thus, the correlation coefficient between features of sets S1 and S2 is $C_{1,2} = 0.87$. Next, we notice that the feature importance does not change drastically when comparing results for sets S3 and S4, where we add nest A44, even though the nest itself is large in size (see Figure S17 in Supporting Information S1). Thus, their correlation coefficient is $C_{3,4} = 0.95$. For the sets in group B (S5, S6, S7, S8), the feature importance is stable with respect to adding or removing nests, except for S8; therefore, we have the
correlation coefficients for C<sub>5,6</sub>, C<sub>5,7</sub>, C<sub>5,8</sub>, C<sub>6,7</sub>, C<sub>6,8</sub>, C<sub>7,8</sub> as 0.94, 0.96, 0.75, 0.97, 0.71, and 0.83, respectively. It seems that nest A40 makes a difference for set S8 when compared to sets S5, S6, and S7. For group C (S9, S10, S11), the feature importance is quite stable, having the correlation coefficients C<sub>9,10</sub>, C<sub>9,11</sub>, C<sub>10,11</sub> as 0.98, 0.98, and 0.99. We notice that for sets which contain ≤4 nests (such as group A), there are usually one or two important features, while for sets with more than 4 nests there is homogenization of the feature importance (as in groups B and C). Having fewer nests might imply that a single lunar orbital parameter is enough to explain the occurrence of the nests, which are unique in nature. Moreover, the feature importance changes for different sets that have unique combinations of nests, which may hint that these nests have different source mechanisms. Mixing more nests suggests that we might be mixing nests with similar temporal patterns, thus learning how to classify them is more challenging. It is worth noting that reported feature importance values are only indicative of the importance of features for the particular RF classifier but not suggesting that they are important in more general sense outside of the machine learning context. Furthermore, since modeling always involves two or more nests, it is ambiguous to give interpretation for a particular nest. Taking this into account, we further investigate the correlation between feature importance and location for the sets in group A. From this group, the nests that belong to sets S13 and S3 are the most separated in terms of travel time calculation, where S13 is closer to the landing location. The RF classifier associated with these nests are feature-wise the most sensitive to different values, δt<sub>AN</sub> versus t<sub>i+1</sub> - t<sub>i</sub>, respectively. The correlation coefficient, between features associated with RF classifiers of S13 and S3 nests, C<sub>13,3</sub> = 0.22 is very low. Therefore, in the context of the feature we selected, the nests we are discussing, and RF classifier, we can hypothesize that the southern (where S13 sets is) and the northern (where S3 nests is) lunar hemisphere respond differently to the tidal forcing. We can perform a similar analysis between S13 and S12 sets, which are closest in travel time distance. We notice that their correlation coefficient is C<sub>13,12</sub> = 0.48, slightly larger than previous one, and the most important feature for the RF classifier of S12 nest is sin(γ). For further comparison, the correlation coefficient of the feature vectors related to S12 and S3 sets is C<sub>12,3</sub> = 0.32, which is consistent with the previous observation. Next, we notice that the largest correlation coefficient related to the feature vector of S13 set is with S2 set C<sub>13,2</sub> = 0.60, while travel time of S12 sets is half-way between S13 and S3. The correlation coefficient between features related to S2 and S3 sets is C<sub>2,3</sub> = 0.67. If we assume that half-way between S13 and S3 nests, there is a region which reflects characteristic of both southern and northern lunar hemisphere, then having correlation coefficients C<sub>13,2</sub> = 0.60 and C<sub>2,3</sub> = 0.67 support this statement. The physical interpretation of nests has been tackled in paper by Weber et al. (2009), where they showed that nests are likely to activate when the stress vector is aligned with the fault orientation. Therefore, even at similar distance, nests may act differently. For the nests A9, A24, A6, A17, A8, A14, A33 associated with the sets S13, S3, S12, we can notice the relative importance of stress rates, assessed through two failure criterion equations, differs between them. According to Table 3 in the paper Weber et al. (2009), the nests A6 and A33 are better fit by Equation 1, while the nests A9, A17, and A8 are better fit by Equation 11, where Equation 1 and 11 can be found in Weber et al. (2009). Therefore, considering just the set S3 with A6 and A17, we notice that these two nests respond differently to stress rate even though they share similar travel times. This may further support our results that is why RF models are able to distinguish between different patterns associated with different nests within sets.

We observe that different ranges of distances are well constrained with this analysis. In Figure 1, we notice that nests considered in this analysis, on average, have travel time distance measurements ranging from 160 to 240 s. Sets from group A (S1, S2, S3, S4, S12, S13, S14, S15, S16) which have good performances also span the same time distance range. While this means that we may not be able to differentiate all nests, we may have a reasonable set of nests to work with to improve our existing lunar models.

4.4. Discussion on the Future Mission Application

This study shows us that with partial training data sets of Apollo DMQ data, we have enough information to
achieve agreeable prediction models of different nests, allowing for preliminary quake source locations. For a
future deployment such as FSS, the base RF model would be trained with the full existing DMQs data set,
meaning that we take into account the whole activity history of each particular nest instead of just 75% of the
original database as we did in this analysis to avoid the overfitting of the models. Additionally, in practice once
when the new data start to arrive the quake waveforms will be also important. As discussed, in this study, we do
not directly consider waveform information as a part of the prediction model, since the recording station would be
placed on the new location on the farside of the Moon with an instrument with a different frequency response,
where presumably the waveforms coming from the known nests wouldn't match the waveforms from the Apollo
locations. Nevertheless, the waveform information is considered here indirectly by acknowledging the existing
quake labels from the study by Nakamura (2003). Therefore, once the new data arrive the cluster analysis using
the events' waveforms will be performed. This would be the first step of the proposed detection pipeline. The
cluster analysis together with the travel time measurement tsp would inform us of the existing sets of nests. If it
turns out that some sets have more nests than shown in Table 1, we would still have a general idea which nests are
more likely to be active due to the existing statistics from the Apollo missions. Meaning, the quakes coming from
the A1 nest should still be statistically the most detected quakes (see nests' size in Figure S1 of the Supporting
Information S1). After cluster analysis and tsp measurements, next we would apply the RF prediction model.
Therefore, one needs to keep in mind that the trained RF models are essentially providing us with the probability
values of the nests location in the context of accepting the waveforms similarity hypothesis tsp values and sets they
have been previously assigned to. As previously described, the probability ranges between 0 and 1. If some
quakes are assigned a probability value that is smaller than 0.5 threshold, we can speculate that this quake does not
belong to either class that we are trying to predict. To assess how our models are sensitive to quakes that are
outside of the trained data set for a particular RF model, we perform a test where we apply an RF set models to all
quakes that do not belong to that particular set and we look what is the maximum probability value. This test might
be biased, since we do not consider tsp value differences; therefore, we might be mixing nests with similar fault
mechanism. Yet, the results shown in Figure S18 of the Supporting Information S1 indicate that for all models,
even with +σ, we have a statistical window where we can speculate that quakes might not belong to the nests they
have been predicted to. Furthermore, we have reason to believe that every time a new quake from the known nests
is detected, our RF predictions should improve with the updated data set. The proposed models are light in
training, and with the new events, we shall be able to retrain them and follow their statistical performance.

There is a possibility to extend the presented approach by taking into account the amplitude of quakes as an additional feature and not only the occurrence time. By this, we would consider the relative size of the quakes within the nest but also the energy released by the individual nests. However, this requires a bit of precaution, since it would require a scaling law between the amplitude from the Apollo missions and the one from the FSS, if we wish to directly apply the RF model trained on the Apollo data on the new data. Also, we might be able to reformulate the approach in the future into more casual problem, where we would consider which quake preceded or succeeded which one. Again, there is precaution of applying this on Apollo data due to the instrument sensitivities and possibility that the lunar seismicity was not well captured. Both extensions might be better suited once the new catalog would be derived from the future missions. In the context of the machine learning technique perspective, it might be interesting to tackle further the multi-class classification by training a single classifier solving one-to-one or one-to-all decomposition. The one-to-one decomposition splits a multiclass problem into one binary classification problem per each pair of classes, while the one-to-all decomposition splits multi-class problem into one binary problem per class (Del Moral et al., 2022).

**5. Conclusion**

In this paper, we propose how to tackle DMQ nest classification during future lunar missions that will likely host
only one station on the farside of the Moon. We propose constraining their location using differential time travel
measurement t<sub>sp</sub> and parameters related to the temporal patterns of the DMQ occurrence. First, in our analysis, we
assume that we cannot differentiate between nests whose differences in travel time are less than 5 s. Thus, we form
set of nests that have similar travel times. Second, for each quake within the nests, we calculate features that are
used to build a random forest model. This model is trained to classify between distinct nests. The features used for
training are built by associating each event in all nests with the time difference between events' origin time and the
time of lunar ascending node, Full Moon phase, perigee, then position of the Moon in its orbit expressed by true
anomaly angle, distance of the Moon from the Earth, rate change of this distance, and the time between two
successive quakes. We show that by training random forest models to differentiate between distinct nests within
sets, we can obtain models with high accuracy (more than half of the models score above 70% accuracy). Yet, the
performances of these models depend on the number of nests within the set. More nests implies that the problem is
more difficult to solve, probably because (a) nests might have similar source mechanisms, (b) the number of
events within nests is unbalanced, and (c) we don't have enough data. Since RF models also arrange features by
their importance to make a final classification decision, we observe that the importance of the features change
with different sets. This complements findings of previous papers, since it signifies that nests do correspond to
different lunar events, which eventually might be connected to the distribution of tidal stresses during these
events. Finally, our model provides a good first approximation of the nest localization. As the catalog of new events shall grow, it should be straightforward to retrain the RF model with the new enlarged data set and consequently improve the nest localization and classification.

## Data Availability Statement

We provide the DMQ RF Classifier script used to train RF models at Majstorović et al. (2024). The deep moonquake catalog used in this study is published in Nakamura et al. (1981) and revisited in Nunn et al. (2020a, 2020b). The models of lunar interior are published in Garcia et al. (2019). Version 2.6 of The TauP Toolkit, used for calculating travel times of P and S waves, is available via LGPL-3.0 license and developed openly at Crotwell et al. (1999). Version of The TauP Toolkit implemented as a part of Python library Obspy has also been used to calculate travel times of P and S waves can be found at Beyreuther et al. (2010). Version 1.17 of the Skyfield, used for calculating Moon's orbital parameters based on JPL ephemeris, is available via MIT license and developed openly at Rhodes (2019). For our implementation of the random forest algorithm, we use Scikit-learn machine learning Python library (Pedregosa et al., 2011). Figures were made with Matplotlib version 3.7.3 available under the Matplotlib license (Hunter, 2007). Version 0.22.0 of the Python package Cartopy, used for plotting map of Moon, is available via LGPL-3.0 and GPL-3.0 licenses and developed openly at Elson et al. (2023).

## Acknowledgments

J. Majstorović would like to thank Léonard Seydoux for his comments, suggestions, and discussions throughout this work and on an early draft of this manuscript. French co-authors thanks the French Space Agency, CNES, for supporting this research in the frame of the French contribution to FSS as well as PSUPC (Pole Spatial de l'Université Paris Cité), LABEX UnivEarthS (ANR-10-LABX-0023), and IDEX Paris Cité (ANR-18- IDEX-0001) for additional supports. MPP was supported by the NASA Farside Seismic Suite project and funds from the Jet Propulsion Laboratory, California Institute of Technology under a contract with the National Aeronautics and Space Administration (80NM0018D0004). We are grateful to three anonymous reviewers for their careful remarks that improved the manuscript.

## References

* Banerdt, W. B., Smrekar, S. E., Banfield, D., Giardini, D., Golombek, M., Johnson, C. L., et al. (2020). Initial results from the insight mission on Mars. *Nature Geoscience*, 13(3), 183-189. https://doi.org/10.1038/s41561-020-0544-y

* Benna, M., Sarantos, M., Schmerr, N., Malespin, C., & Bailey, S. (2020). The lunar environment monitoring station (LEMS). In *Lunar surface science workshop* (Vol. 22s41, p. 5022).

* Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., & Wassermann, J. (2010). Ob: A python toolbox for seismology. *Seismological Research Letters*, 81(3), 530-533. https://doi.org/10.1785/gssrl.81.3.530

* Bills, B. G., Bulow, R., & Johnson, C. L. (2008). Influence of Earth-Moon orbit geometry on deep moonquake occurrence times. *Lunar and Planetary Science*, XXXIX.

* Breiman, L. (1996). Bagging predictors. *Machine Learning*, 24(2), 123-140. https://doi.org/10.1007/bf00058655

* Breiman, L. (2001). Random forests. *Machine Learning*, 45(1), 5-32. https://doi.org/10.1023/a:1010933404324

* Breiman, L., Friedman, J., Stone, C. J., & Olshen, R. A. (1984). *Classification and regression trees*. CRC Press.

* Bulow, R., Johnson, C., Bills, B., & Shearer, P. (2007). Temporal and spatial properties of some deep moonquake clusters. *Journal of Geophysical Research*, 112(E9), E09003. https://doi.org/10.1029/2006je002847

* Bulow, R., Johnson, C., & Shearer, P. (2005). New events discovered in the apollo lunar seismic data. *Journal of Geophysical Research*, 110(E10), E10003. https://doi.org/10.1029/2005je002414

* Cameron, W. S., & Gilheany, J. J. (1967). Operation Moon blink and report of observations of lunar transient phenomena. *Icarus*, 7(1-3), 29-41. https://doi.org/10.1016/0019-1035(67)90044-9

* Cheng, C. H., & Toksöz, M. N. (1978). Tidal stresses in the Moon. *Journal of Geophysical Research*, 83(B2), 845-853. https://doi.org/10.1029/JB083iB02p00845

* Crotwell, H. P., Owens, T. J., & Ritsema, J. (1999). The taup toolkit: Flexible seismic travel-time and ray-path utilities. *Seismological Research Letters*, 70(2), 154-160. https://doi.org/10.1785/gssrl.70.2.154

* Cutler, J. W., Nguyen, T. A., Kano, T., Kumar, Y. R. L., Panning, M., April, S., & Haque, S. (2023). Overview of the avionics design for the farside seismic suite. In *AIAA SciTech 2023 forum*. https://doi.org/10.2514/6.2023-1880

* Dainty, A., Goins, N., & Toksoz, M. (1975). Natural lunar seismic events and the structure of the Moon. In *Lunar science conference*, 6th (Vol. 3, pp. 2887-2897). Pergamon Press, Inc.

* Del Moral, P., Nowaczyk, S., & Pashami, S. (2022). Why is multiclass classification hard? *IEEE Access*, 10, 80448-80462. https://doi.org/10.1109/access.2022.3192514

* Elson, P., Sales de Andrade, E., Lucas, G., May, R., Hattersley, R., Cambell, E., et al. (2023). SciTools/cartopy: v0.22.0 [Software]. Zenodo. https://doi.org/10.5281/zenodo.8216315

* Ewing, M., Latham, G., Press, F., Sutton, G., Dorman, J., Nakamura, Y., et al. (1971). Seismology of the Moon and implications on internal structure, origin and evolution. *Highlights of Astronomy*, 2, 155-172. https://doi.org/10.1007/978-94-010-3102-8_7

* Frohlich, C., & Nakamura, Y. (2009). The physical mechanisms of deep moonquakes and intermediate-depth earthquakes: How similar and how different? *Physics of the Earth and Planetary Interiors*, 173(3), 365-374. https://doi.org/10.1016/j.pepi.2009.02.004

* Gagnepain-Beyneix, J., Lognonné, P., Chenet, H., Lombardi, D., & Spohn, T. (2006). A seismic model of the lunar mantle and constraints on temperature and mineralogy. *Physics of the Earth and Planetary Interiors*, 159(3-4), 140-166. https://doi.org/10.1016/j.pepi.2006.05.009

* Garcia, R. F., Gagnepain-Beyneix, J., Chevrot, S., & Lognonné, P. (2011). Very preliminary reference Moon model. *Physics of the Earth and Planetary Interiors*, 188(1-2), 96-113. https://doi.org/10.1016/j.pepi.2011.06.015

* Garcia, R. F., Khan, A., Drilleau, M., Margerin, L., Kawamura, T., Sun, D., et al. (2019). Lunar seismology: An update on interior structure models. *Space Science Reviews*, 215(8), 1-47. https://doi.org/10.1007/s11214-019-0613-y

* Hempel, S., Knapmeyer, M., Jonkers, A., & Oberst, J. (2012). Uncertainty of Apollo deep moonquake locations and implications for future network designs. *Icarus*, 220(2), 971-980. https://doi.org/10.1016/j.icarus.2012.06.026

* Ho, T. K. (1995). Random decision forests. In *Proceedings of 3rd international conference on document analysis and recognition* (Vol. 1, pp. 278-282).

* Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. *Computing in Science & Engineering*, 9(3), 90-95. https://doi.org/10.1109/MCSE.2007.55
Kawamura, T., Grott, M., Garcia, R., Wieczorek, M., de Raucourt, S., Lognonné, P., et al. (2022). An autonomous lunar geophysical experiment package (ALGEP) for future space missions: In response to call for white papers for the voyage 2050 long-term plan in the esa science program. *Experimental Astronomy*, 54(2-3), 617-640. https://doi.org/10.1007/s10686-022-09857-6

Kawamura, T., Lognonné, P., Nishikawa, Y., & Tanaka, S. (2017). Evaluation of deep moonquake source parameters: Implication for fault characteristics and thermal state. *Journal of Geophysical Research: Planets*, 122(7), 1487-1504. https://doi.org/10.1002/2016je005147

Khan, A., Connolly, J. A. D., Pommier, A., & Noir, J. (2014). Geophysical evidence for melt in the deep lunar interior and implications for lunar evolution. *Journal of Geophysical Research: Planets*, 119(10), 2197-2221. https://doi.org/10.1002/2014JE004661

Khan, A., & Mosegaard, K. (2002). An inquiry into the lunar interior: A nonlinear inversion of the Apollo lunar seismic data. *Journal of Geophysical Research*, 107(E6), 3-1-3-23. https://doi.org/10.1029/2001je001658

Khan, A., Mosegaard, K., & Rasmussen, K. L. (2000). A new seismic velocity model for the Moon from a Monte Carlo inversion of the Apollo lunar seismic data. *Geophysical Research Letters*, 27(11), 1591-1594. https://doi.org/10.1029/1999gl008452

Lammlein, D. R. (1977). Lunar seismicity and tectonics. *Physics of the Earth and Planetary Interiors*, 14(3), 224-273. https://doi.org/10.1016/0031-9201(77)90175-3

Lammlein, D. R., Latham, G. V., Dorman, J., Nakamura, Y., & Ewing, M. (1974). Lunar seismicity, structure, and tectonics. *Reviews of Geophysics*, 12(1), 1-21. https://doi.org/10.1002/9781118728118.ch1

Langer, H., Falsaperla, S., & Hammer, C. (2019). *Advantages and pitfalls of pattern recognition: Selected cases in geophysics*. Elsevier.

Lognonné, P., Banerdt, W. B., Giardini, D., Pike, W. T., Christensen, U., Laudet, P., et al. (2019). SEIS: Insight's seismic experiment for internal structure of Mars. *Space Science Reviews*, 215, 1-170. https://doi.org/10.1007/s11214-018-0574-6

Lognonné, P., Banerdt, W. B., Pike, W. T., Giardini, D., Christensen, U., Garcia, R. F., et al. (2020). Constraints on the shallow elastic and anelastic structure of Mars from insight seismic data. *Nature Geoscience*, 13(3), 213-220. https://doi.org/10.1038/s41561-020-0536-y

Lognonné, P., Gagnepain-Beyneix, J., & Chenet, H. (2003). A new seismic model of the Moon: Implications for structure, thermal evolution and formation of the Moon. *Earth and Planetary Science Letters*, 211(1-2), 27-44. https://doi.org/10.1016/s0012-821x(03)00172-9

Lognonné, P., & Johnson, C. (2015). 10.03—Planetary seismology. *Treatise on Geophysics*, 2, 65-120.

Lognonné, P., Schmerr, N., Antonangeli, D., Bailey, S., Banerdt, B., Banks, M., et al. (2020). Seismology on Artemis III: Exploration and science goals. White paper.

Majstorović, J., Lognonne, P., Kawamura, T., & Panning, M. (2024). Predicting deep moonquake source regions using their temporal and spatial patterns and machine learning [Software]. Zenodo. https://doi.org/10.5281/zenodo.11103122

Matsumoto, K., Yamada, R., Kikuchi, F., Kamata, S., Ishihara, Y., Iwata, T., et al. (2015). Internal structure of the Moon inferred from Apollo seismic data and selenodetic data from GRAIL and LLR. *Geophysical Research Letters*, 42(18), 7351-7358. https://doi.org/10.1002/2015GL065335

Meeus, J. (1991). Astronomical algorithms.

Middlehurst, B. M. (1967). An analysis of lunar events. *Reviews of Geophysics*, 5(2), 173-189. https://doi.org/10.1029/rg005i002p00173

Minshull, T., & Goulty, N. (1988). The influence of tidal stresses on deep moonquake activity. *Physics of the Earth and Planetary Interiors*, 52(1), 41-55. https://doi.org/10.1016/0031-9201(88)90056-8

Moore, P. (1968). Transient lunar phenomena: A review. 1967. *The Journal of the British Astronomical Association*, 78, 138-144.

Nakamura, Y. (1983). Seismic velocity structure of the lunar mantle. *Journal of Geophysical Research*, 88(B1), 677-686. https://doi.org/10.1029/JB088iB01p00677

Nakamura, Y. (2003). New identification of deep moonquakes in the Apollo lunar seismic data. *Physics of the Earth and Planetary Interiors*, 139(3), 197-205. https://doi.org/10.1016/j.pepi.2003.07.017

Nakamura, Y. (2005). Farside deep moonquakes and deep interior of the Moon. *Journal of Geophysical Research*, 110(E1), E01001. https://doi.org/10.1029/2004JE002332

Nakamura, Y., Latham, G., Dorman, H., & Harris, J. (1981). Passive seismic experiment long-period event catalog. *Galveston Geophysics Laboratory Contribution*, 491, 314.

Nakamura, Y., Latham, G. V., & Dorman, H. J. (1982). Apollo lunar seismic experiment—Final summary. *Journal of Geophysical Research*, 87(S01), A117-A123. https://doi.org/10.1029/jb087is01p0a117

Nunn, C., Garcia, R. F., Nakamura, Y., Marusiak, A. G., Kawamura, T., Sun, D., et al. (2020a). Lunar seismology: A data and instrumentation review. *Space Science Reviews*, 216(5), 89. https://doi.org/10.1007/s11214-020-00709-3

Nunn, C., Garcia, R. F., Nakamura, Y., Marusiak, A. G., Kawamura, T., Sun, D., et al. (2020b). Lunar seismology: A data and instrumentation review—Electronic supplement [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.7787610

Panning, M., Kedar, S., Bowles, N., Calcutt, S., Cutler, J., Elliott, J., et al. (2021). Farside seismic suite (FSS): First seismic data from the farside of the Moon delivered by a commercial lander. In *AGU fall meeting abstracts* (Vol. 2021, p. P54C-01).

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., et al. (2011). Scikit-learn: Machine learning in python. *Journal of Machine Learning Research*, 12, 2825-2830.

Pickrell, J. (2022). These six countries are about to go to the moon—Here's why. *Nature*, 605(7909), 208-211. https://doi.org/10.1038/d41586-022-01252-7

Pou, L., Panning, M., Styczinski, M., Melwani Daswani, M., Nunn, C., & Vance, S. (2024). Tidal seismicity in the Moon and implications for the rocky interior of Europa. *Planetary and Space Science*, 5(6), 142. https://doi.org/10.3847/psj.ad47bc

Powers, D. M. (2020). Evaluation: From precision, recall and F-measure to ROC, informedness, markedness and correlation. arXiv preprint arXiv:2010.16061.

Probst, P., Wright, M. N., & Boulesteix, A.-L. (2019). Hyperparameters and tuning strategies for random forest. *Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery*, 9(3), e1301. https://doi.org/10.1002/widm.1301

Quinlan, J. R. (1986). Induction of decision trees. *Machine Learning*, 1, 81-106. https://doi.org/10.1007/bf00116251

Rhodes, B. (2019). Skyfield: Generate high precision research-grade positions for stars, planets, moons, and Earth satellites [Software]. *Astrophysics Source Code Library*. Retrieved from https://rhodesmill.org/skyfield/

Standley, I. M., Pike, W. T., Calcutt, S., & Hoffman, J. P. (2023). Short period seismometer for the lunar farside seismic suite mission. In *2023 IEEE aerospace conference* (pp. 1-9). https://doi.org/10.1109/AERO55745.2023.10115559

Toksöz, M. N., Dainty, A. M., Solomon, S. C., & Anderson, K. R. (1974). Structure of the Moon. *Reviews of Geophysics*, 12(4), 539-567. https://doi.org/10.1029/RG012i004p00539

Toksöz, M. N., Goins, N. R., & Cheng, C. (1977). Moonquakes: Mechanisms and relation to tidal stresses. *Science*, 196(4293), 979-981. https://doi.org/10.1126/science.196.4293.979

Turner, A. R., Hawthorne, J. C., & Gaddes, M. (2022). Stresses in the lunar interior: Insights from slip directions in the a01 deep moonquake nest. *Journal of Geophysical Research: Planets*, 127(12), e2022JE007364. https://doi.org/10.1029/2022JE007364