A potential subsurface cavity in the continuous ejecta deposits of the Ziwei crater discovered by the Chang’E-3 mission

In the radargram obtained by the high-frequency lunar penetrating radar onboard the Chang’E-3 mission, we notice a potential subsurface cavity that has a smaller permittivity compared to the surrounding materials. The two-way travel time between the top and bottom boundaries of the potential cavity is ~ 21 ns, and the entire zone is located within the continuous ejecta deposits of the Ziwei crater, which generally have similar physical properties to typical lunar regolith. We carried out numerical simulations for electromagnetic wave propagation to investigate the nature of this low-permittivity zone. Assuming different shapes for this zone, a comprehensive comparison between our model results and the observed radargram suggests that the roof of this zone is convex and slightly inclined to the south. Modeling subsurface materials with different relative permittivities suggests that the low-permittivity zone is most likely formed due to a subsurface cavity. The maximum vertical dimension of this potential cavity is ~ 3.1 m. While the continuous ejecta deposits of Ziwei crater are largely composed of pre-impact regolith, competent mare basalts were also excavated, which is evident by the abundant meter-scale boulders on the wall and rim of Ziwei crater. We infer that the subsurface cavity is supported by excavated large boulders, which were stacked during the energetic emplacement of the continuous ejecta deposits. However, the exact geometry of this cavity (e.g., the width) cannot be constrained using the single two-dimensional radar profile. This discovery indicates that large voids formed during the emplacement of impact ejecta should be abundant on the Moon, which contributes to the high bulk porosity of the lunar shallow crust, as discovered by the GRAIL mission. Our results further suggest that ground penetrating radar is capable of detecting and deciphering subsurface cavities such as lava tubes, which can be applied in future lunar and deep space explorations. Subsurface cavities of the Moon are the ideal natural shelter for future lunar explorations since they provide ideal insulation from cosmic radiation, meteoroid impacts, temperature fluctuations, etc. In this paper, we first report that a potential subsurface cavity with a height of ~3.1 m is revealed by the lunar penetrating radar onboard Chang’E-3 Yutu Rover Subsurface cavities of the Moon are the ideal natural shelter for future lunar explorations since they provide ideal insulation from cosmic radiation, meteoroid impacts, temperature fluctuations, etc. In this paper, we first report that a potential subsurface cavity with a height of ~3.1 m is revealed by the lunar penetrating radar onboard Chang’E-3 Yutu Rover Subsurface cavities of the Moon are the ideal natural shelter for future lunar explorations since they provide ideal insulation from cosmic radiation, meteoroid impacts, temperature fluctuations, etc. In this paper, we first report that a potential subsurface cavity with a height of ~3.1 m is revealed by the lunar penetrating radar onboard Chang’E-3 Yutu Rover


Introduction
In 2013, the Chinese Chang'E-3 mission (CE-3) carried out high-resolution ground penetrating radar profiling on the Moon and it was the first in situ application of ground penetrating radar on an extraterrestrial body. The landing site and the traverse route of the Yutu rover were both located on the continuous ejecta deposits of the ~ 450-m-diameter Ziwei crater ( Fig. 1a and b).
Using the high-frequency lunar penetrating radar (LPR) data returned along the route of the Yutu rover, subsurface structures were revealed in unprecedented detail (Fa et al. 2015;Feng et al. 2017;Xiao et al. 2015;Zhang et al. 2015;Zhang et al. 2019). Widespread disagreements exist about previous interpretations of possible subsurface structures in the radargram (Fa et al. 2015). From the perspective of echo strengths, there is a general consensus that the top ~ 50 ns of the radargram was returned from the continuous ejecta deposits of Ziwei crater. This interpretation is consistent with the geological context of the landing site (Qiao et al. 2016).
Within the top ~ 50 ns of the radargram, we noticed a small zone where the reflections from different depths feature phase reversals, which were not mentioned before. Phase reversal is a fundamental feature of electromagnetic wave propagation when the wave encounters an interface with materials that have a higher refractive index (Neal 2004). With ground-based validations, phase reveals of radar reflections have been successfully used to identify subsurface cavities, such as lava tubes on Earth (Miyamoto et al. 2005).
Using numerical simulation of electromagnetic wave propagation in media, this study is focused on unlocking the nature of the observed phase reversal zone in the radargram obtained by the Yutu rover. We confirm that the relative permittivity of this zone is most consistent with that of a cavity that is ~ 3.1 m high. In addition to the widely mentioned subsurface lava tubes on the Moon, the existence of impact-induced cavities is indicative of both the lunar crustal evolution and engineering needs of future deep space explorations.

