Evidence of shallow basaltic lava layers in Von Kármán crater from Yutu-2 Lunar Penetrating Radar

Yutu-2


Introduction
Ground-penetrating radar (GPR) is a mature methodology with a wide range of applications in geoscience (Schroeder et al., 2020), civil/environmental engineering (Sonkamble and Chandra, 2021) and landmine detection (Daniels, 2005).GPR has been used in planetary science for orbiter sounders (Orosei et al., 2018;Kaku et al., 2017), and in the last ten years has become a standard part of the scientific payload of planetary rovers (Lai et al., 2019;Giannakis et al., 2023).Currently, there are two active planetary rovers on Mars (Curiosity, Perseverance) and one on the Moon (Yutu-2).By 2028, three additional rovers are planned to be deployed, with one on Mars (Rosalind Franklin) and two on the Moon (Chang'E-7, VIPER).Four of the active and planned rovers (Yutu-2 (Ding et al., 2022), Perseverance (Eide et al., 2021), Rosalind Franklin (Hervé et al., 2020), Chang'E-7 (Zou et al., 2020)) are equipped with GPR as part of their scientific payload.Additionally, in-situ GPR was utilized in Yutu-1 (Su et al., 2014;Fa et al., 2015) and channel with ≈60 MHz central frequency (Li et al., 2020).The interpretation of the low frequency data is debatable because of the low quality of the data corrupted by clutter due to interactions between the un-shielded antenna and the metallic parts of the rover (Li et al., 2018;Zhang et al., 2020).Nonetheless, recent approaches have managed to derive useful information using the low frequency radar data, revealing a complex layered structure with multiple Imbrian basaltic layers (Feng et al., 2023).The data from the high-frequency channels (for both Yutu-1 and Yutu-2) are of sufficient quality providing a unique opportunity to conduct a detailed investigation of the Lunar shallow subsurface.
The high frequency GPR on board of Yutu-2 provided pivotal information regarding the stratigraphy of VK crater, revealing a complex layered structure for the first 30-40 m depth (Li et al., 2020;Zhang et al., 2020).Here we need to mention that such abnormal penetration depths (for terrestrial standards) are due to the absence of liquid water, which is the main contributing factor to electromagnetic losses and signal attenuation in GPR (Daniels, 2004).The absence of any visible horizons in the radagram within the first ≈ 200 ns lead to the early conclusion that the top ≈ 10 m of the landing site is the weathered Lunar regolith (Li et al., 2020), a fine-grained medium with electric permittivity that increases monotonically with depth (Dong et al., 2020).Using the estimated electric permittivity and semi-empirical formulas (Carrier et al., 1991;Hickson et al., 2018;Olhoeft and Strangway, 1975), the FeO + Ti 2 O 3 content and the density of the ejecta were approximated down to 50 m depth (Dong et al., 2020).These early conclusions were challenged by recent papers that underlined the uncertainty of permittivity estimation using conventional hyperbola fitting (Giannakis et al., 2022), and provided evidence for a more complex layered structure within the first 10 m depth (Giannakis et al., 2021(Giannakis et al., , 2023)).Further studies demonstrated a complex distribution of dielectric properties, with paleo-craters (Zhang et al., 2021) and multiple overlaying ejecta blankets (Feng et al., 2022a) with debatable origins (Giannakis et al., 2021).
In this paper, we infer a detailed dielectric structure for the first ≈ 30 m depth at the Chang'E-4 landing site.Using stochastic hyperbola fitting and probabilistic inversion (as described in Giannakis et al. (2023)), we firstly obtain a pseudo-2D permittivity profile along the first ≈ 1000 m travelled by the Yutu-2 rover and down to ≈ 10 m depth.A complex structure is inferred consistent with previous results, Giannakis et al. (2023) with a layer with high relative permittivity ( > 8) at approximately ≈ 10 m depth.Furthermore, mapping the central frequency of the received signal (Ding et al., 2020a) revealed a clear layered structure where the central frequency exhibits a lowrate decrease within the first 200 ns, followed by a rapid decrease from 200-400 ns, and a frequency plateau with no further decrease from 400-500 ns.Via detailed theoretical analysis, we demonstrate that conductivity cannot explain the frequency dispersive phenomena observed within the first 400 ns.Dispersive Cole-Cole (Cole and Cole, 1941) phenomena associated with ilmenite (Boivin et al., 2022a) are necessary to explain the observed shift of central frequency with depth.Based on the laboratory measurements presented in Boivin et al. (2022a), we support the premise that a varying ilmenite content within the first ≈ 30 m, and a high-ilmenite layer with more than 10% ilmenite content starting at ≈ 10 m depth, are necessary to explain the observed frequency dispersion.High-ilmenite content can exist in both hightitanium basalts (Meyer, 2012) and immature basaltic soils (Chambers et al., 1995), although according to Heiken and Vaniman (1990) the ilmenite content of Lunar soils is usually below 10%.Based on the significant ilmenite content inferred at ≈ 10-20 m depth, the estimated -from probabilistic inversion-high permittivity of this layer, distinctive morphological features, and the lack of visible rocks and boulders in the radargram ( although lack of hyperbolic signatures might be due to interferences between different point targets), we propose the presence of a shallow sequence of high-Ti basaltic layers part of the Imbrian basaltic flood.The Imbrian basaltic layers overlay a lowilmenite ejecta blanket, and are situated underneath a complex layered structure consisted of Eratosthenian ejecta.(Lu et al., 2021).The age of Leibnitz crater is constrained by the age of its floor (mare basalt) ≈ 3.3 Ga.Leibnitz crater is considered to be younger than Von Kármán crater but of similar age i.e.Nectarian/pre-Nectarian (Wang et al., 2023).

