The adopted gravity data in this paper are derived from the GRGM1200A spherical harmonic coefficient model from National Aeronautics and Space Administration (NASA)’s GRAIL satellite (Lemoine et al., 2014; Goossens et al., 2016), which is a corrected Bouguer gravity field model with the degree of 1200 and covers the entire Moon.

The Bouguer gravity anomaly Δ g φ , λ in the Mare Serenitatis area can be expressed as (Wang et al., 2009; Liang, 2010): Δ g φ , λ = G M R 2 ∑ l = 2 ∞ ∑ m = 0 l l + 1 R r l + 2 C ¯ l m cos m λ + S ¯ l m sin m λ P ¯ l m cos φ (1) where G is the universal gravitational constant, M is the lunar mass, R represents mean lunar radius, r is the distance from the calculation point to the mass center of the Moon, C ¯ l m and S ¯ l m are the spherical harmonic coefficients provided by gravity field model, l and m are the respective degree and order, λ and φ are the respective longitude and colatitude, and P ¯ l m cos φ is the regularized Legendre function. Δ g φ , λ reflects the sum of all subsurface anomalous materials. In order to extract the corresponding gravity signals of targeted materials at different depths, Δ g φ , λ needs to be further decomposed.

According to the principle of wavelet multi-scale analysis (Mallat, 1989), Δ g φ , λ can be decomposed into a low-frequency part and high-frequency parts of different orders, as shown in Eq. 2. Δ g φ , λ = A j φ , λ + ∑ j = 1 J D j φ , λ (2) where A j φ , λ is the j th order wavelet approximation, reflecting the low frequency signal; D j φ , λ is the j th order wavelet detail, reflecting the high frequency signal; and J is the maximum order of decomposition.

Gravity anomalies in various frequency bands can be regarded as the signals of anomalous materials at different depths. The average field source depth H j of D j φ , λ is further estimated by the radial power spectrum method (Spector and Grant, 1970) H j = Δ ln P j 4 π Δ k (3) where ln P j is the logarithm of the power spectrum for D j φ , λ , k is the wave number, and Δ is the variability.