Data
The high-frequency LPR onboard the Yutu rover Su et al. 2014) obtained 2308 traces of valid data from navigation points N101 to N208, a travel distance of ~ 107.4 m (Feng et al. 2017). The routinely processed radargram is used as our observational data (Fig. 1c), and the detailed processing procedure has been repeated in multiple previous studies (Dong et al. 2017;Feng et al. 2017;Ding et al. 2020); thus, it is not repeated here.
Phase reversals in reflections of ground penetrating radar (GPR) can be applied to identify subsurface cavities (Miyamoto et al. 2005;Park et al. 2018). When the radar wave is propagated from either a low-to-high and then highto-low permittivity media or from a high-to-low and then low-to-high permittivity media, the phase of radar echoes before and after entering this media is reversed in the A-scan (i.e., the relationship between the amplitude and two-way travel time of the received signal). For the definition of phase reversals that are used to detect cavities in the subsurface (e.g., Miyamoto et al. 2005;Park et al. 2018), phase reversals are for echoes that are retuned directly from the upper and lower bounds of the subsurface object, and the echo polarity at the top of the subsurface object should have the same polarity as the incident pulse. The incident pulse of the high-frequency LPR has a positive polarity , so phase reversals in the radargram that feature positive and negative echo polarities at the upper and lower reflected echoes, respectively, would indicate the existence of subsurface low-permittivity zones.
For known subsurface cavities on Earth, phase reversals of radar reflections are manifested in the A-scan (Miyamoto et al. 2005). Using GPR to detect subsurface cavities that have no surface expressions is mainly based on the pattern of reflections shown in radargrams (Park et al. 2018), which occur as parallel positive and negative signals at different depths (see Figs. 6 and 7 of Miyamoto et al. 2003). On the Moon, the existence and exact locations of subsurface cavities cannot be deduced from surface observations unless collapsed skylights (Haruyama et al. 2009) or regolith leakage pits (Shoemaker and Morris 1970) are observed on the surface.
A zone of phase reversals formed by subsurface low-permittivity materials is observed in the radargram obtained by the high-frequency LPR onboard the Yutu rover. During the initial search, four criteria were used to identify possible phase reversals in the radargram: (1) reflections from the upper and lower bounds of the potential zone of phase reversals are parallel with each other; (2) the distances between the upper and lower bounds are larger than the resolution of the high-frequency LPR (i.e., ~ 3.5 ns in the two-way travel time; Su et al. 2014); (3) the echo polarity at the top of the candidate zone has the same polarity as the incident pulse; and (4) the phases at the upper and lower bounds of the A-scan are reversed. In the radargram, an ~ 4-m-wide zone between navigation points N106 and N108 is the only location that is qualified as a zone of phase reversal (Fig. 1c). The observed phase reversals are recorded in the A-scan ( Fig. 1f-h).

Method
In this study, we carried out numerical simulations of the propagation of electromagnetic waves in different media to (1) confirm the existence of phases revealed in the radar data and prove that the phase reversals are caused by a pocket of low-permittivity materials; (2) study the possible geometry of this low-permittivity zone, and (3) deduce the possible range of permittivity for this zone so that its nature can be constrained. The 2D Finite-Difference Time-Domain (FDTD) code for ground penetrating radar is used for the simulations (Irving and Knight 2006). The FDTD is a widely used code that performs differential calculations of Maxwell equations during the propagation of point-source electromagnetic waves in media.
The simulated electromagnetic wave in our models is identical to that transmitted by the CE-3 high-frequency LPR, i.e., a 500-MHz bandwidth and a quasi-Ricker wavelet . For the relative permittivity of lunar materials used in our model, the background dielectric constant is assumed to increase with depth as constrained by Apollo core samples, and the background dielectric constant is a function of the bulk density of lunar regolith (Carrier et al. 1991). For the loss tangent of materials used in our model, we employ the empirical relationship between the TiO 2 + FeO content and loss tangent (Carrier et al. 1991). The TiO 2 + FeO content at the Chang'E-3 landing site is constrained to ~ 25.5% , and the loss tangent is assumed to be ~ 0.0184. Subsurface structures that have different geometries and dielectric constants are used as our model inputs.