Geology of the landing site
The VK crater is within the Mg-rich annulus (Moriarty and Pieters, 2018) and its age was estimated Nectarian/pre-Nectarian (Huang et al., 2018;Losiak et al., 2009) at ≈ 4.2 Ga (Lu et al., 2021), very close to the formation of SPA (Hiesinger et al., 2012;Lu et al., 2021).The stratigraphy of the VK crater is speculated to consist of crater ejecta layers and multiple basaltic layers from the Imbrian basaltic floods (Huang et al., 2018;Pasckert et al., 2018;Lu et al., 2021;Gou et al., 2021).In particular, from morphological characteristics (see Fig. 1) it is evident that Leibnitz crater provided the north part of VK crater with ejecta layers prior to the Imbrian basaltic floods (Huang et al., 2018).Imbrian craters in the near proximity such as the Alder crater (dated at ≈ 3.5 Ga (Lu et al., 2021)) might also have contributed to the ejecta prior to the Imbrian basaltic layers (Huang et al., 2018).The basaltic floods are dated at ≈ 3.2-3.3Ga, and is believed that they took place prior to the Eratosthenian craters Finsen, VK L and VK L' (Losiak et al., 2009;Zhang et al., 2020;Paskert et al., 2018).Consequently, the top post-Imbrian layers in the VK crater are expected to be ejecta of the Eratosthenian craters Finsen, VK L and VK L' (Giannakis et al., 2021), while there are some theories that suggest that ejecta from Orientale basin might have contributed as well (Sun et al., 2021).Using LROC NAC images, the thickness of the weathered top soil is estimated at ≈ 2.5-7.5 m (Huang et al., 2018).Below the weathered top soil -based on the  3 reflectance data-it is speculated that there is a low-calcium pyroxene (LCP) layer (Gou et al., 2019;Li et al., 2019) ranging from ≈ 8-13 m, on top of a high-calcium pyroxene (HCP) layer from ≈ 13-53 m (Huang et al., 2018).

Radar dataset
The data used in this paper are the published lunar-penetrating radar (LPR) data from the second channel of the Yutu-2 rover (Feng et al., 2022a).The data were collected using a common-offset configuration, and a typical processing pipeline was applied to enhance signal and reduce unwanted clutter and ringing noise.In particular, a finite impulse response filter with range 250-750 MHz was initially applied (Feng et al., 2022a); subsequently a background removal was used to remove the direct coupling and ringing noise; and lastly, electromagnetic losses and geometrical spreading are mitigated via a quadratic gain.For illustration purposes, the envelope of the signal (Hilbert transform (Daniels, 2004)) can also be used to highlight reflectors and reveal hidden structures.Fig. 2 (A) illustrates the envelope of the processed radargram along the 1000 m trajectory of Yutu-2 rover.The quality of the signal is sufficient throughout the scan, and the signal is reliably above the noise level for the first ≈ 450 ns.
There are six distinct formations that can be identified from Fig. 2. Formation 1 is the top Lunar soil consisting of fine-grained materials with numerous buried rocks/boulders as indicated by the plethora of clear hyperbola visible in this segment (Feng et al., 2022a;Giannakis et al., 2023).Formation 2 is speculated to be a paleo-crater (Zhang et al., 2021) responsible for the observed discontinuity in the formations 3-5.Formations 3-5 are distinct layers with no visible hyperbolic reflections, indicating lack of buried rocks/boulders, although this might be due to interferences between various point targets within these formations.Lastly, formation 6 is another layered structure with no visible hyperbola, and higher frequency content compared to the previous layers.This is most-likely due to the reduction of the amplitude of the signal at the interface between the formations 5 and 6, which drops the signal below the noise level.
In the next sections we will explore the dielectric properties of these formations, their origins and mineralogical content.In Section 4 we will focus on formation 1, while in Section 5 we will explore the nature of formations 3-5.