According to the mean field source depth H j   j = 1,2,3 . . . , J , the lunar crust and upper mantle of the Mare Serenitatis is firstly layered. Subsequently, each layer is further gridded by Tesseroids. The function between the decomposed gravity anomaly D j φ , λ and the density anomaly Δ ρ j φ 0 , λ 0 of each Tesseroid at corresponding layer can be expressed as (Heck and Seitz, 2007): D j φ , λ = ∑ φ 0 ∑ λ 0 G Δ ρ j φ 0 , λ 0 Δ r j Δ φ Δ λ L 000 + 1 24 L 200 Δ r 2 + L 020 Δ φ 2 + L 002 Δ λ 2 (4) where φ 0 and λ 0 are the respective central latitude and longitude of the Tesseroid, Δ φ and Δ λ are the respective latitude and longitude intervals, and Δ r j is the thickness of the Tesseroid. L 000 , L 200 , L 020 and L 002 are the correlation coefficients as shown in Eqs 5–8 respectively. L 000 = r 0 2 R − r 0 cos ψ 0 cos φ 0 l 0 3 (5) L 200 = R cos φ 0 l 0 3 { 2 − 3 r 0 l 0 2 5 r 0 − 2 R + 3 r 0 cos Ψ 0 cos Ψ 0 + 15 r 0 3 l 0 4 sin 2 ψ 0 r 0 − R cos ψ 0 } (6) L 020 = r 0 l 0 3 cos φ 1 − 2 sin 2 φ 0 cos δ λ + r 0 2 l 0 5 { − R R 2 + r 0 2 cos φ 0 + r 0 sin φ [ − R r 0 sin φ cos φ 0 − cos φ sin φ 0 cos δ λ + sin φ 0 cos φ 0 2 R 2 + 4 r 0 2 − 3 R r 0 sin φ sin φ 0 ] + r 0 2 cos φ cos δ λ 1 − 2 sin 2 φ 0 × r 0 + R cos φ cos φ 0 cos δ λ + R r 0 2 cos φ sin φ 0 cos φ 0 cos δ λ × 3 sin φ cos φ 0 − 4 cos φ sin φ 0 cos δ λ + 5 R r 0 3 l 0 7 { − R R 2 + r 0 2 sin φ 0 + r 0 2 cos φ sin φ 0 cos φ 0 cos δ λ × r 0 + R cos φ cos φ 0 cos δ λ + r 0 sin φ 2 R 2 − r 0 2 − R r 0 cos ψ 0 + sin 2 φ 0 × R 2 + 2 r 0 2 − R r 0 sin φ sin φ 0 } × sin φ cos φ 0 − cos φ sin φ 0 cos δ λ (7) L 002 = r 0 l 0 3 cos φ cos 2 φ 0 × { cos δ λ − 3 R l 0 2 [ 2 r 0 cos φ cos φ 0 sin 2 δ λ + R − r 0 cos ψ 0 cos δ λ ] + 15 R 2 r 0 l 0 4 cos φ cos φ 0 R − r 0 cos ψ 0 sin 2 δ λ } (8) where r 0 = R − H j , δ λ = λ 0 − λ , l 0 = R 2 + r 0 2 − 2 R r 0 cos Ψ 0 , cos Ψ 0 = sin φ sin φ 0 + cos φ cos φ 0 cos δ λ . Simplify and write Eq. 4 into the matrix form, as shown in Eq. 9 D j = B Δ ρ j (9) where B is the kernel function matrix. At last, Δ ρ j φ 0 , λ 0 is solved by using the Tikhonov regularization method (Tikhonov and Arsenin, 1977): Δ ρ j = B T B + α I − 1 B T D j (10) where I is the unit matrix and α is the regularization factor, which is determined by L-curve method (Hansen and O’Leary, 1993).

It is a key factor for obtaining refined crust-mantle interface depth to accurately extract the gravity anomaly signal caused by the crust-mantle interface relief, which needs to remove the effect of sediments, consolidated crystalline lunar crust and lower lunar mantle (Wan et al., 2019). The common methods for these corrections depend on existing lunar crustal and mantle models. However, at present the lunar crustal and mantle models are not accurate, which may cause large errors during the corrections. Hence, in this paper the gravity anomaly signal g B A caused by the crust-mantle interface relief of the Mare Serenitatis is extracted by combining D j φ , λ , its corresponding field source depth H j , and the a priori information on the available crust-mantle interface depth. We compared the estimated H j with the average depth of crust-mantle interface provided by a priori information to determine which order of D j φ , λ was the gravity anomaly signal caused by the crust-mantle interface undulation.

Subsequently, the crust-mantle interface depth h φ , λ of the study area is inverted from g B A using the effective power topography method (Wieczorek and Phillips, 1998): h φ , λ = w l g B A M 2 l + 1 4 π Δ δ Z 2 ( R / Z ) l − Z ∑ n = 2 l + 3 h φ , λ n Z n n ! ∏ s = 1 n l + 4 − s l + 3 ) (11) where Δ δ is the density contrast at the crust-mantle interface, l is the degree of spherical harmonic coefficients, Z is the distance from the mean depth of crust-mantle interface to mass center of the Moon, and w l is the downward topological filter, as shown in Eq. 12 w l = { 1 + μ [ M 2 l + 1 4 π Δ δ Z 2 R / Z l ] 2 } − 1 (12) in which μ is the Lagrange multiplier. The determination of μ is a subjective process, where the larger μ is, the more the high-frequency signal will be filtered. According to Wieczorek and Phillips (1998), here we choose μ such that w 30 = 0.5 . The final crust-mantle interface depth h φ , λ will be determined by the iteration of Eqs 11, 12.

Because the vertical density variation of the layered density inversion is not continuous, the profile density inversion is further conducted in order to recover a more refined 3-D structure of the lunar crust and upper mantle in the Mare Serenitatis region. The compact gravity inversion method proposed by Last and Kubik (1983) is adopted to construct the profile density model in this paper. The profile gravity anomaly g i can be written as: g i = ∑ j = 1 M a i j v j + e i , i = 1,2,3 . . . , N (13) where v j is the profile density anomaly of j t h block, e i is the noise associated with i th data point, having an initial value of 0.1 g i . a i j is the coefficient, as shown in Eq. 14 a i j = 2 G [ ( x i − x j + d / 2 ) log r 2 r 3 / r 1 r 4 + d × log r 4 / r 3 − ( z j + h / 2 ) θ 4 − θ 2 + ( z j − h / 2 ) θ 3 − θ 1 ] (14)