Confirmation of phase reversals in the radargram
The zone of phase reversals in the LPR radargram ( Fig. 1d) is found based on the reversely occurring parallel positive and negative reflections (Fig. 1e). The upper hyperbolic echo pattern exhibits more or less symmetric wings, while the summit is slightly inclined to the left (Fig. 1e). The echo polarity at the upper hyperbolic pattern is positive, which is consistent with the polarity of the incident pulse ). Compared to distinct phase reversals formed by subsurface low-permittivity zones (Miyamoto et al. 2005), those observed in the LPR A-scan ( Fig. 1f-h) are not as obvious in terms of the magnitudes of the positive and negative signals. The different radar systems employed are a possible reason for the unobvious signal magnitudes, i.e., a low-frequency stepped-frequency radar system by Miyamoto et al. (2005) versus a high-frequency time-domain radar system of CE-3 . This is because stepped-frequency radar can substantially reduce the strength of both direct and coupled waves (Hamran et al. 1995), so more power of the transmitted radar waves can be delivered to the subsurface, causing larger system dynamic ranges (Iizuka et al. 1984). On the other hand, the different subsurface structures are another possible reason for the different signal magnitudes observed, i.e., the lava tube studied by Miyamoto et al. (2005) and regolith interior structures identified by CE-3. Notably, Miyamoto et al. (2003) and Miyamoto et al. (2005) used two different GPR systems to detect a given lava tube, and the observed phase reversals exhibit different strengths in the radargrams. Indeed, for the known lava tube studied by Miyamoto et al. (2005), the A-scan does not always exhibit obvious phase reversals in terms of signal magnitudes, e.g., compare Fig. 7 in Miyamoto et al. (2003) and Fig. 3 in Miyamoto et al. (2005).
To confirm that the somewhat atypical phase reversals observed in Fig. 1f-h are real, we build a regolith model that contains a subsurface cavity. The model setup is not designed to reconstruct the radargram observed (Fig. 1d), but rather to show the characteristics of the phases recorded by the CE-3 high-frequency LPR. Figure 2 shows the model setup, the simulated radargram, and the extracted A-scans. The regolith models are 10 × 6.5 m, and the space interval is set to 0.01 m. The time window is set to 100 ns, and the polarity of the incident pulse is set to positive. The cavity (relative permittivity ε = 1) is (See figure on next page.) Fig. 1 The Chang'E-3 landing site and a zone of phase reversals exist in the top 50 ns of the radargram. a Image returned by the Lunar Reconnaissance Orbit (M102285549LE + RE; 1.66 m/pixel) showing that the Chang'E-3 landing site (red star) and the traverse route of the Yutu rover are both on the continuous ejecta deposits of the Ziwei crater. b Traverse route of the Yutu Rover from navigation points N101 to N208 (black dots). Navigation points that start with '1' represent the first lunar day, while those beginning with '2' represent the second lunar day. The red star represents the position of the lander. The segment between N106 and N108 is marked in green, which represents the surface projection of the observed low-permittivity zone. The background image was acquired by the descent camera of the Change'-3 Lander (ID: CE3_BMYK_LCAM-3006). c Radargram obtained by the high-frequency LPR channel. The radargram is processed following a standard routine outlined by Feng et al. (2017). The left y-axis is the apparent depth that corresponds to an assumed dielectric constant of 1, and the two-way travel time of the electromagnetic wave is shown in the right y-axis. The surface projection of the low-permittivity zone is denoted as green lines. d Regional radargram for the zone of phase reversals of reflections. The horizontal distance is identical to those used in panel c, and the y-axis is the two-way travel time. e Radargram with annotations of phase reversal locations. The negative peak of the reflected signal is marked in red, and the positive peak is marked in green. The three yellow dashed lines f-h correspond to traces of data where the A-scan is extracted. Red and green arrows point to the negative and positive peaks, respectively, where the phase reversals are identified implanted at a depth of ~ 3.2 m, and the shape is generated using the Gaussian stochastic equation.
For the regolith model that contains a subsurface cavity (Fig. 2a), phase reversal is not as obvious in the A-scan in terms of the echo strengths ( Fig. 2d and e). However, phase reversals are well manifested in the radargram, as the reversely occurring positive and negative signals are roughly parallel with each other (Fig. 2b and c). Similar atypical phase reversals in the A-scan have been observed for regolith models that contain subsurface cavities with different geometries (Additional file 1: see Text S1 and Figure S1). In the simulation results, the echo polarity at the top of the cavity is positive ( Fig. 2d and e), which is consistent with that of the simulated incident pulse. In conclusion, it is highly likely that the observed phase reversals in the CE-3 LPR radargram are true (Fig. 1d), and they should be caused by a pocket of subsurface materials that have a lower permittivity compared to the surrounding materials.