Pseudo-2D permittivity inversion
We first focus on the first ≈ 200 ns of the processed radargram.We will use probabilistic inversion based on stochastic hyperbola fitting as presented in Giannakis et al. (2023).Conventional hyperbola fitting assumes an ideally spherical target, with a known radius (often assumed to be zero), buried in a homogeneous half-space (Giannakis et al., 2022).In stochastic hyperbola fitting (Giannakis et al., 2023), the radius of the target is considered unknown, and the permittivity of the medium varies smoothly with respect to depth (Giannakis et al., 2023).In contrast to conventional hyperbola fitting where a single set of bulk electric permittivity (  ) and depth () are fitted to a given set of arrival times (), stochastic hyperbola fitting takes into account the inherit non-uniqueness of the problem (Giannakis et al., 2022) resulting in a conditional probability  (  , |) for every hyperbola.This is very important because conventional hyperbola fitting has very high uncertainty that often spans along the permittivity range of interest in Lunar and planetary radar (Giannakis et al., 2023).This is often undermined in planetary radar, which results in biased and unreliable estimates (Giannakis et al., 2022(Giannakis et al., , 2023)).Using conventional hyperbola fitting assuming a zero-radius will systematically overestimate the bulk velocity (and consequently underestimate the bulk permittivity), especially in the presence of rocks/boulders with large radius (Giannakis et al., 2022).Moreover, using Dix conversion to transform bulk to actual permittivity will result to unreliable estimates, since Dix conversion was developed for seismic exploration in areas with distinct layers and not gradational complex formations such as the ones expected in a Lunar setup (Giannakis et al., 2021).Here we need to mention that the bulk permittivity is estimated from the bulk velocity, i.e. the term ''bulk permittivity'' used throughout this paper does not refer to the average permittivity of a multi-phase complex medium, but it refers to the equivalent permittivity of a homogeneous medium that would give rise to the same velocity observed in the investigated multi-phase complex medium.
We have manually picked 101 hyperbola along the 1000 m travelled by the Yutu-2 rover.Fig. 3 shows the distribution of the selected hyperbola along the processed radargram.Notice that all hyperbola are within formation 1, since no clear hyperbola were observed in the other formations.Fig. 4 shows the conditional probabilities  (  , |) for the 101 selected hyperbola.Despite the wide uncertainty range due to the inherit non-uniqueness of the problem (Giannakis et al., 2022(Giannakis et al., , 2023)), a trend ( shown with red line in Fig. 4) can be observed where the bulk permittivity increases down to ≈ 3-4 m depth, then decreases monotonically until ≈ 9 m depth, and then it mildly increases until ≈ 11 m.This trend was also observed by Feng et al. (2022b); and is consistent with the results shown in Giannakis et al. (2023) where stochastic hyperbola fitting was applied to the first 100 m travelled by Yutu-2 rover.The trend in Fig. 4 (red line) was calculated by fitting a 3rd order polynomial to the permittivity values with maximum conditional probability for each hyperbola (green markers).
Next, we use the 1D probabilistic inversion as described in Giannakis et al. ( 2023) to derive a pseudo-2D permittivity profile.The 1D probabilistic inversion uses the conditional probabilities  (  , |) of numerous hyperbola to infer the 1D permittivity structure of the investigated area.In this paper we use the 1D probabilistic inversion in overlapping sliding windows with 70 m width and 20 m step.The probabilistic inversion (Giannakis et al., 2023) has the ability to infer the uncertainty range of the inverted permittivity, and therefore provide indications of how reliable the results are.The reliability is related to the number, distribution and quality of the selected hyperbola.For example, if there is an area with lack of hyperbola the uncertainty will be high.Based on that, we kept only the results for which two times the standard deviation of the estimated permittivity is below 2. In other words, all the estimated relative permittivity values are within an uncertainty range of ±2.
The results are shown in Figs. 5 and 6.The estimated permittivity at the surface is  ≈ 2, which is in good agreement with recent papers that estimated the surface permittivity using reflection coefficients (Ding et al., 2021).The permittivity increases to  ≈ 6 at approximately 2-4 m depth, and then starts to decrease monotonically to  ≈ 3 at around 8-9 m depth, where it shows a relatively sudden increase to  ≈ 11 ± at approximately 11 m depth at the boundaries between formation and formations 2, 3 and 4.
The high permittivity layer at ≈ 2-4 m is in good agreement with the results in Giannakis et al. (2023), and also consistent with the trend of the bulk permittivity shown in Fig. 4. As numerically demonstrated in Giannakis et al. (2023), the lack of clear reflectionhorizons from these layers are due to the smooth variation of the permittivity that can decrease the overall reflection coefficients of an interface and lead to transparent structures that are not directly visible in the radargram (Diamanti et al., 2014;Giannakis et al., 2021Giannakis et al., , 2023)).
At ≈ 4 m depth, the permittivity decreases rapidly throughout the scan apart from the area around  ≈ 250 m.Notice that this is the area right on top of the speculated paleo-crater (Zhang et al., 2021).This indicates that the impact that created the paleo-crater at  ≈ 250 m, most-likely penetrated the low permittivity layer that starts at ≈ 4 m depth.The crater was later on filled with the same high permittivity ejecta that created the high permittivity layer at ≈ 2-4 m depth (see Fig. 7).
The uncertainty of the estimated permittivity at depths below ≈ 7 m does not allow us to reliably infer the permittivity throughout the interface between the formation 1 and the formations 2, 3 and 4 .This is due to the small number of hyperbola at depths greater than 140 ns, which increases the overall uncertainty of the results.Moreover, lack of shallow hyperbola can also increase the uncertainty of permittivity estimation at deeper depths.Hyperbola fitting estimates the conditional probability of the bulk/averaged permittivity down to a given depth.To infer the actual permittivity from its bulk permittivity at a given depth, we need to have a reliable knowledge of the permittivity down to that depth.Conclusively, a sufficient amount of hyperbola throughout the  3. Notice the high uncertainty due to the inherit non-uniqueness of hyperbola fitting when the radius of the body is unknown (Giannakis et al., 2022).Green colour markers depict the maximum probability for each hyperbola.In order to showcase the trend of the bulk permittivity, a 3rd order polynomial (red line) is fitted to the maximum probabilities (green markers).Above 3 m the distribution of the conditional probability is sparse and primarily concentrated in low permittivity values, while at ≈ 3-5 m depth the conditional probability spreads from   ≈ 3-8.Subsequently, we see a concentration of the probability at low permittivity values at ≈ 7-9 m depth followed by a mild increase at ≈ 10 m depth.2023), there is a layer with high permittivity at ≈ 2-3 m between two layers with low permittivity.Notice that there are indications of a very high permittivity layer at the boundary between the formation 1 and the formations 2, 3 and 4 (see Fig. 6).
profile needs to be present in order to have a reliable estimation of the permittivity down to a given depth.
There are only three areas (within the 70 m range of the overlapping windows where the inversion takes place) where the uncertainty in the estimated permittivity is within a reasonable range (±2).At  ≈ 250 and 700 m there are two areas with high permittivity, while at  ≈ 600 m is an area with low permittivity.Notice that the high permittivity area at ≈ 700 m is at the apex of the anticline of formation 4 (see Figs. 2 and 6), while the low permittivity area lays at the left limb of the anticline.As shown in Fig. 6, the two areas with high permittivity are at the boundaries between formation 1 and the formations 2, 3 and 4 , while the low permittivity area is within formation 1, and approximately ≈ 50 ns above formation 4. This indicates the existence of a thin high permittivity layer at the boundaries between the formation 1 and the formations 2, 3 and 4 , as illustrated in Fig. 7.
We speculate that this thin high permittivity layer might be the paleo-regolith of formations 3 and 4. Prior to the deposition of the ejecta of formation 1, micrometeorites and space weathering acting upon formations 3 and 4 led to a high permittivity paleo-regolith with  ≈ 11 ± 2 (see Fig. 7).The lack of visible reflections in the permittivity transitions shown in Fig. 6 is due to the gradational smooth transition between these layers leading to small reflection coefficients and negligible reflected signals (Diamanti et al., 2014).This phenomenon was discussed in Diamanti et al. ( 2014) and was also   (Giannakis et al., 2021;Diamanti et al., 2014).We speculate that formation 1.d is the high permittivity paleo-regolith of formations 3 and 4; and the formation 1.b is a high permittivity layer with higher rock abundance (Giannakis et al., 2023) compared to the low permittivity formations 1.a and 1.c.demonstrated numerically in Giannakis et al. (2021Giannakis et al. ( , 2023)).A smoothgradational transition is expected between different ejecta layers due to the reworking of the materials during the deposition of a new ejecta layer (Giannakis et al., 2021).
The high permittivity paleo-regolith implies that the permittivities of the formations 3 and 4 must be higher since space weathering leads to an increased porosity (Horz et al., 1991) and vitrification (Nash and Conel, 1973) that reduce the overall bulk permittivity of the host medium (Giannakis et al., 2021).A permittivity above  > 11 ± 2 can be explained by the presence of basalt.Based on samples collected in Apollo 11, 12 and 14, the relative permittivity of Lunar basalts range from  ≈ 9-15 (Chung et al., 1970(Chung et al., , 1972)), and the Lunar breccia from  ≈ 6-9 (Chung et al., 1970(Chung et al., , 1972)).The high permittivity of Lunar basalts (compared to terrestrial ones) is due to their high ilmenite content (Chung et al., 1970(Chung et al., , 1972)).Ilmenite is the most important ore of titanium, with very high permittivity  ≈ 35-80 (Chung et al., 1970;Parkhomenko, 1967), leading to the increased bulk permittivity of Lunar high-Ti basalts.
Fig. 7 shows a conceptual model for the first ≈ 300 ns of the radargram.The formation 1 is divided into four sub-sections 1.a, 1.b, 1.c and 1.d.We speculate that the formation 1.d is the high bulk permittivity paleo-regolith of formations 3 and 4. Formations 1.a and 1.c are low permittivity layers while formation 1.b has permittivity ≈ 5-6, most likely due to the higher rock abundance in that depth as shown in Giannakis et al. (2023).Fig. 7 implies that after the development of the paleo-regolith 1.d, an impact event deposited the ejecta 1.c.Then the paleo-crater at ≈ 200 m penetrated 1.c, 1.d and formation 3. The paleo-crater was filled with reworked materials from its rims and weathered materials from formation 3. Subsequently another impact event deposited the ejecta 1.b.Lastly, space weathering and micrometeorite showering developed the fine-grained low permittivity layer 1.a.
The high bulk permittivity of the formations 1.d and 1.b do not agree with the low permittivity values of Lunar soil-samples from the Apollo missions (Olhoeft and Strangway, 1975).This is due to the fact that Lunar soil samples (Olhoeft and Strangway, 1975) and Lunar analogues (Boivin et al., 2022a,b) consist of very fined grainedpowdered-materials with air-fraction > 50% that represent the shallow weathered top Lunar soil.From the cluttered radargram in the first 200 ns, and the migration results provided in Giannakis et al. (2023), it is evident that formation 1 is a rock-soil mixture with an increased rock abundance at ≈ 3-4 and ≈ 9-11 m depth (Giannakis et al., 2023).Deeper non-weathered ejecta are expected to consist of a mixture of soil and rocks/boulders that give rise to a bulk permittivity higher than the permittivity of powdered fine-grained soils.To illustrate that, we conducted the numerical experiment shown in Fig. 8.A homogeneous half-space with permittivity  = 4 represents a high-density Lunar soil as expected at depth ≈ 10 m.The medium is then filled with fragments of high-Ti basalt with permittivity  = 11 (Chung et al., 1970(Chung et al., , 1972)).This is a 2D model simulated using the finite-difference time domain (FDTD) method (Taflove and Hagness, 2000).The spatial step equals with  =  = 3 mm, and the time-step is 0.99 of the Courant limit (Taflove and Hagness, 2000).The boundaries are truncated using the time-synchronized convolutional perfectly matched layer (Giannakis and Giannopoulos, 2015).The source is an ideal Hertzian dipole with 500 MHz central frequency, and the receiver is placed at 1 m depth.From the received signal we can estimate the bulk velocity of the medium, and consequently its bulk permittivity.Fig. 8 shows the estimated bulk permittivity for different volumetric fractions of high-Ti basaltic fragments.The bulk permittivity reaches values above  > 7.5 when basaltic rock fragments occupy more than 30% of the medium; and  > 9 when basaltic fragments occupy 50% of the overall volume.Consequently, the high permittivity values derived from the probabilistic 1D inversion can be explained via a mixture of Lunar soil and high-Ti basaltic fragments.Notice that the basaltic fragments do not have to be significantly large to affect the overall bulk permittivity.Due to the resolution constrains from the frequency range of the antenna, small basaltic fragments will not result to clear hyperbolic reflections, but they will affect the overall propagation velocity (and consequently the bulk permittivity) as clearly demonstrated in the numerical experiment shown in Fig. 8.