Among them, the parameters are: r 1 = ( z j − h / 2 ) 2 + ( x i − x j + d / 2 ) 2 (15) r 2 = ( z j + h / 2 ) 2 + ( x i − x j + d / 2 ) 2 (16) r 3 = ( z j − h / 2 ) 2 + ( x i − x j − d / 2 ) 2 (17) r 4 = ( z j + h / 2 ) 2 + ( x i − x j − d / 2 ) 2 (18) θ 1 = arctan ( x i − x j + d / 2 ) / ( z j − h / 2 ) (19) θ 2 = arctan ( x i − x j + d / 2 ) / ( z j + h / 2 ) (20) θ 3 = arctan ( x i − x j − d / 2 ) / ( z j − h / 2 ) (21) θ 4 = arctan ( x i − x j − d / 2 ) / ( z j + h / 2 ) (22) where x i is the coordinate of g i , x j , z j is the coordinate of j t h block, d and h are the length and width of each block. According to the least square method, the profile density anomaly vector V is solved by: V = W v − 1 A T A W v − 1 A T + W e − 1 − 1 g (23) where A is the coefficient matrix of a i j , W v is the density weighting matrix, and W e is the noise weighting matrix. These two weight matrices are shown in Eqs 24, 25, respectively. W v − 1 = V 2 + ε (24) W e − 1 = l 0 2 A W v − 1 A T (25) where ε is taken as 10⁻⁸ and l 0 is the ratio of estimated noise to signal (initial value is 0.1).

Figure 2 shows the Bouguer gravity anomaly at the elevation of 0 km with the spatial resolution of 0.25 ∘ × 0.25 ∘ in the Mare Serenitatis area using Eq. 1, where GM product is 4902.801374 × 10⁹ km³/s² and R is 1738 km (Liang, 2010).

In Figure 2, the Bouguer gravity anomalies range from −186 mGal to 336 mGal. Compared to the surrounding area, the gravity anomaly of the Mare Serenitatis Basin (9°E–29°E and 17°N–36°N) is positive and high, which is pear-shaped. The gravity anomaly is the highest in the center of the Mare Serenitatis Basin and gradually decreases in all directions, showing a clear Mascon character. In addition, there is a negative correlation between the Bouguer gravity anomaly and the topography (Figure 1). The gravity anomaly highs in the center of the Mare Serenitatis correspond to lowlands with the elevation of approximately −5 km. However, the gravity anomaly lows are consistent with the mountains, such as the MC, the MHL and the MH. Moreover, the ridges are primarily located at the gravity high-low transitional zones, such as the DG, the DS and the SR.

The Bouguer gravity anomaly is the comprehensive reflection of all subsurface anomalous materials. It cannot directly reflect the distribution of the material at different depths. Thus, the Bouguer gravity anomaly is further separated by wavelet multi-scale method. According to the previous research on the selection of the optimal wavelet basis (Xu et al., 2017), this paper adopts the “coif3” wavelet basis and performs a 2D wavelet decomposition for the Bouguer gravity anomaly of the Mare Serenitatis using Eq. 2. The decomposed results are shown in Figure 3, and their mean field source depths estimated by Eq. 3 are shown in Figure 4 and Table 1. In the subplots of Figure 3, D1–D8 represent the order one to eight wavelet details, respectively. The subplots in Fig. 4 are the corresponding radial logarithmic power spectrum of wavelet detail D1–D8, respectively.