Roof shape of the low-permittivity zone
The roof shape of this low-permittivity zone can be roughly deduced by comparing the observed radargram with simulations for subsurface cavities that have different shapes (Fig. 3). Note that in this step, we are not certain that the low-permittivity zone is indeed a cavity, but all the models assume a subsurface cavity instead of other low-permittivity materials (e.g., more porous regolith). Using subsurface cavities in the regolith models only qualitatively constrains the roof shape of the low-permittivity zone. Large cavities are not mechanically stable in shallow porous materials unless heterogeneous materials in terms of strength exist around the cavity. Therefore, the regolith models assume that the subsurface cavities are supported by platy basaltic boulders that have a relative permittivity of 7.5 . Based on more than 150 models that consider subsurface cavities with different shapes and/or geometries (Additional file 1: see Text S2 and Figures S2-S11), it can be ascertained that the geometry of the upper bound of phase reversals (i.e., the upper hyperbolic echo pattern shown in Fig. 1d) is mainly controlled by the roof shape of the cavities. The comparison between the model results and the observations suggests that a slightly convex roof is the best approximation for this low-permittivity zone (Fig. 3f ). The criteria used to determine the similarity between the upper hyperbolic echoes in the radargram (Fig. 1d) and simulated reflections caused by the roof include the following: (1) the hyperbolic echo pattern exhibits more or less symmetric wings, and (2) the summit of the hyperbolic echo pattern is slightly inclined to the left. We must mention that this quantitative morphologic comparison cannot constrain the detailed shape of the roof, since possible combinations of various subsurface structures (e.g., the unknown distribution of small scatters in regolith such as multilayered structures) cannot be enumerated. The actual subsurface morphology and geometry of the low-permittivity zone are definitely much more complicated than the simplified models have sketched ( Fig. 3; see also Additional file 1: Text S3 and Figure S12). However, as a solution to recover the upper hyperbolic echo pattern, our further modeling work is built based on the best approximation shown in Fig. 3f.