Frequency dispersion
In the previous section, using 101 hyperbola from formation 1, we inferred the dielectric properties of the landing site down to ≈ 10 m I. Giannakis et al. Fig. 8.A numerical 2D experiment to evaluate the effect of high-Ti basaltic fragments to the bulk velocity (and consequently bulk permittivity) of the medium.The permittivity of the soil (grey colour) is  = 4, and the permittivity of the high-Ti basaltic fragments (white colour) is  = 11 (based on Lunar samples from the Apollo missions (Chung et al., 1970(Chung et al., , 1972))).Left diagram illustrates how the model looks for 10% volumetric fraction of basaltic fragments.Right diagram shows the bulk permittivity of the medium with respect to the volumetric fraction of basaltic fragments.The bulk permittivity is estimated from the bulk velocity of the medium calculated from the received signal at 1 m depth.The source is an ideal Hertzian dipole with 500 MHz central frequency.depth.The lack of hyperbola in formations 2-6, does not allow us to use probabilistic inversion to get a direct insight on the dielectric properties of these layers.Indications were given in the previous section that supported the premise that formations 3 and 4 might be high permittivity layers with permittivity  > 10, pointing to high-Ti Lunar basalts (according to the dielectric properties of Lunar samples from the Apollo missions (Chung et al., 1970;Parkhomenko, 1967)).Nonetheless, these indications need to be backed up with more robust and coherent evidence to further support the existence of shallow basaltic layers in the VK crater.
In this section we focus on the frequency content of the reflected signals, and in particular on the central frequency.We first calculate the Short-Time Fourier Transform (STFC) for each trace using a Hamming window with length 31.25 ns and step 0.3125 ns.The maximum frequency of the clutter for every distance () and time () is saved in a matrix (, ).Subsequently the matrix (, ) is smoothed using a Gaussian filter.Fig. 9 shows both the raw and smoothed (, ).From Fig. 9 it is apparent that there is a clear correlation between the central frequency of the received signal and the formations as defined in Fig. 2. In particular, as shown in Fig. 10 there are 3 areas with different rates of frequency shift aligned with the formations 1, 3/4 and 5 respectively.In the next three subsections we explore the source of these dispersive phenomena.We first showcase that a simple conductivity term cannot explain the shift of the central frequency.Subsequently, we demonstrate that the frequency-dependent electric permittivity of ilmenite (Boivin et al., 2022a) can explain the frequency shift, and based on that we try to get an insight into the ilmenite content of the formations 3 and 4.