According to Figure 3, the gravity anomalies in D1 are very weak and may be interference signals. The gravity anomaly values of D2 vary slightly from −0.3 mGal to 0.3 mGal, whose average source depth is approximately 4.7 km corresponding to the distribution of shallow sediments. Gravity anomaly values of D3 with the average source depth of 13.4 km range from −3.0 mGal to 3.0 mGal. Small positive and negative alternating gravity anomaly circles are beginning to appear, obviously in the south of DS. The average source depth of D4 is approximately 23.1 km and the gravity anomaly values range from −12.0 mGal to 12.0 mGal. The gravity anomaly circles are further expanded, which are mainly distributed along the outer ring of the Mare Serenitatis, such as DO, DS and DAD. It indicates that tectonic structure begins to become complex. In addition, there are apparent gravity anomaly signals during 15°E–25°E and 35°N–42°N, which is highly consistent with the location of the center of the “D ridge system” mentioned by Maxwell et al. (1975). This could be the signals of anomalous material left over from the meteorite event that formed the “D ridge system”. Based on the lunar crust thickness model of Wieczorek et al. (2013), the mean source depth of D4 should be roughly at the interface between the lunar crust and mantle in the Mare Serenitatis region. The gravity anomaly values of D5 with the average source depth of 33.4 km range from −40 mGal to 40 mGal. The gravity anomaly circles are larger. The mean source depth of D6, with gravity anomaly values ranging from −30 mGal to 30 mGal, is approximately 51.7 km. During the 20°E–27°E and 23°N–30°N corresponding to the inner ring region between DA and SR (Maxwell et al., 1975), there is the first large anomaly in the Mare Serenitatis Basin. The center of this large anomaly is at approximately 23°E and 25°N. The D7 shows the distribution of gravity anomalies with the average source depth of 66 km. The values range from −50 mGal to 100 mGal, which is the largest variation ranges for D1–D8. There is another obvious gravity anomaly high, the second large anomaly, whose center is located at 15°E and 25°N. Compared with the first large anomaly mentioned in D6, the second large anomaly in this area is stronger, deeper and more widespread. Therefore, this paper agrees with Watters and Konopliv (2001) that the Mare Serenitatis consists of two overlapping basins. The western basin corresponds to a larger impact forming the outer ring and the eastern basin is in agreement with the inner ring structure of the Mare Serenitatis Basin. The mean source depth of D8 is approximately 86.4 km and the gravity anomaly values range from −10 mGal to 20 mGal. It can be seen that there are no longer any obvious gravity anomaly circles, and the anomaly signals gradually weaken and disappear. It indicates that the tectonic structure at this depth tends to be stable.

According to above analysis, D1–D3 primarily reflects the gravity anomaly distribution of the lunar crust in the Mare Serenitatis area. With depth increasing, the gravity anomaly circles and magnitude gradually become larger. D4—D5 show obvious positive-negative alternating gravity anomalies, indicating the complex tectonic structure. D6—D7 present the morphological characteristics of the two large gravity anomaly highs in the southeast and southwest, which are presumed to be closely related to the inner and outer ring structure of the Mare Serenitatis. D8 indicates that there are no more obvious anomalies at the depth of 86.4 km.