Dielectric constant of the low-permittivity zone
Based on the best-fit roof geometry of the low-permittivity zone (Fig. 3f ), different relative permittivities (ε) are assigned for materials in this zone to deduce the possible range of ε. Since the absolute magnitude of reflections cannot be directly compared between the observation and simulations, for each trace of the raw A-scan, we normalize the echo strength at the lower hyperbolic echo pattern against the maximum echo strength of the entire A-scan. The normalized echo strengths of reflections are extracted along the lower hyperbolic echo patterns (i.e., B-scan), and we compare the general trend of the B-scan between the observation and simulations. The logic supporting this approach is that assuming materials in the low-permittivity zone have different values of ε compared to the maximum echo strength along a given A-scan, the relative echo decay at the lower hyperbolic echo pattern shows the energy decay pattern after the electromagnetic wave has propagated through the media, which is only related to the relative permittivity used in our model (see "Method" section for an introduction of the assumed constant loss tangent). Therefore, for the entire low-permittivity zone, the normalized echo strength along the lower hyperbolic echo pattern can be used to deduce the possible range of dielectric constants. For the model setups, the observed phase reversals have a time-delay difference of ~ 21 ns (see Additional file 1: Text S4 and Figure S13), and the height of the subsurface low-permittivity zone is varied according to the assumed dielectric constant (Fig. 4). It has been constrained that materials within the top ~ 50 ns of the LPR data have similar physical properties with typical regolith (Fa et al. 2015), so that the low-permittivity zone should feature a dielectric constant less than that of typical regolith (ε ≈ 3) (Carrier et al. 1991). The dielectric constant of the low-permittivity zone is set to 1-3 with an interval of 0.5, i.e., the 1st to 5th models shown in Fig. 4. Similar to the models shown in Fig. 3, the subsurface low-permittivity zones in our models are assumed to be supported by platy basaltic blocks. This assumption is supported by the fact that Ziwei crater was formed on the Eratosthenian-aged mare surface (Qiao et al. 2016). The width of the subsurface cavity is a free parameter in our models because although the general roof shape can be constrained, the exact geometry of this zone cannot be deciphered in greater detail using the single radar profile. In principle, the width of this zone does not affect the comparison of the normalized echo strength along the lower hyperbolic echo pattern. This is because the shape or geometry of the lower hyperbolic echo pattern is not used to constrain which assumed ε is correct, but the normalized echo strengths along the lower hyperbolic echo pattern are the evaluator. Indeed, because the low-permittivity zone has different heights and ε values, the lower hyperbolic echo patterns in our models exhibit different widths and geometries compared with the observation (denoted as red curves shown in Fig. 4g-k). On the other hand, up to 5 layers were observed in the  Comparison between the observed bottom radargram and those simulated for subsurface low-permittivity zones that have different dielectric constants. The regolith models use the depth profile of the dielectric constant that is adapted from the Apollo core samples, and the geometry of the regolith models is introduced in "Method" section. The implanted low-permittivity zone is assumed to be supported by basaltic plates (relative permittivity of 7.5). The thickness of the basaltic plates is assumed to be 0.2 m. The roof shape is adapted from the modeling results shown in Fig. 3f. The heights of the low-permittivity zones vary according to the assumed relative permittivity (i.e., denoted as "per") among the models, which ranges from 1-3 with an interval of 0.5. The width of the low-permittivity zone is not constrained (e.g., set to 2 m here), but this parameter does not affect our purpose ("Dielectric constant of the low-permittivity zone") section. Multiple-layered structures at a shallower regolith (Fa et al. 2015) are implanted, which have an assumed relative permittivity of 5. The first row shows the 5 model setups from the 1st to the 5th, the second row lists the observed and simulated radargrams, and the third row lists the annotated lower hyperbolic echo pattern (red curves) top ~ 15 ns of the radargram ( Fig. 1c; Fa et al. 2015). Multilayered structures within lunar regolith are likely caused by discontinuous materials that have a slightly larger bulk density than typical lunar regolith but a lower bulk density than competent rocks. Accounting for the multilayered structures observed in the top ~ 15 ns of the radargram (Fig. 1d), we insert similar structures in the regolith models, where the dielectric constant of the layers is set to 5, and the thickness of each layer is assumed to be 4 cm. Figure 4 shows the model setups. Compared to earlier simplified models (Figs. 2 and 3), the simulated radargrams shown in Fig. 4 are more similar to the observations. Again, for a fixed roof shape of the subsurface cavity, the upper hyperbolic echo pattern has a more or less constant shape regardless of the different cavity heights (Fig. 4a-e), width and existence of shallower structures (Fig. 3f ).
Regarding the distribution pattern of the normalized echo strengths along the lower hyperbolic echo pattern, the observation (black curve in Fig. 5) is most consistent with the model result that assumed ε = 1 for the low-permittivity zone (red curve in Fig. 5). A Kolmogorov-Smirnov test (Lilliefors 1967) is employed to compare the general trend of the normalized echo strengths along the lower hyperbolic echo pattern between the observation and each of our simulations. The Kolmogorov-Smirnov test is suitable for this purpose because this nonparametric statistical method does not need a prior distribution of the target (Lilliefors 1967;Razali et al. 2011). The null hypothesis used for our test is acceptable if the p-value ≥ 0.05. Table 1 shows the results of the Kolmogorov-Smirnov test. Only the first model shown in Fig. 4 fulfills the null hypothesis (p-value = 0.11). Therefore, our results support that the low-permittivity zone is most likely a cavity. It is also notable that the normalized amplitudes derived from the first model are not identical to those observed (e.g., from distances ~ 3.4 to 4 m), suggesting that differences exist between the constructed model and the actual subsurface structure.
Using ε = 1 as the relative permittivity for the low-permittivity zone, the ~ 21 ns difference in the two-way travel time between the upper and bottom phase reversals corresponds to a height of ~ 3.1 m for the cavity. The relationship between the depth of signal propagation (H) and the dielectric constant (ε) is H = c × t × ε −0.5 , where c is the speed of electromagnetic wave propagation in a vacuum and t is the two-way travel time (Ulaby et al. 1981). The top of this subsurface cavity is ~ 1.9 m below the surface. The result suggests that within the ~ 4.2-m-thick continuous ejecta deposits of Ziwei crater (Xiao et al. 2015), there is an ~ 3.1-m-tall cavity.