Effect of conductivity on frequency shift
The propagation of a monochromatic electromagnetic plane wave can be expressed via (Balanis, 2012): where  is the electric field (V/m),   is the initial amplitude of the wave for frequency  (),  is the direction of the propagation (),  is time (s),  = 2 is the angular frequency, and  = √ −1 is the imaginary number.The variables  and  are related to the velocity and the losses of the wave respectively and they are equal to Balanis (2012): where  (∕ 2 ) and  ( ∕) are the magnetic permeability and electric permittivity respectively, and  (∕ ) is the electric conductivity.From Eq. ( 1) it is clear that the wave attenuates as it propagates in the  axis, and the attenuation of the wave is related to the term   ⋅  − .If we know the initial   for different frequencies, we can assess the affects of conductivity on the frequency content of the reflections.The only modification that we need to do for GPR is to account for the two way travel time of the wave i.e.   =   ⋅  −2 , where   is the frequency spectrum of the reflected wave from a target at depth .
The frequency spectrum of the initial pulse can be approximated with a Gaussian function with 500 MHz mean and 180 MHz standard deviation i.e. 500 MHz central frequency with 500 MHz bandwidth, similar to the frequency content of the second channel on board of Yutu-2 rover (Feng et al., 2022a).Based on that we calculated   for different dielectric properties.Even for extreme values of conductivity the central frequency of the reflected wave practically stays unaffected, indicating that conductivity cannot be responsible for the frequency shift observed in Fig. 9.This is also illustrated in Giannakis (2016), where the attenuation factor  (described in Eq. ( 2)) was plotted against frequency for various values of conductivity, indicating that the attenuation factor  is practically the same for frequencies > 250 MHz.