Firstly, the lunar crust and upper mantle of the Mare Serenitatis is layered according to the corresponding mean field source depths of the decomposed Bouguer gravity anomalies, see Table 1. Then, each layer is gridded and modelled using Tesseroids and the size of each Tesseroid is 0 . 5 ∘ × 0 . 5 ∘ ( Δ φ = 0 . 5 0 ; Δ λ = 0 . 5 0 in Eq. 4. Lastly, the density distribution at the corresponding depth was calculated using Eq. 10 and the results are shown in Figure 5.

The results of D1 and D2 are not presented in Figure 5 because of their small density fluctuations. D3 presents the density distribution at the depth of 13.4 km underneath the Mare Serenitatis. The density values range from 2.88 g/cm3 to 2.92 g/cm3. The circles of density highs and lows begin to appear mainly around the Mare Serenitatis basin. D4 corresponds to the density distribution at the depth of 23.1 km, which is approximately at the interface between the lunar crust and mantle in the Mare Serenitatis Basin. The values of density distribution range from 2.90 g/cm³ to 3.10 g/cm³. Compared with D3, the circles of density highs and lows become larger. It is noteworthy that in the 15°E–25°E and 35°N–42°N (the center of the aforementioned “D ridge system”), there are several obvious density high circles, which are presumed to be meteorite impact residue. Besides, there is a small high-density circle at (30°E, 27°N), which may be caused by the meteorite impact corresponding to LS in Figure 1. It indicates that meteorite impacts may induce density heterogeneity, which has been also discussed by Kierfer et al. (2012). The density distribution of D5 reflects the materials at the depth of 33.4 km. The high-density areas are discrete, implying the complex tectonic structure. The D6 corresponds to the density distribution at the depth of 51.7 km beneath the Mare Serenitatis. There is an obvious high density at (23°E, 25°N), which may be the result of a meteorite impact event forming the inner ring structure. D7 corresponds to the density distribution at the depth of 66 km. The whole Mare S. basin is nearly covered by a high-density material, whose center is at (15°E, 25°N). This is the second high density circles within the Mare Serenitatis, and it is presumed to be closely related to the outer ring structure of the Mare Serenitatis. The density distribution of D8 with the depth of 86.4 km has become smooth, indicating that the tectonic structure is beginning to stabilize.

In summary, there are two high-density bodies with different locations, depths and sizes in the Mare Serenitatis Basin. The center and mean depth of the southeastern high-density body is (23°E, 25°N) and 51.7 km, respectively. The southwestern one is centered at (15°E, 25°N) and has the mean depth of 66 km. These two high-density materials are presumed to be closely related to the inner and outer ring structure of the Mare Serenitatis (Maxwell et al., 1975; Watters and Konopliv, 2001).

According to the a priori information (Hikida and Wieczorek, 2007; Wieczorek et al., 2013), it is known that the average crust-mantle interface depth in the Mare Serenitatis region is approximately 25 km, which is close to the average field source depth of D4 in Figure 3. Thus, the 4th-order wavelet approximation A4 (see Figure 6) of the Bouguer gravity anomaly is regarded as the signal of crust-mantle interface relief, and then the crust-mantle interface depth of the Mare S. region is determined by the iterative inversion of Eq. 11, where the average depth and density contrast of crust-mantle interface are set as 25 km and 0.56 g/cm³, respectively (Hikida and Wieczorek, 2007; Wieczorek et al., 2013). Figure 7 shows the determined crust-mantle interface depth underneath the Mare Serenitatis.

On the whole, the crust-mantle interface depth gradually shallows from the outside to the inside of the Mare Serenitatis, see Figure 7. The depth in MC, MHL and MH is approximately 35 km, and that in LS, PO and LM is between 20 and 26 km. The depths of the marginal area (green part in Figure 7) of the Mare Serenitatis are approximately 22 km. Further, the light blue area has the average depth of approximately 15 km, and the dark blue area has the average depth of about 10 km. The shallowest depth, located at the south of the Mare Serenitatis Basin (east and west of BE), is only approximately 4 km, whose positions are correspondence with the two high-density materials aforementioned.

According to above analysis, the crust-mantle interface of the Mare Serenitatis is obviously uplifted, which is in good agreement with the result of Wieczorek et al. (2013). The reason may be the post-impact rebound after a large meteorite impact event, which generated thermal energy perturbation leading to lateral temperature differences in the lunar mantle, followed by the rapid and uniform uplift of the mantle (Wieczorek and Phillips, 1999; Chen et al., 2009). In addition, there is a new finding that the uplift of crust-mantle interface falls back in the center of the Mare Serenitatis, see the white dashed lines in Figure 7. The deepest depth surrounded by the white dashed lines is approximately 17 km, which is deeper than the depth of the dark blue area. It suggests that the uplift of the crust-mantle interface is irregular. Combined with the results of layered density inversion (Figure 5), it further demonstrates that the Mascon in the Mare Serenitatis area may be formed by the combination effect of the high-density materials and crust-mantle interface uplift.

Further, we extracted the lunar seismic inversion results of P1and P2 (Figure 7) from Chenet et al. (2006) and compared them with the crust-mantle interface depth in this study, as shown in Table 2. The difference of crust-mantle interface depths between Chenet et al. (2006) and this study is 0.1 km at P1 and 3.75 km at P2, respectively. The inverted crust-mantle interface depths in this study are in agreement with those provided by Chenet et al. (2006), which verifies the correctness of our results.

Firstly, a simple synthetic test is designed to verify the correctness of the compact gravity inversion method proposed by Last and Kubik (1983), which will be applied to invert the profile density of the lunar crust and upper mantle in the Mare Serenitatis. In the synthetic test, we set up several modules with different shapes and depths, as shown in the lower subplot of Figure 8, where each module (deep red areas) has a relative density value of 0.3 g/cm³ and other areas (white regions) have a relative density value of 0 g/cm³. Based on this model, we calculate the gravity anomaly (Gmodel), as shown at the upper subplot (red line) of Figure 8. Then, this calculated gravity anomaly is considered as input (Gobs in the upper subplot of Figure 9), and the compact gravity inversion method is employed to inverted the profile density. The inverted results using Eq. 23 are shown in the lower subplot of Figure 9. The inverted density distribution in Figure 9 is in good agreement with the original model of Figure 8, which demonstrates that the primary density anomalous bodies can be recovered effectively by the compact gravity inversion method. In addition, the RMS of the difference between Gmodel and Gobs is 0.0004 mGal. Thus, it indicates that the compact gravity inversion method is correct according to the synthetic results. It is also worth noting that there are some differences between the lower subplot of Figures 8, 9, such as the light blue areas in Figure 9. The reason may be that the shape of synthetic modules in Figure 8 is regular and the values of density setting are not continued, which may cause signal leakage and distortion during the inversion.

Subsequently, the compact gravity inversion method is used to invert the density distribution of two intersecting profiles, E-F and Q-O (as shown in Figure 5), in order to further reveal the detailed morphological characteristics of two high-density materials in the Mare Serenitatis. E-F and Q-O just go through the centers of two high-density materials, whose profile density anomalies (anomalies relative to the average surrounding density) are inverted as shown in Figures 10, 11. During the profile density inversion, the length and width of each cell module are set as 20 km and 5 km respectively, and the total depth of the inversion is 100 km (the density anomalies at deeper depths are no longer apparent). Two subplots are included in Figures 10, 11. The upper subplot shows the comparison between Gmodel and Gobs, and the RMS of the difference between Gmodel and Gobs for E-F and Q-O is 3.886 mGal and 2.4706 mGal, respectively. It indicates that the inverted profile density anomalies are reliable. The lower subplot shows the profile density anomaly distribution.

In Figure 10, there is an obvious large-scale high-density anomalous material, which is located at the distance of 430 km–530 km from Point E with the depth of 30 km–60 km. It is the first high-density anomalous body mentioned in the Figure 5. In Figure 11, there are two obvious larger-scale high-density anomalous materials. The left one is located at 110 km–260 km from Point Q with the depth of 50 km–80 km, which is the second high-density anomalous body mentioned in the Figure 5. The right one is consistent with the first high-density anomalous body presented in Figure 10. In addition, the falling phenomenon of the lunar crust-mantle interface uplift in the center Mare Serenitatis can be observed clearly both in the Figures 10, 11. Above the two high-density anomalous materials, there are apparent low-density materials, which may be formed by the filled magma after meteorite impacts.

Lastly, we compare the inverted profile density (Figures 10, 11) with the results of the same profiles from the inverted layered density (Figure 5) as shown in Figure 12 and Figure 13. Comparing the results in Figures 10, 12, the locations of the first high-density anomalous body are similar, which are both at the distance of approximately 500 km from Point E with the depth of approximately 50 km. Above the first high-density anomalous body, there are both low-density materials in Figures 10, 12. Comparing the results in Figures 11, 13, the locations of the second high-density anomalous body are also consistent, which are both at approximately 200 km from Point Q with the depth of approximately 70 km. Above the second high-density anomalous body, there are also both low-density materials in Figures 11, 13. Thus, on the whole, the inverted profile density (Figures 10, 11) and the inverted layered density (Figures 12, 13) can match well at the large-scale features. However, there are also some apparent differences in details between the inverted profile density (Figures 10, 11) and the inverted layered density (Figures 12, 13). For example, there is an obvious low-density anomaly at the distance of 0–110 km from Point E with the depth of 25–40 km in Figure 12, while it is not existed in Figure 10. The reason may be that the adopted methods with the assumptions of the profile density inversion and layered density inversion are all different. According to above analysis, the inverted profile density and the inverted layered density can complement each other, which may provide more information for understanding the lunar crust and upper mantle structure in the Mare Serenitatis.