Flaws of the models
Flaws exist in our models, and they should be noted to avoid overinterpretation of the results. The demonstration that the observed phase reversals in the radargram are likely caused by a subsurface cavity is based on both geometrically simplified regolith models and various assumed relative permittivities. The background relative permittivity within the top 50 ns at the landing site was estimated to be ~ 3 on average (Fa et al. 2015), which is similar to that of typical lunar regolith (Carrier et al. 1991). However, if this assumption is invalid, the estimated relative permittivity for this low-permittivity zone might be slightly changed. Nevertheless, in the regolithlike ejecta deposits of Ziwei crater (Fa et al. 2015), a cavity is the most likely interpretation for this low-permittivity zone, considering that the ε = 1 model yields the most consistent fit with the observation (Fig. 5).
Basaltic plates are used to construct the low-permittivity zone (Figs. 3, 4), but we are not certain that the lowpermittivity zone must be supported by basaltic plates that have regular shapes rather than other less competent  Fig. 4. The calculation method of the normalized echo strength is introduced in "Dielectric constant of the low-permittivity zone" section. The black, red, green, blue, yellow and cyan curves are the extracted normalized echo strengths for the observation (Fig. 4l) and the 1st to 5th models (Fig. 4), respectively materials, e.g., discontinuous melt-rich breccia. Phase reversals are formed when the electromagnetic waves are propagated from high-to-low and then from low-to-high permittivity materials (Neal 2004), so changing the supporting materials of the low-permittivity zone would not affect the existence of the observed phase reversals, e.g., the simulations shown in Fig. 2 and Additional file 1: Text S1 of the supporting information.
The detailed architecture and geometry of the low-permittivity zone cannot be deciphered based on the simplified models shown in Fig. 3 and the Additional file 1. Assuming that the low-permittivity zone is supported by regularly shaped basaltic plates, the general shape of the roof can be estimated for the low-permittivity zone. However, the detailed shape (e.g., width) of the cavity along the route of Yutu cannot be constructed using the single radar profile. Furthermore, it should not be excluded that for a subsurface cavity in regolith, other debris materials might exist within the walls and floors of the cavity, which could cause further scattering of the radar waves. For example, using the observed plane morphology of ejected boulders around fresh impact craters (Krishna and Kumar 2016), we construct a more realistic regolith model that contains small rock piles on the floor of the cavity (Fig. 6a). The irregularly shaped small blocks are generated using a Gaussian random function, and their sizes are also randomly set within a range of 0.1-0.3 m. Indeed, the simulated radargram is more similar to the observation in terms of the shape of both the upper and lower hyperbolic echo patterns (Fig. 6b). Nevertheless, this model may still be simplified in comparison to the actual subsurface structure.