Ilmenite and Cole-Cole dispersion
Ilmenite is a titanium mineral that is abundant on the Moon, with higher concentrations in Lunar Maria (Prettyman et al., 2006;Jackson and Carter, 2007).Its high relative permittivity  ≈ 35-80 (Chung et al., 1970) makes it very important in planetary radar since it can greatly affect the overall dielectric properties of a medium and consequently the performance of GPR.Ilmenite is a linear dispersive material, meaning that its electric permittivity is complex and frequency dependent.Early evidence of this were given in Nelson et al. (1989), where the electric permittivity of pulverized ilmenite was calculated for frequencies 1-22 GHz exhibiting a clear dispersive behaviour.Recent papers (Boivin et al., 2022a,b) have measured the permittivity of powdered ilmenite and mixtures of it with powdered bytownite for frequencies between 100 MHz-8.5 GHz.According to the experimental work in Boivin et al. (2022a,b), ilmenite mixtures are Cole-Cole media (Cole and Cole, 1941) and their electric permittivity can be expressed via: where  is the angular frequency,  ∞ is the permittivity at infinite frequency,  is the difference between the static permittivity and  ∞ ,  0 is the relaxation time, and  is a positive constant.The Cole-Cole parameters of different mixtures of powdered ilmenite and powdered bytownite are given in Table 1.We can use these parameters and Eq. ( 1) to assess the effects of ilmenite mixtures on the frequency shift of a given pulse.Due to the fact that the permittivity is now complex, the frequency spectrum of a plane wave   after propagating in a Cole-Cole medium equals   =   ⋅  ⋅ℜ(−2−2) , where ℜ denotes the real part of a complex number.Notice that if the permittivity is dispersion-less and real, this equation reduces to the one used in Section 5.1.Fig. 11 shows the shift of the central frequency for the ilmenite mixtures shown in Table 1.Similar to Section 5.1, the frequency content of the initial pulse is a Gaussian function with 500 MHz central frequency and 500 MHz bandwidth.From Fig. 11 it is evident that ilmenite content is directly related to the frequency shift, and in particular, ilmenite content is proportional to the rate of frequency shift.

Ilmenite content
The average relative permittivity of formation 1, based on the results of probabilistic inversion (see Fig. 5) is approximately  ≈ 4.7.
Based on that, and as shown in Fig. 5, the two way travel time for the first 10 m is  ≈ 145 ns.This is aligned with the boundary shown in Fig. 10 between the two areas with different rate of frequency change.Therefore in the first 10 m of formation 1, the frequency shifts linearly from 500 MHz to ≈ 370 MHz.This is in good agreement with the bytownite-ilmenite powdered mixture with 10% ilmenite as shown in Fig. 11.Both the formation 1 and the 10% ilmenite mixture are low dielectric media with permittivity  ≈ 4, and both media result in a linear shift of the central frequency with respect to depth, both reaching a central frequency ≈ 350-370 MHz at depth 10 m.
Despite the fact that we do not have any direct information regarding the permittivity of formations 3 and 4, nonetheless, from Fig. 10 it is evident that the formations from  ≈ 200 − −350 ns have higher ilmenite content compared to formation 1, since the rate of frequency shift is clearly higher as shown in Figs. 9 and 10.From Fig. 10 is also evident that formation 5 exhibits no frequency shift, indicating the lack of ilmenite in that layer.Lastly, the increase of central frequency in formation 6, is most-likely due to the fact that the amplitude of the refracted waves from the interface of formations 5 and 6 falls below the noise threshold.