Origin of the subsurface cavity
The potential subsurface cavity was formed during the emplacement of the continuous ejecta deposits of the ~ 450-m-diameter Ziwei crater. This cavity is located ~ 50 m from the southeast-east rim of Ziwei crater (Fig. 1a). The thickness of ejecta deposits at this location is ~ 4 m, as estimated using an empirical function (Fa et al. 2015;Qiao et al. 2016;Xiao et al. 2015), which is consistent with most interpretations for the radargram (Fa et al. 2015;Xiao et al. 2015).
The subsurface cavity was likely formed due to the stacking of ejected basaltic boulders. Referring to the empirical relationship that the final crater diameter of simple craters is ~ 1.3 times that of the transient crater, we can infer that materials excavated from the Ziwei crater had a maximum original depth of ~ 11.5 m (Sharpton 2014). For comparison, the pre-impact target of Ziwei contained a regolith layer and competent basalts. This paleoregolith layer was mainly developed in the Eratosthenian-aged mare basalts. Based on the morphology and reflectance spectra of small impact Comparison between the observed radargram and that derived from a more realistic regolith model that contains a subsurface cavity. a A complex yet still simplified subsurface cavity model in lunar regolith. Blocks supporting the subsurface cavity have shapes adapted from ejected blocks around fresh impact craters (Krishna and Kumar 2016). The blocks have an assigned relative permittivity of 7.5 and thicknesses of ~ 0.35 m. The width of the cavity is assumed to be ~ 4 m here, and its height is set to ~ 3 m. Small irregularly shaped rocks are added at the bottom of the cavity. Multilayered structures are implanted within the shallow regolith, as shown in Fig. 4. The background relative permittivity is cited from the empirical equation of the Apollo core samples (Carrier et al. 1991). b Simulated radargram based on the model shown in panel a. c Annotated phase reversals in the simulated radargram, and the positive and negative amplitudes of the echo patterns are shown in green and red, respectively. d and e are the observed radargram and annotated phase reversals, respectively Ding et al. Earth, Planets and Space (2021) 73:53 craters, previous studies have suggested that the median thickness of this paleoregolith layer is ~ 8 m (Fa et al. 2014) and that the Eratosthenian-aged mare basalts are constrained to 40-43 m thick (Qiao et al. 2016). While the thickness of the regolith layer at the landing site shows large variations, the median thickness of ~ 8 m (Qiao et al. 2016) is considered the thickness of the regolith layer at the pre-impact site of Ziwei crater. Therefore, the ejecta deposits at the landing site were mainly from this regolith layer. This is consistent with the similar relative permittivity of subsurface materials within the top ~ 50 ns of the radargram with that of typical lunar regolith (Fa et al. 2015). On the other hand, the impact excavation model predicts that deeper materials are deposited closer to the crater rim, which are from larger depths (Melosh 1989). Therefore, as the excavation depth of Ziwei crater is larger than the depth of the pre-impact regolith layer, numerous meter-scale large blocks are visible on the wall and rim of Ziwei crater ( Fig. 7; Xiao 2014). It is highly likely that the subsurface cavity is supported by large basaltic blocks that were stacked during energetic emplacement.

Indications of the evolution of lunar shallow crust
The discovery of subsurface cavities caused by ejecta emplacement partly explains the high bulk porosity of the lunar shallow crust. High-resolution gravity data reveal that the bulk porosity of the lunar anorthositic crust is ∼12% (Wieczorek et al. 2013), and the upper 20 km of lunar crust features an even higher porosity. Ishiyama et al. (2013) reported high porosity within a depth of several hundred meters in the lunar mare and discussed the contributions of macro-cracks produced by meteorite impacts. It has been inferred that for the lunar uppermost crust, a portion of the porosity should occur as large-scale voids that cannot be sampled in hand specimen-sized rocks (Besserer et al. 2014). On the other hand, meter-scale boulders are common ejecta around fresh impact craters (Bart and Melosh 2010), so that large cavities formed during ejecta emplacement should not be a special case for the ~ 450-m-diameter Ziwei crater. Therefore, macro-cavities formed during the emplacement of block-containing ejecta could have contributed to the high bulk porosity of the lunar shallow crust. As the first indirect evidence for impact-induced large voids on the Moon, this indication further supports the hypothesis that impact cratering could increase near-surface porosity (Alejano and Alonso 2005;Pilkington and Grieve 1992).