Formations 1 and 2
Via probabilistic inversion (using 101 hyperbola) and frequency attributes we derived a complex layered structure within formation 1 with bulk permittivity  ≈ 4.7 (see Fig. 5), and ilmenite content < 10% (see Section 5).Within formation 1 there are 4 detected layers (see Fig. 5) with smooth interfaces between them resulting in no direct reflections (Giannakis et al., 2021(Giannakis et al., , 2023;;Diamanti et al., 2014).The deepest layer of formation 1 is a high permittivity layer  > 9, most-likely the paleo-regolith of formations 3 and 4, developed via micrometeorite showering and space weathering acting upon formations 3 and 4. Overlaying this layer is a low permittivity formation followed by a rocky formation at approximately ≈ 2-4 m depth with relatively high bulk permittivity  ≈ 5-6.In Giannakis et al. (2023), clear indications were given that supported the premise that the rock abundance at ≈ 2-4 m is increased, which is probably the reason for the increased bulk permittivity of this layer (Giannakis et al., 2023) (further supported by the numerical experiment shown in Fig. 8).The top layer of formation 1 is a fine grained weathered layer with low permittivity  ≈ 2, which is in good agreement with independent measurements based on the reflections coefficients of the surface permittivity of the landing site (Ding et al., 2021).Fig. 7 shows a conceptual model for the first ≈ 300 ns of the radargram.The formation 1.d is the paleoregolith of the formations 3 and 4 formed prior to the deposition of the low permittivity ejecta 1.c.Then the paleo-crater at ≈ 200 m penetrated 1.c, Fig. 11.The shift of central frequency for different percentages of ilmenite for the ilmenite/bytownite mixtures shown in Table 1.The initial pulse has 500 MHz central frequency and 500 MHz bandwidth, similar to the pulse used in the high frequency radar on board of Yutu-2 rover (Feng et al., 2022a).
1.d and formation 3, creating formation 2 from the reworked materials that were deposited from the rims at the floor of the paleo-crater.Subsequently another impact event deposited the rocky ejecta 1.b.The top layer of 1.b was subjected to space weathering and micrometeorite showering resulting to the low permittivity top layer 1.a.

Formations 3 and 4
The inferred information regarding the formations 3 and 4 is outlined below: • There are no clear hyperbola within these formations (see Section 3), meaning that most-likely there are no distinct point targets (rocks/boulders) in these layers.Although lack of hyperbola might be due to interferences from multiple point targets.• The ilmenite content of formations 3 and 4 is high (significantly higher than 10%) based on the rate of frequency shift within these formations (see Section 5.3 and Fig. 9).• If layer 1.d (the last layer of formation 1) is the paleo-regolith of formations 3 and 4, then the high permittivity of 1.d indicates that the permittivity of formations 3 and 4 are also high and most-likely above  > 10 (see Section 3).This is also consistent with the estimated high ilmenite content of formations 3 and 4, since ilmenite is a mineral with a unique and characteristic high permittivity ( ≈ 30-80 (Chung et al., 1970;Parkhomenko, 1967)).
A high permittivity layered structure, with high ilmenite content and no visible rocks/boulders points to high-Ti basalts.This premise is supported by the fact that high-Ti basalts have permittivity values (based on samples from the Apollo missions)  ≈ 10 (Chung et al., 1970(Chung et al., , 1972)); ilmenite content ≈ 0-20% (Meyer, 2012); higher ilmenite content than Mare soil (according to Heiken and Vaniman (1990) the ilmenite content of Mare soils is < 10%, while for high-Ti basalts can go up to ≈ 20%); and a relatively homogeneous texture with no distinct dielectric contrasts from buried rocks/boulders as expected in ejecta blankets.

Formation 5
The only information that we can infer for formation 5 is the lack of dispersion (see Fig. 10), which according to Section 4 implies that the ilmenite content in this layer should be negligible.This also points to a relatively low permittivity (Chung et al., 1970;Parkhomenko, 1967), despite the high density (density is proportional to permittivity (Olhoeft and Strangway, 1975)) that is expected at that depth.The low ilmenite content indicates that the origin of the ejecta of formation 5 might be from the Lunar highlands.This is based on the fact that titanium oxides are almost entirely concentrated in the Mare and SPA regions, while Lunar highlands have negligible amounts (Prettyman et al., 2006;Jackson and Carter, 2007).A potential source of highlands' ejecta within the VK crater can be the ejecta from the Orientale basin, as speculated by Sun et al. (2021).The size of the Orientale ejecta in the VK crater was estimated at ≈ 8 m (Sun et al., 2021;Anon, 1974).Based on that, the two way travel time within formation 5 is ≈ 100 ns (see Fig. 2), therefore the relative permittivity of this layer is  ≈ 3.7, a low permittivity layer as expected from the absence of ilmenite.

Conclusions
Using stochastic hyperbola fitting, probabilistic inversion and frequency analysis we derive a detailed stratigraphic column for the landing site of the Chang'E-4 mission.Multi-source radar evidence are given to support the existence of a sequence of shallow basaltic lava layers starting at ≈ 10 m depth, overlaying an ejecta blanket with low ilmenite content.The current study utilizes frequency attributes to get an insight into the ilmenite content of the Lunar subsurface.This can potentially contribute to the detection of subsurface areas with high ilmenite resources that have been suggested as a potential source of oxygen for supporting future Lunar bases.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The landing site of the Chang'E-4 mission is illustrated with red cross in the Von Kármán crater at 44.45 • S, 176.3 • E. The dates of the depicted craters are based on(Lu et al., 2021).The age of Leibnitz crater is constrained by the age of its floor (mare basalt) ≈ 3.3 Ga.Leibnitz crater is considered to be younger than Von Kármán crater but of similar age i.e.Nectarian/pre-Nectarian(Wang et al., 2023).

Fig. 2 .
Fig. 2. (A) Processed radar data from the second channel of Yutu-2 rover.The processing pipeline consists of a band-pass filter, background removal, quadratic gain and Hilbert transform.(B) Conceptual sketch of the 6 distinct formations observed in the processed radagram.The paleo-crater defined as formation 2 was first introduced in Zhang et al. (2021).

Fig. 3 .
Fig. 3. (A) The selected 101 hyperbola used in the probabilistic inversion.Notice that all the selected hyperbola are within formation 1, since there are no clear hyperbola in the other formations.(B) 3 enlarged areas of (A) showing clearly visible hyperbolic targets present in formation 1.

Fig. 4 .
Fig.4.The conditional probabilities  (  , |) for all the 101 hyperbola shown in Fig.3.Notice the high uncertainty due to the inherit non-uniqueness of hyperbola fitting when the radius of the body is unknown(Giannakis et al., 2022).Green colour markers depict the maximum probability for each hyperbola.In order to showcase the trend of the bulk permittivity, a 3rd order polynomial (red line) is fitted to the maximum probabilities (green markers).Above 3 m the distribution of the conditional probability is sparse and primarily concentrated in low permittivity values, while at ≈ 3-5 m depth the conditional probability spreads from   ≈ 3-8.Subsequently, we see a concentration of the probability at low permittivity values at ≈ 7-9 m depth followed by a mild increase at ≈ 10 m depth.

Fig. 5 .
Fig. 5.The results of the pseudo-2D probabilistic inversion excluding the areas with high uncertainty range, where the two times standard deviation of the estimated relative permittivity is above 2. Consistent with the results presented in Giannakis et al. (2023), there is a layer with high permittivity at ≈ 2-3 m between two layers with low permittivity.Notice that there are indications of a very high permittivity layer at the boundary between the formation 1 and the formations 2, 3 and 4 (see Fig.6).

Fig. 6 .
Fig. 6.The results of the pseudo-2D probabilistic inversion shown in Fig. 6 overlaid on the processed radargram for the first ≈ 200 ns.Notice that the two areas with high permittivity at ≈ 150 ns are at the boundaries between formation 1 and the formations 2, 3 and 4, while the low permittivity area at the same depth is further away from the boundary due to the anticline structure of formation 4 at  ≈ 550-900 m.

Fig. 7 .
Fig. 7. Conceptual model for the first ≈ 300 ns of the radargram.The formation 1 as defined in Fig. 2 is sub-divided into four sections 1.a, 1.b, 1.c and 1.d.The boundaries between these layers are smooth and gradational leading to decreased reflection coefficients and the absence of clear horizons in the radargram(Giannakis et al., 2021;Diamanti et al., 2014).We speculate that formation 1.d is the high permittivity paleo-regolith of formations 3 and 4; and the formation 1.b is a high permittivity layer with higher rock abundance(Giannakis et al., 2023) compared to the low permittivity formations 1.a and 1.c.

Fig. 9 .
Fig. 9. (A) The central frequency of the received signal (, ).(B) Smoothed (, ).Notice that there are some clearly visible structures correlated with the formations in Fig. 2. The formations observed in the radagram (Fig. 2) seems to be aligned with the interfaces where the rate of frequency shift changes i.e. the green areas in Fig. 9 (B).

Fig. 10 .
Fig. 10.The smoothed distribution of the central frequency (, ) averaged over (A) x=0-250 m and (B) x=600-850 m.Notice that the three areas with different rate of frequency shift -illustrated with red, green and blue colour-are aligned with the formations 1, 3&4 and 5 respectively.