Indications of future lunar exploration
Subsurface voids are ideal natural shelters for future lunar surface explorations since they provide ideal insulation from cosmic radiation, meteoroid impacts, and temperature fluctuations. Therefore, large lava tubes are high-priority targets for future lunar explorations (Kaku et al. 2017). Impact-induced subsurface voids may not be as good as lava tubes for future lunar exploration since such voids would be irregular in shape and difficult to access. However, this work demonstrates the advantage of ground penetrating radar in detecting meter-scale lava tubes on the Moon. Indeed, Carrer et al. (2018) proposed using a multifrequency radar sounder to map the distribution of lava tubes. Future explorations using ground penetrating radar should also consider the possibility of widespread impact-induced large voids in the subsurface.

Conclusion
In the high-frequency LPR data returned by the Yutu rover, we noticed a zone that features phase reversals of reflections, and the polarity of the upper echo pattern is the same as that of the incident pulse. This feature is located within the top 50 ns of the radargram, and it was not noticed before. We carried out numerical simulations of the propagation of electromagnetic waves in materials that have different relative permittivities and subsurface structures to investigate their identity. Based on the morphology of the upper hyperbolic echo pattern of the zone, the general shape of the roof of the low-permittivity zone is constrained as convex upward, but the south-north width of this zone cannot be constrained. Based on the distribution pattern of the normalized echo strengths along the lower hyperbolic echo pattern, we compared the simulated and observed radargrams to estimate the possible range of relative permittivity. The results suggest that a relative permittivity of 1 for the low-permittivity zone is most consistent with the observation. This conclusion is based on a large number of simplified models, and the existence of a subsurface cavity is supported by the fact that lower permittivity in typical lunar regolith is most likely caused by cavities. The ~ 3.1-m-high cavity in the continuous ejecta deposits of Ziwei crater was most likely formed due to stacking of large excavated basaltic boulders. This interpretation is also consistent with the fact that large blocks excavated by Ziwei crater are mainly deposited close to the crater rim. The results suggest that large voids that formed during impact cratering contribute to the high porosity of the lunar shallow crust. Lunar penetrating radar is a valuable technique to detect subsurface cavities, such as buried lava tubes.
Additional file 1. Text S1. FDTD simulations for the propagation of Chang'e-3 high-frequency LPR in regolith models that contain variousshaped subsurface cavities. In all the models simulated, untypical phase reversals occur in both the radargram and A-scan. Text S2. FDTD simulations for the propagation of Chang'e-3 high-frequency LPR in regolith models that contain subsurface cavities with different shapes and geometries. Text S3. Other structures of lunar regolith that contain subsurface cavities could form radargrams similar with the observation. Text S4.
When building models to deduce the possible relative permittivity of the low-permittivity zone, the height of the low-permittivity zone needs to be adjusted according to the assumed relative permittivity. The basic geometry of the low-permittivity zone is referred from that shown in Figure 3f. Figure S1. Regolith models contain subsurface cavities with different roof geometries. Figure S2. Modeling radargram for typical lunar regolith that has an internal cavity with a shape followed by triangle. Figure S3. Phase reversals formed due to a subsurface cavity (inverse triangle in shape, but with various sizes) in regolith. Figure S4. Phase reversals formed due to a subsurface cavity (rectangular in shape, but with various sizes) in regolith. Figure S5. Phase reversals formed due to a subsurface cavity (trapezoid in shape, but with various sizes) in regolith. Figure S6. Phase reversals formed due to a subsurface cavity (inversed trapezoid in shape, but with various sizes) in regolith. Figure S7. Phase reversals formed due to a subsurface cavity (pentagon in shape, but with various sizes) in regolith. Figure S8. Phase reversals formed due to a subsurface cavity (inversed pentagon in shape, but with various sizes) in regolith. Figure  S9. Phase reversals formed due to a subsurface cavity (ellipse in shape, but with various sizes) in regolith. Figure S10. Phase reversals formed due to a subsurface cavity (semi-ellipse in shape, but with various sizes) in regolith. Figure S11. Phase reversals formed due to a subsurface cavity (irregular polygon in shape, but with various sizes) in regolith. Figure S12.
Alternative regolith models that contain subsurface cavities that can form similar radargrams with the observation. Figure S13. Schematic diagram showing depth estimation for the low-permittivity zone according to the assumed dielectric constant.