Unveiling magmatic structures and connectivity beneath the lunar Oceanus Procellarum region from GRAIL gravity data

Domain: earth_science
Published: 2024-11-27
Editorial outcome: major_revision

Abstract
The Oceanus Procellarum region, characterized by its vast basaltic plains and pronounced volcanic activity, serves as a focal point for understanding the volcanic history of the Moon. Here we present density models of the magmatic structures beneath Oceanus Procellarum, derived from Gravity Recovery and Interior Laboratory (GRAIL) mission data. The models uncover pronounced linear magmatic structures along the Procellarum’s western border and significant intrusions within the northern and southern Marius Hills. Crucially, they reveal three narrow near-horizontal sheeted magmatic conduits, 80-150 km long, extending from near-surface to 6-7 km depth. These magmatic conduits connect the Marius Hills’ northern and southern intrusions and bridge them with the Procellarum western border structures, suggesting that they likely served as central pathways facilitating magma transport across various volcanic systems. These discoveries reveal widespread magmatic connectivity beneath Oceanus Procellarum and underscore the critical role of lateral magma transport processes in shaping the Moon’s volcanic evolution.

Introduction
Oceanus Procellarum, the largest mare on the Moon’s near side, is characterized by extensive basaltic plains and significant volcanic activity1. Uncovering its subsurface structures is pivotal for understanding lunar volcanism. Gravity data, which have played a crucial role in studying magma plumbing systems on rocky planets and satellites2–6, provide the ideal tool for such an investigation of Oceanus Procellarum. Forward modeling and inversion of gravity anomaly data enable the creation of density models of various magmatic structures that exhibit a density contrast with the surrounding rocks, such as shallow intrusions7–9, magmatic-related deposits10, and plutons11,12. Data from NASA’s Gravity Recovery and Interior Laboratory (GRAIL) mission have generated a high-resolution global map of the Moon’s gravity field. By modeling the gravity field to a spherical harmonic degree and order of 1200, a spatial resolution of 4.5 km by 4.5 km at the Moon’s equator was achieved13,14. Combined with topographic data from the Lunar Orbiter Laser Altimeter (LOLA), these data aid the detection of small-scale features such as dykes7, ancient igneous intrusions3, buried impact craters15, and empty lava tubes16,17.

The volcanic legacy of the Oceanus Procellarum is exemplified by the Marius Hills volcanic complex (MHVC), a 35,000 km2 plateau, which exhibits the densest concentration of extrusive volcanic features on the Moon in the form of volcanic domes/shields and cones18,19. This region has consistently emerged as a prominent site for lunar exploration, primarily due to its distinctive geological characteristics and the valuable insights it offers into lunar volcanism20–23. Utilizing gravity data from the Lunar Prospector spacecraft24 and the GRAIL mission25,26, two quasi-circular, positive gravity anomalies within the MHVC (north and south) were identified and modeled as two massive sill/dyke groups that are responsible for the volcanic activity in this region. Investigation into the surface distribution of volcanic vents confirms the presence of the zones of magma accumulation27. However, magmatic structures beneath the MHVC and surrounding areas and processes by which magma is transported to the shallow crustal depths remain enigmatic to date28,29.

GRAIL gravity data revealed a pattern of narrow linear anomalies that border the Procellarum KREEP Terrane (PKT)5. The source of anomalies is interpreted to be solidified lava within rifts and the feeder dykes beneath them, which likely served as the magma plumbing system for much of the nearside mare volcanism. The MHVC is positioned roughly 200 km east of the PKT’s western border structures. Previous research has primarily examined these volcanic systems independently5,24–27,30, overlooking the critical magmatic structures between them and their collective impact on volcanism in the Oceanus Procellarum. This oversight underscores a significant gap in our understanding of the magmatic architecture beneath the MHVC and the PKT’s western border, highlighting the importance of considering these systems as interconnected entities to comprehend the extensive volcanic activity fully in the Oceanus Procellarum region.

In this work, to better understand the magmatic architecture beneath the Oceanus Procellarum, we apply a three-dimensional (3D) inversion technique to the GRAIL gravity data to construct detailed 3D density models. The magmatic structures, characterized by higher density than the host rocks, can be recovered through the inversion of the Bouguer gravity anomaly, which plays a key role in understanding the magmatic sources and processes2,7,9,31. This study encompasses the southern and northern MHVC, the western border structures of the PKT, and their neighboring areas, enabling a comprehensive understanding of their interconnections, and facilitating their integration into the wider PKT region. The detailed modeling of these magmatic structures provides valuable information regarding the formation of volcanic landforms, magma transport, and energy exchange across various lunar systems.

Results
Bouguer gravity data and topography
Figure 1a shows the topography of the lunar near side derived from the MoonTopo2600p model based on data collected by the Lunar Reconnaissance Orbiter mission32. GRAIL gravity data extracted from the GRGM1200B model13 were used for this study. The Bouguer gravity anomaly (Fig. 1b) used for imaging the subsurface structures was derived from SHTools33 with a reference radius of 1738 km. To remove the gravitational contribution of the surface topography, a crust density of 2800 kg/m3 was used for Bouguer correction34. To eliminate long-wavelength variations in the crustal structure and short-wavelength noise, we filtered the Bouguer gravity data to include data between degree and order 7 and degree and order 660, yielding a half-wavelength surface resolution of ~8 km. A cosine taper was implemented between degrees 550 and 660 to reduce high-frequency ringing35.

Study area A (Fig. 1) covers the MHVC and partial western border structures of the PKT5. The Bouguer gravity anomaly (Figs. 1b and 2b) reveals two quasi-circular mass anomalies, marked as Marius-NH (northern half) and Marius-SH (southern half) based on their locations relative to the Marius crater26. The two anomalies exhibit maximum amplitudes of 165 and 169 mGal, respectively. Marius Hills plateau, corresponding to the northern anomaly, contains most domes and cones situated in a localized highland. In contrast, the southern anomaly covers fewer volcanic domes without elevated topography. A narrow linear positive gravity anomaly, ranging from ~40 to ~70 mGal, connects these two quasi-circular anomalies, forming a “bridge” that has been identified in an earlier study24. Furthermore, a significant high-gravity anomaly belt, spanning ~800 km with a NNW–SSE orientation and peaking at 183 mGal, is identified. This belt has been interpreted to be part of the magmatic border structures of the PKT region5. The linear gravity anomaly extends NNW and diverges into three thinner branches at the northwestern border of the PKT (Figs. 1b and 3b). Study area B (Figs. 1 and 3a, b) includes the western and middle branches of the northwestern border, each marked by distinct linear gravity anomalies with a general NNW–SSE orientation, spanning over 500 km and displaying values from ~30 to 138 mGal. The western edges of the belt gravity anomalies (Figs. 2b and 3b) are generally coincident with steep topographic gradients leading up to the adjacent highlands (Figs. 2a and 3a), while their eastern sides show less topographical correlation. Prominent wrinkle ridges (Supplementary Fig. 1) run parallel to the belt gravity anomaly in study area A and the belt gravity anomalies of the middle and eastern branches in study area B, exhibiting a NNW–SSE orientation. Also, a distinct set of wrinkle ridges (Supplementary Fig. 1) traversing the Marius Hills region predominantly adheres to this NNW–SSE trend36.

3D inversion density models from Bouguer gravity data
The 3D inversion technique applied here is based on the probabilistic method37 to generate 3D density models by inverting gravity data. The method was assessed and proven successful in identifying the density structures of both the upper and lower crust of the Earth9,38–44. A significant advantage of this inversion approach is its capacity to integrate a priori information regarding physical properties, spatial extent, depth, and orientation of the source bodies via the covariance matrix to construct an accurate and reliable solution. To address the non-uniqueness inherent in the inversion process, we performed a sensitivity analysis on several critical parameters that could affect the inversion results, including the covariance matrix parameter, the depth weighting function, and the depth extent of the inversion domain (see “Methods” section for more details).

Figures 2c–e and 3c–e display the 3D density models for the two study areas, which represent optimal models obtained from inversions using the best-fitting parameters (see “Methods” section for more details). The 3D density models recovered from the Bouguer gravity data with a half-wavelength surface resolution of ~8 km facilitate the investigation of complex subsurface magmatic structures in detail. Given the geological evidence for potential intrusive activity24,25,27, the 3D continuous density inversions are volume-rendered using discrete cutoffs to delineate the extent of possible intrusive bodies. The red surfaces indicate the density contrast iso-surface of 160 kg/m3, which outlines high-density bodies that we interpret as magmatic bodies. The yellow surfaces depict a higher density contrast iso-surface of 350 kg/m3 to delineate even denser magmatic bodies. The top depth corresponds to either the shallowest limit of the magmatic bodies or the deepest extent of impact brecciation, where the density contrast fades. The bottom depth denotes the base of the magmatic bodies or where the density contrast with surrounding rocks significantly diminishes.

Figures 4 and 5 show horizontal and vertical cross-sections through the inverted density models for the two study areas. The blue dashed lines indicate the 160 kg/m3 density contrast contour lines. The assumed density contrast iso-surface of 160 kg/m3 for the magmatic body does not consider potential density changes amongst the magmatic rocks, and could over- or underestimate the size and depth of the magmatic bodies and thus presents a limitation and potential source of error in interpretation. When using the iso-surfaces of 150 and 170 kg/m3 as the interpreted magmatic bodies, the total volume increases by 5% and decreases by 6%, respectively, compared to the 160 kg/m3 iso-surfaces.

Figures 2c, 3c, 4, and 5 exhibit distinct belt structures oriented NNW–SSE with density contrasts ranging from 160 to ~700 kg/m3. These belt structures correspond to the belt gravity anomalies depicted in Figs. 2b and 3b and match the PKT border magmatic intrusions as identified in a previous study5. We observe that the density contrast of the intrusions from the inversion gradually decreases from the center to the edges. This is mainly due to the smoothing nature of the inversion method used44; the actual density contrast variation within the intrusions may not be as pronounced. The density model for area A (Figs. 2c and 4) shows that the border intrusion spans from near the surface down to depths of 10–13 km with a length of ~800 km, a thickness of 8–11 km, and a width of 130–270 km. In study area B, the western branch intrusion is ~600 km in length, ~140 km in width, and ~8–10 km in thickness, and the middle branch intrusion is about 550 km long, 80 km wide, with a thickness ranging from 5 to 9 km. The western PKT border intrusions become thinner and narrower, and the depth of the top surface increases toward the north.

The inversion density model indicates that the Marius-NH intrusion (Figs. 2d and 4f, h) and Marius-SH intrusion (Figs. 2e and 4e, h), reported in the literature24,25, extend down to depths of 11 and 12 km, with thicknesses of 10.2 and 9.7 km, respectively. The Marius-NH body spans ~100 km in diameter, while the Marius-SH body spans ~130 km. The density contrasts of these two intrusions lie between 160 and 720 kg/m3, with the maximum value in the upper center of the intrusions. Both intrusions have a bowl-like shape. Moreover, a significant magmatic body (located at 55.7°W, 12.9°N; see Figs. 2c, d and 4a, b, f) characterized by a quasi-circular configuration is located adjacent to the west of the Marius-NH intrusion. This body has a diameter of ~90 km, a thickness of ~6 km, and is positioned beneath the more subdued western inclines of the elevated Marius Hills plateau. A magmatic body (located at 55.3°W, 6.3°N; see Figs. 2c, e and 4a, b, e) of similar dimensions, with a diameter of ~80 km and a thickness of ~5 km, is identified. The intrusion is located directly beneath the Reiner crater at a shallow depth (3.5 km; Fig. 4a) and tilts slightly southwestward at a deeper depth (5.5 km; Fig. 4b). Both of these two magmatic bodies were neither discovered nor mentioned in previous studies15,24–26.

It is noteworthy that the inversion density model for study area A (Figs. 2c and 4) reveals a complex network of connectivity among the various volcanic systems under investigation. The Marius-NH and Marius-SH intrusions are interconnected with a narrow corridor of dense rocks extending ~80 km in length, 50 km in width, and 5 km in thickness (Figs. 2c and 4a, b, h). The narrow corridor is sill-like in structure and nearly horizontal. Furthermore, another two sill-like narrow corridors connect the two intrusions with the western border structures. The southern magmatic conduit, which traverses the intrusion beneath the Reiner crater, is ~130 km long, 50 km wide, and 5 km thick (Figs. 2c and 4a, b, e). The northern magmatic conduit, crossing the quasi-circular intrusion adjacent to the west side of the Marius-NH intrusion, measures ~110 km in length, 55 km in width, and 6 km in thickness (Figs. 2c and 4a, b, f).

A lower density of 2560 kg/m3 and a higher density of 3000 kg/m3 were also used for Bouguer corrections to derive the Bouguer anomalies, which were then inverted to construct density models, as detailed in Supplementary Fig. 2. The impact of conducting the Bouguer correction using different density values on the resultant Bouguer gravity anomaly is especially pronounced in regions with significant terrain variations, such as the Reiner crater (Supplementary Fig. 3a) and Marius-NH (Supplementary Fig. 4a). The higher density of 3000 kg/m3 yields the best-fit Bouguer correction for the craters (Supplementary Fig. 3a), which is similar to the mean density estimated for the mare45 (~3090 kg/m3) and the bulk density estimated for the Marius Hills surface basalts (3150 kg/m3) based on observations from the Lunar Prospector Gamma Ray Spectrometer46,47. However, it also significantly enhances the western belt gravity anomalies (Supplementary Fig. 3a) due to the relatively low topography in that area. Ultimately, we opted for a density of 2800 kg/m3 as a more representative choice for the entire region. The comparison of inversion models (Supplementary Fig. 2) shows that as the Bouguer correction density is increased, there is a slight increase in the volume of the western PKT border intrusion in study area A and the western branch in study area B. Conversely, the volumes of the middle branch intrusion in study area B, the Marius-NH intrusion, and the sill-like narrow corridors identified above decrease slightly. However, these effects are limited, and the main structural characteristics of the high-density bodies in this region exhibit highly consistent geometries.

To minimize the potential influence of crustal thinning on the inversion models, which is characterized by relatively long-wavelength positive gravity signals, we filtered signals with wavelengths <11 degrees (Supplementary Fig. 5a, c), and a subsequent 3D inversion was performed. The inversion results (Supplementary Fig. 5b, d) indicate a modest reduction in the lateral spread of the PKT border structures and a 2–3 km shallowing of their bases. These two inversion models may more accurately reflect the border structures by removing the effect of crustal thinning. Nevertheless, we observed reduction in both the volume and maximum density contrast of the intrusions beneath Marius-SH and Marius-NH as well as the sill-like narrow corridors (Supplementary Figs. 3e, 4e, and 5b), where mantle uplift is deemed improbable24,25,48,49. Moreover, the middle branch structures of the northwestern border become discontinuous (Supplementary Fig. 5d). Thus, it seems that filtering out gravity signals lower than degree and order 11 has simultaneously affected the long-wavelength signals associated with shallow sources.

Beyond study areas A and B, we also examined the Flamsteed region (labeled as C in Supplementary Fig. 1) and the Aristarchus plateau region (labeled as D in Supplementary Fig. 1) within the Oceanus Procellarum. These features were also analyzed, though the interpretation of magmatic structures remains uncertain. However, neither the southern Aristarchus plateau nor the Flamsteed region appears to connect with the western PKT border intrusions. As a consequence, our discussion is focused on the inversion models (Figs. 2–5) in study areas A and B, which are derived from the Bouguer gravity data spanning from degree and order 7 to degree and order 660, following Bouguer correction with a crust density of 2800 kg/m3.

Methods
The 3D inversion technique used in this research leverages a probabilistic method to generate 3D density models from gravity data37, which has been shown to be effective in identifying crustal structures on Earth9,38,39,94. The objective function employed for the inversion of gravity data is as follows:

\Phi = (G m - d_{obs})^T C_D^{-1} (G m - d_{obs}) + (m - m_{apr})^T C_M^{-1} (m - m_{apr}) (1)

Here, m represents the model parameter vector, values of the densities of the cells, d_{obs} is the gravity data vector, and G is the forward modeling matrix where Gm = d. The first term of Eq. (1) signifies the data misfit function, weighted by the data error covariance matrix, C_D. The second term corresponds to the stabilizing function, which accounts for the deviation of model parameters m from a reference model m_{apr}, weighted by the model covariance matrix, C_M.

The a priori information, including physical properties, spatial extent, depth, and orientation of the source bodies, can be incorporated via C_M to constrain the inversion and achieve an accurate and reliable solution.

To minimize the objective function, the solution is derived below37

m_e = m_{apr} + C_M G^T (G C_M G^T + C_D)^{-1} (d_{obs} - G m_{apr}) (2)

To counteract the natural decay of the forward modeling matrix with increasing depth, a depth weight function with the following form was used:

w(z) = 1 / (z_{max} - z)^\beta (3)

where z_{max} is the maximum depth of the inversion domain and z is the depth of each cell. \beta is the exponent of the depth weight function and \beta = 1 was used for this study.

The solution with the depth weighting function is:

m_e = m_{apr} + W^{-1} C_M W^{-1} G^T (G W^{-1} C_M W^{-1} G^T + C_D)^{-1} (d_{obs} - G m_{apr}) (4)

The parameters of the model covariance matrix C_M for areas A and B are listed in Table 1. We utilized a nested covariance model, as detailed in the second and third rows of Table 1, since the Bouguer gravity anomaly (Figs. 2b and 3b) within this region displays characteristics that span multiple scales. Specifically, the border structures of the PKT present larger linear formations trending NNW, while anomalies like the Marius-SH and Marius-NH appear to be relatively smaller quasi-circular in nature. Data error covariance matrix C_D is a diagonal matrix where each diagonal element is set at 2 mGal, representing an approximate estimate of the noise level in the Bouguer gravity data. A zero-reference model (m_{apr} = 0) was used for the inversions.

Before conducting the inversion, the Bouguer gravity data were transformed from spherical to Cartesian coordinates. The inversion domain was divided into numerous prism cells with dimensions of 1.5 km \times 1.5 km \times 0.25 km in the x, y, and depth directions, respectively. The depth extent of the inversion domain was set at 20 km. To mitigate edge effects during the inversion, the spatial dimensions of the inversion domain were expanded by 30 km in all lateral directions beyond the data coverage. This parameterization results in a model with 534 \times 526 \times 80 cells in study area A, 433 \times 467 \times 80 cells in study area B. The top of the inversion domain was constrained by the topographic data. Therefore, the cells above the terrain were removed and not inverted. By constraining the top of the terrain, it is possible to mitigate the effects caused by inaccurate terrain-corrected densities and fully take into account the rapid changes in depth within the craters, thereby improving the accuracy of the inversion model. We observed positive topographic elevations in specific localized sections, mainly the rims of large craters located in the lower left corners of areas A (Fig. 2a) and B (Fig. 3a). During the inversion process, the topographic elevations of these regions were set to zero, as their very limited area did not impact our inversion results.

Supplementary Fig. 9a, b shows the difference maps between the observed gravity anomalies (Figs. 2a and 3a) and the predicted gravity anomalies calculated from the 3D inversion models (Figs. 4 and 5). The difference maps show how well the inverted models resemble the observed gravity. The standard deviations of the differences for study areas A (Supplementary Fig. 9a) and B (Supplementary Fig. 9b) are 1.3 and 1.2 mGal, respectively, with a mean of 0.0 mGal. In Supplementary Fig. 9a, we observe short-wavelength fluctuations appearing as narrow strips, with widths of 6–10 km and lengths spanning 20–100 km, exhibiting a distinct S-N orientation. The short-wavelength jitters presented in Supplementary Fig. 9 black a clear directional pattern. The short-wavelength jitters are mainly related to the orbit-parallel striping in the GRAIL data at high degree17. These fluctuations are not anticipated to disrupt the primary magmatic structures in this study, as the latter possess scales significantly larger than those of the strips.

The inversion method introduced above was proposed in the Cartesian coordinate system and can be directly extended to the spherical coordinate system. However, the inversion process can be much more efficient in the Cartesian coordinate system, where the inversion space can be divided into a series of adjacent prisms of the same size, and the FFT technique can highly speed up the inversion process94. To investigate the effect of the lunar curvature for this study, we inverted the Bouguer gravity in study areas A and B in the spherical coordinate system separately (Supplementary Fig. 10a, b). The resulting density models show a high degree of similarity to the models obtained in the Cartesian coordinate system (Figs. 2c and 3c) in terms of the density contrasts and the spatial extent of the inverted anomalies. However, due to the utilization of finer cells (1.5 km \times 1.5 km \times 0.25 km), the model in the Cartesian coordinate system reveals more detailed structures than in the spherical coordinate system (0.15° \times 0.15° \times 0.25 km), in which the cell size was limited by computational performance. Therefore, this study presents and discusses the inversion results using the Cartesian coordinate system.

Several factors can influence inversion results when using the probabilistic method described above, including the parameters of the model covariance matrix (C_M), the depth weighting function, and the depth extent of the inversion domain. To assess the impact of these factors, we conducted a sensitivity analysis and compared the inversion results obtained with varying parameter settings (refer to the Supplementary for details; Supplementary Figs. 11, 12, and 13). These factors primarily have a modest impact on the thicknesses of the high-density bodies. The connectivity between the western PKT border intrusion and the northern and southern intrusions remains distinctly apparent. After comparing the inversion density models with previous studies on the northern and southern intrusions of the Marius Hills24,25 and the middle branch intrusion5 in study area B, through forward modeling, we selected the optimal models for study areas A and B shown in Figs. 2 and 3, as they demonstrate better agreement with these earlier studies. These models are derived using a \beta value of 1, an inversion domain depth extent of 20 km, and the model covariance parameters are listed in Table 1.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have combined gravity and topography data to explore the subsurface density structure of a major volcanic province on the Moon in an attempt to throw light on the volcanic plumbing system. Extremely interesting results are obtained that support and significantly expand on previousy-published ideas about the sub-surface volcanic systems in the areas studied.\nI am not an expert on the mathematical techniques used in the analysis. However, I am happy to comment on the authors' geophysical interpretation of their results, and find these to be entirely plausible.\nI have line-numbered a copy of the Word file of the text and linked the following questions and minor editorial suggestions to it.\nline 20. I suggest you remove the \"'s\" from the end of \"Procellarum's\".\nline 81. Why did you select the specific value 2,800 kg/m3 for the higher crustal density used to derive the Bouguer anomaly? And looking at Fig. 6a, how does the presence of the surface volcanics, presumably significantly denser than the anorthosite crust unless they were very vesicular, influence this part of the analysis?\nline 128. The word \"which\" makes no sense here and should be deleted.\nlines 131-140. The shortest dimension of the bodies in Area B is 5 to 7.5 km. Have you given any thought to the likely survival time against conductive cooling of magma bodies with these dimensions? A calculation based on\n(distance of travel of thermal wave) = ~2 x SQRT(thermal diffusivity x time)\nwould suggest that with a thermal diffusivity of say 5 times 10 to the power -7 meters-squared per second, a 5 km thick body\nwould cool significantly in ~400 thousand years.\nline 140. It is better to change \"as the range goes\" to \"toward the\".\nline 168. Use \"less than\" rather than \"under\".\nline 257. I suggest you add \"horizontal\" after \"long\" for clarity.", "reviewer_2": "This manuscript presents analyses of gravity data for the Marius Hills region of the Moon, one of the few large volcanic constructs on the Moon that is an area of high interest. Previous work found evidence for large intrusive bodies below the Marius Hills. Separately, a previous study identified a large set of linear gravity anomalies surrounding the nearside\nProcellarum KREEP terrane interpreted as the magma plumbing system of the nearside maria. This study makes the case\nthat the intrusions below the Marius hills are connected to the linear anomalies in the region as part of a more complex and\nintegrated magma system related to the magma plumbing of the maria and Marius hills.\nOverall, I find this to be a very interesting and well-written manuscript. This is a novel and potentially important interpretation\nthat is well supported by the manuscript. Overall, I find the connection between the Marius Hills volcanism and the PKT\nborder structures to be an interesting and important result. This connection clearly has important implications for the volcanic\nevolution of the Moon. However, the manuscript could do better in laying out exactly what the broader implications of this\nnew interpretation are. How does this new interpretation change our overall understanding of lunar volcanism?\nBelow I include a number of line-by-line comments for the authors to consider\n[82-86] I am glad to see you tested two density values for the Bouguer correction. Since your region of interest is within the\nmaria, a higher density would be appropriate (~3000 kg/m3; see Goossens 2022). This choice could have a noticeable\nimpact for your analysis of the Marius Hills themselves, which appear dominantly basaltic and do have considerable relief.\nThis is particularly true for the Marius NH anomaly, in which your Bouguer map shows a positive Bouguer anomaly beneath\nan area of positive topography - using an overly low density in the Bouguer correction contributes to this positive Bouguer\nanomaly. Given the relief of 2-3 km, a Bouguer correction density of 2500 rather than 3000 would contribute an extra 40-60\nmGal, which is not trivial in comparison to the observed value of 150 mGal in your Bouguer map. I would recommend\nfocusing on the higher density value that you tested (2800 rather than 2550).\n[90] It would be best not to justify the choice of density based on the density of feldspathic rocks given your study site in the\nmaria.\n[Figure 2] This is a matter of personal preference, but I find these sorts of perspective views draped on topography to\nobscure rather than reveal information. There is no new information in the topographic drape that is not also in a flat simple\ncolor map. The perspective adds distortion (including areas in which the data cannot be seen behind a topographic high).\nSimple map views as in Figure 4 are preferable. The perspective views of the isosurfaces in Figure 3 are a useful way to\nsupplement a simple map view (though figure 2c is in a different perspective than 2b, causing some confusion).\n[120] Please include a brief description of the methodology in the main text, which can refer to the supplement for more\ndetail. Please also include a brief statement/discussion of non-uniqueness in the main text, with more detail in the appendix.\nInversions of gravity data like this are fundamentally non-unique - this is not a weakness of this study which is using an\naccepted and rigorous approach for inverting the gravity data, but it should be acknowledged. This data could also be (and\nhas also been) inverted assuming that these same gravity anomalies arise from thickness variations in the maria, or\nvariations in dense layers below the crust. 3D inversions are particularly non-unique and rely on added assumptions and\nconstraints, and this requires more discussion in the manuscript and supplement. The inversions and interpretations in this\nstudy are sound and are very interesting, but some discussion of this non-uniqueness and alternative interpretations is\nneeded. It would also be good to include geological evidence supporting the inversions - for example, given the surface\nvolcanism associated with the marius hills and the strong anomalies that have been previously interpreted as sill-like\nintrusions, the interpretation of the other linear anomalies as large intrusive bodies is geologically reasonable.\n[147] “The densities of these intrusions lie between…” This statement depends on the density of the surrounding medium. At\nthese depth ranges, these intrusions are presumably surrounded by feldspathic crust (unlike the basaltic surface material,\nsee above), thus the lower density surrounding medium is reasonable. The reduction in porosity at depth may have an effect\nalso. Some discussion of these matters (perhaps in the supplement) is warranted.\n[145] “reported in the literature” Added discussion of how these inversion results compare with the previous simpler study\nwould help.\n[149] Please give the coordinates of the centers of these newly discovered features. Could these two “bowl-like” features\nactually be impact craters that were filled in by the mare volcanism (e.g., ref 15)?\n[152-155] The intrusion inferred beneath Reiner crater is interesting, but I have some concerns. In the Bouguer maps, it does\nnot appear that there is a positive Bouguer anomaly centered on the crater. Because of the different perspective view for the\ndensity model, it is difficult to see how the density anomaly interpreted as an intrusion compares to the crater. This is also an\narea where the Bouguer correction density matters - using too low of a density for the Bouguer correction will result in a relict\nnegative gravity anomaly over the crater itself (this can be seen in figure 2b), which will affect the inversion (though cannot\nexplain the positive density anomaly discussed here).\n[167] “To minimize the potential influence of crustal thinning…” I appreciate the inclusion of this section. In practice, any\ngravity anomalies below about degree 80 can be ascribed to crustal thickness variations, though they may also be explained\nby intrusive bodies. Some added discussion of this non-uniqueness and the merits of the different interpretations would be\nhelpful.\n[255-257] “However, the origin of the magmatic conduits, their role in the regional magmatic transport, and their impacton the\nvolcanic evolution in the Procellarum Oceanus region demand further inquiry.” This is a very important point, and the lack of\na more substantive discussion of this point is a shortcoming of this manuscript. The manuscript shows that the Marius Hills\nmay be connected to larger regional elongated anomalies that may represent the magma plumbing of the region - but what\ndoes it all mean? Previous studies interpreted the linear anomalies as related to the magma plumbing system of this region,\nand interpreted the Marius Hills gravity as related to subsurface intrusions. How does this new study connecting the two\nchange our understanding of the evolution of this region and the Moon? The manuscript could go further in addressing the\nbroader significance.\n[264-265] This discussion of vertical dikes vs lateral sills raises a question - previous work interpreted the linear structures\nas related to dike-like structures based on their linearity. In the new interpretation in this study, why are the elongated gravity\nanomalies that connect the MH anomalies to each other and to the PKT border structures linear? The gravity inversions are\nconsistent with sill-like intrusions, but those sill like intrusions are greatly elongated in one direction.\n[370] The methods describe a depth weighting as is commonly invoked in these types of models. Is wz the only component\nof the overall weighting matrix W? Other models have also included model smoothness criteria, and used a regularized\ninversion - is this also done in this model?\n[374-375] The choice of different values of beta (exponent in the depth weighting) for the two study areas seems odd. Please\nshow inversions for both values of beta for both areas and provide more discussion in the supplement.\nSome of the references are incomplete, giving only authors, title, and year without journal information.", "reviewer_3": "This study makes use of GRAIL data to investigate the subsurface magmatic structures of two regions within the Oceanus\nProcellarum: (1) Procellarum’s western border and (2) Northern and Southern Marius Hills. To achieve this objective, 3D\ninversion technique is employed to translate the GRAIL gravity data and construct detailed 3D density models of the\nsubsurface to understand the density distribution beneath the specified regions. Examining the 3D density models, authors\nsuggest an interconnected sill-like magmatic conduits’ network beneath the Oceanus Procellarum connecting the northern\nand southern halves of the Marius Hills as well as the border structures of the Procellarum region. The proposal of a ‘bridge’\nbetween the northern and southern halves of the Marius Hills was proposed by Kiefer (2013). This research, for the first time,\nprovides evidence for the connected subsurface structures not only for the MHVC but also for the adjacent regions.\nFurther, this finding of interlinked subsurface magmatic structures is also suggested to be accountable for the existing ages\nand compositional differences between the older domes of the MHVC and the younger mare flows in the region. This paper\nproposes that the younger mare flows of the MHVC could have been originated from the deep feeder dykes beneath the\nwestern border of the Procellarum terrain as these regions exhibit a connected subsurface magmatic structure in the derived\n3D density models.\nTwo previously undiscovered magmatic intrusions situated west of the northern Marius intrusion and to the southwest of the\nReiner crater have also been detected for the very first time from the derived density models of this paper. These newly\nidentified magmatic bodies are inferred to serve as potential filled and buried impact craters.\nOverall, this study and observations are good as they play a crucial role in understanding the magma plumbing system and\ntheir subsurface behaviour as the norm. It highlights the importance of looking at the regional geological context while\nexamining magmatism in a region on the Moon. However, the arguments presented in the study needs to be substantiated\nby the compositional and temporal relationships to build a complete picture of volcanism in the Oceanus Procellarum. In\nthis view, ‘Major revisions’ are suggested in the manuscript. I suggest the following comments be incorporated within the\nmanuscript:\n(1) The title of the manuscript implies that the magmatic structures and interconnections existing in the entire Oceanus\nProcellarum have been examined, which would have been a very significant contribution. However, the study emphasizes\nprimarily on the MHVC of the Oceanus Procellarum. Examination of the surrounding regions of Aristarchus, Flamsteed etc.\nis suggested to observe if such interconnections also exist in these areas.\n(2) The Results and Discussion section needs reorganization having separate section on ‘Study Area A’ and ‘Study Area B’\nfor clarity and emphasis on the importance of both these regions. Further, the importance of the ‘Study Area B’ needs to be\ndiscussed in detail.\n(3) Carry out a detailed mineralogical analysis of the region. This is essential to substantiate the proposed hypothesis of\ninterconnections, as shown in Fig. 6. The smaller craters in the area may provide important insights regarding the same.\n(4) The temporal relationships between the formation of the domes, surrounding basalts and the wrinkle ridges need to be\ndiscussed.\nAdditionally, the following concerns are needed to be addressed:\n(1) Authors are suggesting that the magma that accumulates beneath the craters need not to have deep crustal or mantle\norigin. Rather, they can migrate through the fractures present in the lunar crust. To make such a statement, it is necessary to\naddress the factors that control the magma transport at shallow crustal levels and the local structural features that exert\ncontrol over the magma migration.\n(2) Line 132, What is the level of uncertainties or errors?\n(3) Line 194, It is not clear from the description that how it can be inferred from the 3D inversion model that initially 3 discrete\nmagma chambers gradually interconnected with the progress in the filling processes. This needs to be elaborated.\n(4) Line 222-224, In order to establish that the buried magma intrusion in the region exerted structural controls on the wrinkle\nridges, it is important to constrain to timescale of these two processes. In case of a younger magmatic episode, it is also\npossible that the wrinkle ridges controlled the emplacement of the magmatic unit. Consider this, especially since this region\nhas experienced prolonged magmatic activity.\n(5) Lines 266-272, Is megaregolith the only cause for the formation of the interconnections? Do subsurface structural\nfeatures such as faults associated with the wrinkle ridges have a role to play?\n(6) Lines 272-273 It is mentioned that higher heat flux in the PKT region could have played a role in these magmatic\nactivities. How?\n(7) Line 322, Deep feeder dykes below the PKT border are not shown in the schematic fig. 6.\n(8) Line 323-325. It is essential to discuss the implications of the proposed scenario on the mineralogy and the surface age\nof the basalts in the region to substantiate the findings. I suggest a brief discussion on the same may be provided.\n(9) Line 333, “100nT” looks typo error; instead, it should be “100mGal”\n(10) The paper is generally well-written, presenting some interesting and significant results, but the novelty of the study\nneeds to be brought out and highlighted more clearly and prominently."}, "author_response": "Dear Reviewers,\nWe sincerely appreciate the positive feedback and evaluations from the reviewers regarding our\nresearch. Your insightful comments and suggestions have been invaluable in enhancing the quality\nof our manuscript, and we are grateful for your contributions to our work.\nWe have responded to each of the questions and made corresponding revisions and improvements\nto the manuscript. The following is a detailed response to all of the points raised during the review\nprocess. Please note that all changes to the manuscript can be readily seen in the marked-up version\nthat accompanies the updated submission. Unless otherwise stated, line number references below\nare in relation to the marked-up version of the manuscript (“manuscript_markedup.pdf”). The\nreviewer’s comments are in black, and our responses are in blue.\nWe hope the revisions meet the requirements. Once again, we greatly appreciate the constructive\nfeedback and look forward to providing a revised submission that aligns with the high standards\nof your journal.\nREVIEWER COMMENTS\n-------------------------------------- Reviewer #1 (Remarks to the Author) -----------------------------\nComments: The authors have combined gravity and topography data to explore the subsurface\ndensity structure of a major volcanic province on the Moon in an attempt to throw light on the\nvolcanic plumbing system. Extremely interesting results are obtained that support and significantly\nexpand on previously-published ideas about the sub-surface volcanic systems in the areas studied.\nResponse: Thank you very much for your positive feedback and interest in our research. We\nappreciate your valuable comments.\nComments: I am not an expert on the mathematical techniques used in the analysis. However, I\nam happy to comment on the authors' geophysical interpretation of their results, and find these to\nbe entirely plausible.\nResponse: Thank you very much for your feedback. We greatly appreciate your insights on\ngeophysical interpretation, which are highly valuable and have contributed to improving the\nquality of our work.\nComments: I have line-numbered a copy of the Word file of the text and linked the following\nquestions and minor editorial suggestions to it.\nResponse: Thank you for your detailed suggestions on the text modifications. We have accepted\nthem and made the corresponding changes.\nComments: line 20. I suggest you remove the \"'s\" from the end of \"Procellarum's\".\nResponse: It was removed accordingly.\nComments: line 81. Why did you select the specific value 2,800 kg/m3 for the higher crustal\ndensity used to derive the Bouguer anomaly? And looking at Fig. 6a, how does the presence of the\ns of the\nsurface volcanics, presumably significantly denser than the anorthosite crust unless they were very\nvesicular, influence this part of the analysis?\nResponse: Thank you for your valuable comments. In previous studies of the Oceanus Procellarum,\na crust density of 2,550 kg/m³ was used for Bouguer correction when analyzing the border linear\nstructures (Andrews-Hanna et al., 2014). A higher crust density of 2,800 kg/m³ was applied for\nBouguer correction in the study of four positive Bouguer gravity anomalies, which are located\nsouth of the Aristarchus Plateau, north of Flamsteed crater, and within the Marius Hills volcanic\ncomplex, both in the northern and southern regions (Deutsch et al., 2019). Therefore, the 2,560\nkg/m³ and 2,800 kg/m³ density values were both used in this manuscript. In the revised manuscript,\nwe carefully compared the Bouguer gravity anomalies derived with 2,560 kg/m³, 2,800 kg/m³, and\n3,000 kg/m3 and corresponding inversion results (refer to Supplementary Figs. 2-4). As you can\nsee, the density models derived with corrections of 2,560, 2,800, and 3,000 kg/m³ show high\nconsistency. The impact of Bouguer correction using different densities is most pronounced in\nregions with significant terrain variations, such as the Reiner crater and Marius-NH\n(Supplementary Figs. 3a, 4a). However, these changes are localized and do not alter the main\nstructural features or the observed magmatic connectivity.\nWe added text to address the impact of the Bouguer correction with different crust densities in the\nrevised manuscript; please refer to Lines 172-189.\nComments: line 128. The word \"which\" makes no sense here and should be deleted.\nResponse: It was deleted accordingly.\nComments: lines 131-140. The shortest dimension of the bodies in Area B is 5 to 7.5 km. Have\nyou given any thought to the likely survival time against conductive cooling of magma bodies with\nthese dimensions? A calculation based on (distance of travel of thermal wave) = ~2 x\nSQRT(thermal diffusivity x time) would suggest that with a thermal diffusivity of say 5 times 10\nto the power -7 meters-squared per second, a 5 km thick body would cool significantly in ~400\nthousand years.\nResponse: Thank you very much for providing the idea of estimating cooling time for magma\nbodies. This could make for a fascinating study, especially in the context of the Moon's\nProcellarum KREEP Terrane (PKT). The simple model used here assumes that cooling occurs\nsolely through thermal conduction, disregarding other heat transfer mechanisms such as\nconvection and radiation. To make a more reasonable estimation of the cooling time of magma\nbodies, other factors like the initial temperature of the magma, ambient temperature, size and shape\nof the magma body, magma composition, crystallization during cooling, and pressure conditions\nneed to be considered. PKT is a geochemically distinct region on the Moon, characterized by a\nthin lunar crust and a high concentration of radiogenic heat-producing elements that are considered\nto have driven the thermal evolution of the lunar mantle in this region. How these factors affect\nmagma evolution and cooling in the PKT region needs careful simulation. Therefore, using this\nsimple model to estimate cooling times presents many challenges. For these reasons, we did not\ndiscuss the survival time of the magmatic bodies. However, our 3D magmatic modelswill provide\nvaluable reference and support for future studies on magma evolution. Thank you again for\nproposing this interesting comment.\nComments: line 140. It is better to change \"as the range goes\" to \"toward the\".\nResponse: It was revised accordingly.\nComments: line 168. Use \"less than\" rather than \"under\".\nResponse: It was revised accordingly.\nComments: line 257. I suggest you add \"horizontal\" after \"long\" for clarity.\nResponse: It was added accordingly.\n---------------------------------- Reviewer #2 (Remarks to the Author) ---------------------------------\nComments: This manuscript presents analyses of gravity data for the Marius Hills region of the\nMoon, one of the few large volcanic constructs on the Moon that is an area of high interest.\nPrevious work found evidence for large intrusive bodies below the Marius Hills. Separately, a\nprevious study identified a large set of linear gravity anomalies surrounding the nearside\nProcellarum KREEP terrane interpreted as the magma plumbing system of the nearside maria. This\nstudy makes the case that the intrusions below the Marius hills are connected to the linear\nanomalies in the region as part of a more complex and integrated magma system related to the\nmagma plumbing of the maria and Marius hills.\nOverall, I find this to be a very interesting and well-written manuscript. This is a novel and\npotentially important interpretation that is well supported by the manuscript. Overall, I find the\nconnection between the Marius Hills volcanism and the PKT border structures to be an interesting\nand important result. This connection clearly has important implications for the volcanic evolution\nof the Moon. However, the manuscript could do better in laying out exactly what the broader\nimplications of this new interpretation are. How does this new interpretation change our overall\nunderstanding of lunar volcanism?\nResponse: We appreciate your encouraging feedback and interest in our research. We also thank\nyou for your valuable comments. Regarding your suggestion to better elaborate on the broader\nimplications of our interpretation, we agree that this is an important aspect to address.\nWe have expanded our discussion to address the broader implications of our findings in the\nOceanus Procellarum region (see Lines 438-446). Additionally, the abstract has been revised to\nreflect these points. In particular, we emphasize the crucial role that lateral magma transport\nprocesses play in shaping the volcanic evolution of the Moon. Understanding these lateral transport\nmechanisms is crucial for advancing our knowledge of lunar volcanic processes and should be an\nimportant consideration in future research efforts.\nBelow I include a number of line-by-line comments for the authors to consider.\nComments: [82-86] I am glad to see you tested two density values for the Bouguer correction.\nSince your region of interest is within the maria, a higher density would be appropriate (~3000\nkg/m3; see Goossens 2022). This choice could have a noticeable impact for your analysis of the\nMarius Hills themselves, which appear dominantly basaltic and do have considerable relief. This\nis particularly true for the Marius NH anomaly, in which your Bouguer map shows a positive\nBouguer anomaly beneath an area of positive topography - using an overly low density in the\nBouguer correction contributes to this positive Bouguer anomaly. Given the relief of 2-3 km, a\nBouguer correction density of 2500 rather than 3000 would contribute an extra 40-60 mGal, which\nis not trivial in comparison to the observed value of 150 mGal in your Bouguer map. I would\nrecommend focusing on the higher density value that you tested (2800 rather than 2550).\nResponse: Thank you for your insightful comments. In the revised manuscript, we have\nincorporated your suggestion by using the Bouguer gravity anomalies corrected with a density of\n2800 kg/m³, along with the corresponding inversion results. Moreover, we compared the Bouguer\ngravity anomalies derived with 2,560 kg/m³, 2,800 kg/m³, and 3,000 kg/m3 and corresponding\ninversion results (refer to Supplementary Fig. 2). We also conducted detailed profile comparisons\nof these Bouguer gravity anomalies in the Marius-NH, Marius-SH, Reiner crater, and western PKT\nborder regions, analyzing how variations in the crust density used for Bouguer correction affected\nthese anomalies and the corresponding inversion results (see Supplementary Figs. 3 and 4).\nWe added text to address the impact of the Bouguer correction with different crust densities in\nthe revised manuscript; please refer to Lines 172-189.\nComments: [90] It would be best not to justify the choice of density based on the density of\nfeldspathic rocks given your study site in the maria.\nResponse: Thank you for the comment. We have removed the related text.\nComments: [Figure 2] This is a matter of personal preference, but I find these sorts of perspective\nviews draped on topography to obscure rather than reveal information. There is no new information\nin the topographic drape that is not also in a flat simple color map. The perspective adds distortion\n(including areas in which the data cannot be seen behind a topographic high). Simple map views\nas in Figure 4 are preferable. The perspective views of the isosurfaces in Figure 3 are a useful way\nto supplement a simple map view (though figure 2c is in a different perspective than 2b, causing\nsome confusion).\nResponse: Thank you for the comments. We have adjusted the perspective view of Figure 2 so\nthat Figures 2a, 2b and 2c are in the same perspective. Regarding the issue of model data being\nobscured by topographic features, we increased the transparency of the topographic map, which\nhas improved the visibility of the model data. Please refer to the new Figure 2 to check if it meets\nthe requirements.\nComments: [120] Please include a brief description of the methodology in the main text, which\ncan refer to the supplement for more detail. Please also include a brief statement/discussion of non-\nuniqueness in the main text, with more detail in the appendix. Inversions of gravity data like this\nare fundamentally non-unique - this is not a weakness of this study which is using an accepted and\nrigorous approach for inverting the gravity data, but it should be acknowledged. This data could\nalso be (and has also been) inverted assuming that these same gravity anomalies arise from\nthickness variations in the maria, or variations in dense layers below the crust. 3D inversions are\nparticularly non-unique and rely on added assumptions and constraints, and this requires more\ndiscussion in the manuscript and supplement. The inversions and interpretations in this study are\nsound and are very interesting, but some discussion of this non-uniqueness and alternative\ninterpretations is needed. It would also be good to include geological evidence supporting the\ninversions - for example, given the surface volcanism associated with the marius hills and the\nstrong anomalies that have been previously interpreted as sill-like intrusions, the interpretation of\nthe other linear anomalies as large intrusive bodies is geologically reasonable.\nResponse: Thank you for your thoughtful comments and suggestions. We made corresponding\nrevisions based on your comments and suggestions, which have greatly improved the quality of\nthe manuscript.\nWe have included a brief description of the methodology in the main text (refer to Lines 109-118).\nWe added a concise discussion on the non-uniqueness of gravity data inversions in the main text\n(refer to Lines 114-118), with more detailed exploration provided in the section on Sensitivity\nanalysis of inversion in Method (refer to Lines 520-560) where we tested and discussed the crucial\nfactors that can impact the inversion results. The inversion models derived using different\nconstraints, depth weighting functions, and inversion domains are shown in Supplementary Figs.\n11-13).\nWe also incorporated statements on the geological evidence to support the plausibility of the\ninterpretation of linear anomalies as large intrusions (refer to Lines 275-280). Thank you again for\nyour valuable insights.\nComments: [147] “The densities of these intrusions lie between…” This statement depends on\nthe density of the surrounding medium. At these depth ranges, these intrusions are presumably\nsurrounded by feldspathic crust (unlike the basaltic surface material, see above), thus the lower\ndensity surrounding medium is reasonable. The reduction in porosity at depth may have an effect\nalso. Some discussion of these matters (perhaps in the supplement) is warranted.\nResponse: Thank you for your insightful comments. We agree that the density of the surrounding\nmedium, particularly the feldspathic crust, plays a critical role in interpreting the densities of the\nintrusions. Therefore, more factors needed to be considered and discussed when using absolute\ndensity to interpret the magmatic bodies. In the revised manuscript, we used density contrast\ninstead of absolute density to show and interpret magmatic intrusions. This will simplify the\ndensity interpretation. This is also consistent with the way magmatic bodies have been interpreted\nin previous literature (e.g., Andrews-Hanna et al., 2014; Deutsch et al., 2019; Kiefer, 2013). The\neffect of the reduction in porosity at depth was discussed when we considered the limited depth of\nthe inversion models (refer to Lines 281-301).\nComments: [145] “reported in the literature” Added discussion of how these inversion results\ncompare with the previous simpler study would help.\nResponse: Thank you for your comments. These inversion results were carefully compared with\nthe previous simpler study in the Discussion (refer to Lines 238-251).\nComments: [149] Please give the coordinates of the centers of these newly discovered features.\nCould these two “bowl-like” features actually be impact craters that were filled in by the mare\nvolcanism (e.g., ref 15)?\nResponse: Thank you for your suggestion. We added the coordinates for the three inferred buried\ncraters identified from our inversion models.\nEvans et al. (2016) (ref. 15) analyzed GRAIL Bouguer gravity anomalies and proposed that the\ncircular features at Marius Hills may be buried impact craters, caused by the density contrast\nbetween the feldspathic crust and mare basalt infill. Zhang et al. (2018) suggested that volcanic\nfeatures in the area could be intruded sills beneath craters and concluded that impact craters may\ncreate crustal weaknesses favorable for sill and dike intrusions. Our inversion models are\nconsistent with these findings. Besides the two “bowl-like” features in Marius Hills, intrusions\nwith various volumes have been discovered beneath many craters, including but not limited to the\ninferred filled and buried crater adjacent to the Marius-NH (Fig. 2d), Reiner crater (Fig. 2e),\nMarius crater (Fig. 2d), Briggs B crater (Fig. 3d), Lichtenberg crater (Fig. 3d), Seleucus crater (Fig.\n3e), a potential filled and buried impact crater (located at 55.7 W, 3.6 N; see Figs. 4b, 4c and 4i)\nsituated to the southwest of the Reiner crater. We have added text to address this point, please see\nLines 251-255.\nComments: [152-155] The intrusion inferred beneath Reiner crater is interesting, but I have some\nconcerns. In the Bouguer maps, it does not appear that there is a positive Bouguer anomaly\ncentered on the crater. Because of the different perspective view for the density model, it is difficult\nto see how the density anomaly interpreted as an intrusion compares to the crater. This is also an\narea where the Bouguer correction density matters - using too low of a density for the Bouguer\ncorrection will result in a relict negative gravity anomaly over the crater itself (this can be seen in\nfigure 2b), which will affect the inversion (though cannot explain the positive density anomaly\ndiscussed here).\nResponse: Thank you for your insightful comments. In the revised manuscript, we carefully\ncompared the Bouguer gravity anomalies across the Reiner crater using different crust densities\nfor Bouguer correction (2,560, 2,800, 3,000 kg/m3; see Supplementary Fig. 3a) and corresponding\ninversion results (see Supplementary Figs. 3b-3d). The higher density of 3,000 kg/m³ yields the\nbest-fit Bouguer correction for the Reiner crater (Supplementary Fig. 3a), which is similar to the\nmean density estimated for the mare (Neumann et al., 2018) (~3090 kg/m³) and the bulk density\nestimated for the Marius Hills surface basalts (3150 kg/m3) based on observations from the Lunar\nProspector Gamma Ray Spectrometer (Kiefer et al., 2012; Prettyman et al., 2006). However, it\nalso significantly enhances the western belt gravity anomalies (Supplementary Fig. 3a) due to the\nrelatively low topography in that area. Ultimately, we opted for a density of 2,800 kg/m³ as a more\nrepresentative choice for the entire region.\nIt appears that with the increasing Bouguer gravity amplitude over the Reiner crater, it, along with\nthe surrounding gravity anomaly signals, form longer-wavelength gravity signals (see\nSupplementary Fig. 3a). Consequently, the inverted density model shows an increase in the density\ncontrast at depth while the density contrast of the magmatic intrusion beneath the Reiner crater\ndecreases slightly (see Supplementary Figs. 3b-3d). However, the connectivity between the PKT\nborder structures and the Marius-SH remains unchanged.\nIn the density horizontal slice maps (Fig. 4a-4c), we added black circles to indicate the location of\nthe Reiner crater and Marius crater so that the spatial relationship between craters and intrusions\ncan be observed more clearly. As can be seen from Fig. 4, the intrusion is located directly beneath\nthe Reiner crater at a shallow depth (3.5 km; see Fig. 4a) and slightly tilts southwestward with\nincreasing depth (5.5 km; see Fig. 4b).\nComments: [167] “To minimize the potential influence of crustal thinning…” I appreciate the\ninclusion of this section. In practice, any gravity anomalies below about degree 80 can be ascribed\nto crustal thickness variations, though they may also be explained by intrusive bodies. Some added\ndiscussion of this non-uniqueness and the merits of the different interpretations would be helpful.\nResponse: Thank you for your valuable suggestions. To address the non-uniqueness inherent in\nthe inversion process, we performed a sensitivity analysis on several critical parameters that could\naffect the inversion results, including the covariance matrix parameter, the depth weighting\nfunction, and the depth extent of the inversion domain (see Method for more details). The inversion\nmodels underwent slight modifications with changes in the settings, particularly affecting the\ndensity distribution in the depth direction due to the low resolution of the gravity data. This resulted\nin alterations to the size of the magmatic intrusions and the values of density contrast. This reflects\nthe inherent non-uniqueness of gravity data inversion, where anomalous bodies of different sizes\nand densities can produce the same gravity anomalies. However, as evident, the main structures of\nthe magmatic bodies remain consistent with one another. Moreover, the model of crustal thickness\nvariations was discussed in Lines 275-280. Given the potential magma chamber compartments\nrevealed by the 3D inversion model, as well as intrusions identified beneath the Marius Hills that\ncould have contributed to facilitating surface volcanism, the interpretation of linear anomalies as\nlarge intrusions appears to be geologically reasonable.\nComments: [255-257] “However, the origin of the magmatic conduits, their role in the regional\nmagmatic transport, and their impact on the volcanic evolution in the Procellarum Oceanus region\ndemand further inquiry.” This is a very important point, and the lack of a more substantive\ndiscussion of this point is a shortcoming of this manuscript. The manuscript shows that the Marius\nHills may be connected to larger regional elongated anomalies that may represent the magma\nplumbing of the region - but what does it all mean? Previous studies interpreted the linear\nanomalies as related to the magma plumbing system of this region, and interpreted the Marius\nHills gravity as related to subsurface intrusions. How does this new study connecting the two\nchange our understanding of the evolution of this region and the Moon? The manuscript could go\nfurther in addressing the broader significance.\nResponse: Thank you for your insightful comments and suggestions. The discussion primarily\nfocuses on these issues. We have revised the article to make our arguments clearer. We explain the\norigin and formation of the magmatic conduits identified by drawing comparisons to long-distance\nlateral magma transport mechanisms observed on Earth. Pressure gradients caused by the\nprolonged and extensive volcanic activity in the Marius Hills region compared to the plumbing\nsystem beneath the PKTborder (Lines 347-361) and horizontal pressure gradients along the PKT\nrift (Lines 362-371) can account for the origin and formation process. Moreover, fragmented Lunar\nMegaregolith (Lines 340-346), reactivation of pre-existing fractures (Lines 357-361), impact-\ninduced fractures (Lines 373-378), and high heat flow (Lines 390-392) also contribute to this\nprocess. The magmatic conduits between the PKT border structures and the Marius Hills intrusions\naccount for the compositional distinction between the domes and surrounding and embaying mare\nin the Marius region, which might also contribute to the floor-filling process in the Marius crater\n(Lines 393-415), highlighting the role of the magmatic conduits in their impact on the volcanic\nevolution in the Procellarum Oceanus region.\nWe expanded the discussion to cover the broader implications of our findings in the Oceanus\nProcellarum region (see Lines 440-446) and updated the abstract to reflect these changes. We\nparticularly highlight the critical role of lateral magma transport in the Moon's volcanic evolution,\nwhich was not discovered before. Understanding these lateral transport mechanisms is crucial for\nadvancing our knowledge of lunar volcanic processes and should be an important consideration in\nfuture research efforts.\nComments: [264-265] This discussion of vertical dikes vs lateral sills raises a question - previous\nwork interpreted the linear structures as related to dike-like structures based on their linearity. In\nthe new interpretation in this study, why are the elongated gravity anomalies that connect the MH\nanomalies to each other and to the PKT border structures linear? The gravity inversions are\nconsistent with sill-like intrusions, but those sill like intrusions are greatly elongated in one\ndirection.\nResponse: Thank you for your comments. The magmatic conduits resulting from long-distance\nlateral transport on Earth commonly exhibit linear characteristics; please refer to the literature (e.g.,\nGudmundsson et al., 2016; Ishizuka et al., 2017; Sigmundsson et al., 2015; Wilding et al., 2023).\nRegarding the linearity exhibited by the magmatic conduits connecting the Marius-SH, Marius-\nNH intrusions, and PKT border structures, as we explained in the Discussion, this is mainly\nassociated with the gradient pressure and impact-induced fractures. The prolonged and extensive\nvolcanic activity in the Marius Hills region could have caused pressure gradients between the\nplumbing system beneath the PKT border and the Marius Hills region. These pressure gradients\ncould have facilitated rapid magma propagation through the fractures from the PKT border to the\nMarius-SH and Marius-NH areas, resulting in linear magmatic conduit structures between these\nvolcanic systems. Moreover, the magmatic conduit connecting Marius-SH and Marius-NH shows\na NNW–SSE trend that aligns with the regional trend in the Marius Hills region, which is\nconsidered to be related to pre-existing weak zones/fractures (Arya et al., 2018; Raitala, 1980).\nPlease refer to the Discussion for details.\nComments: [370] The methods describe a depth weighting as is commonly invoked in these types\nof models. Is wz the only component of the overall weighting matrix W? Other models have also\nincluded model smoothness criteria, and used a regularized inversion - is this also done in this\nmodel?\nResponse: Thank you for your comments on the Method. The is the only\ncomponent of the overall weighting matrix. 𝛽\n𝑤(𝑧) = 1/(𝑧max−𝑧)\nThe following shows our objective function,\n(1)\n𝑇\n𝑇 −1 −1\nThe seco𝚽nd= te(r𝐆m𝐦 co−rre𝐝sopbosn)d𝐂s Dto (th𝐆e𝐦 sm−oo𝐝tohbnse)s+s c(o𝐦nst−ra𝐦inta purs)ed𝐂 inM th(𝐦e tr−ad𝐦itiaopnra)l. inversion method\n(e.g., Li and Oldenburg, 1998). The traditional inversion method usually applies the derivatives of\nthe density model in three directions to control the smoothness of the inversion. In contrast, the\nmodel covariance matrix C in equation (1) is established using geostatistical theory (Geng et al.,\nM\n2019a; Geng et al., 2019b). An advantage of this inversion method is that by using the model\ncovariance matrix C calculated by geostatistical theory, a priori information, including physical\nM\nproperties, spatial extent, depth, and orientation of the source bodies, can be incorporated to\nconstrain the inversion and achieve an accurate and reliable solution. This can be seen in Table 1,\nwhich shows the parameters of the model covariance matrix. Moreover, the parameter of the Cd\ncan adjust the contribution of the first term (the data misfit function) and the constraint term (the\nsecond term). Please refer to Method for more details.\nComments: [374-375] The choice of different values of beta (exponent in the depth weighting)\nfor the two study areas seems odd. Please show inversions for both values of beta for both areas\nand provide more discussion in the supplement.\nResponse: Thank you for your comments on the Method. We compared the two inversion results\nby using and for study areas A and B (see Supplementary Fig. 12) in the Sensitivity\nanalysis. As you can see, the results show that densities with β=2 slightly converge at a shallower\ndepth than𝛽 t=ho1se with𝛽 β==12. In the revised manuscript, we chose the results with β=1 for both study\nareas A and B. These results still show good consistency with the results from previous studies\nusing forward modeling.\nComments: Some of the references are incomplete, giving only authors, title, and year without\njournal information.\nResponse: Thank you for your feedback. We corrected the references and ensured that all entries\nincluded complete information, such as the journal name and other relevant details.\n------------------------------------ Reviewer #3 (Remarks to the Author) -------------------------------\nComments: This study makes use of GRAIL data to investigate the subsurface magmatic\nstructures of two regions within the Oceanus Procellarum: (1) Procellarum’s western border and\n(2) Northern and Southern Marius Hills. To achieve this objective, 3D inversion technique is\nemployed to translate the GRAIL gravity data and construct detailed 3D density models of the\nsubsurface to understand the density distribution beneath the specified regions. Examining the 3D\ndensity models, authors suggest an interconnected sill-like magmatic conduits’ network beneath\nthe Oceanus Procellarum connecting the northern and southern halves of the Marius Hills as well\nas the border structures of the Procellarum region. The proposal of a ‘bridge’ between the northern\nand southern halves of the Marius Hills was proposed by Kiefer (2013). This research, for the first\ntime, provides evidence for the connected subsurface structures not only for the MHVC but also\nfor the adjacent regions.\nFurther, this finding of interlinked subsurface magmatic structures is also suggested to be\naccountable for the existing ages and compositional differences between the older domes of the\nMHVC and the younger mare flows in the region. This paper proposes that the younger mare flows\nof the MHVC could have been originated from the deep feeder dykes beneath the western border\nof the Procellarum terrain as these regions exhibit a connected subsurface magmatic structure in\nthe derived 3D density models.\nTwo previously undiscovered magmatic intrusions situated west of the northern Marius intrusion\nand to the southwest of the Reiner crater have also been detected for the very first time from the\nderived density models of this paper. These newly identified magmatic bodies are inferred to serve\nas potential filled and buried impact craters.\nOverall, this study and observations are good as they play a crucial role in understanding the\nmagma plumbing system and their subsurface behaviour as the norm. It highlights the importance\nof looking at the regional geological context while examining magmatism in a region on the Moon.\nHowever, the arguments presented in the study needs to be substantiated by the compositional and\ntemporal relationships to build a complete picture of volcanism in the Oceanus Procellarum. In\nthis view, ‘Major revisions’ are suggested in the manuscript. I suggest the following comments be\nincorporated within the manuscript:\nResponses: Thank you for your thorough and insightful comments on our manuscript. We greatly\nappreciate your recognition of the study’s contribution to understanding the subsurface magmatic\nstructures of the Oceanus Procellarum. Your comments are invaluable in helping us refine our\nwork. Thank you again for your constructive comments and support.\nFollowing your suggestions, we incorporated additional analysis of the mineralogical analysis\n(Marius Hills, refer to Lines 393-410; Lichtenberg crater and surrounding area; refer to Lines 423-\n432) and temporal relationships between the formation of the Marius Hill volcanic complex,\nwestern PKT border structures and the wrinkle ridges (refer to Lines 302-326). Additionally, we\nhave provided detailed responses to the comments, which you can refer to in the following sections\nof comments and responses.\nAdditionally, we have added a relevant introduction in the revised manuscript regarding the\nproposal of the \"bridge\" between the northern and southern halves of the Marius Hills, which was\nfirst mentioned in Kiefer's (2013) study. We sincerely thank you for reminding us.\nComments: (1) The title of the manuscript implies that the magmatic structures and\ninterconnections existing in the entire Oceanus Procellarum have been examined, which would\nhave been a very significant contribution. However, the study emphasizes primarily on the MHVC\nof the Oceanus Procellarum. Examination of the surrounding regions of Aristarchus, Flamsteed\netc. is suggested to observe if such interconnections also exist in these areas.\nResponse: Thank you for your insightful comments and suggestions on our manuscript. In the\nrevised manuscript, we added the inversion results in Flamsteed and Aristarchus as well as the\nsurrounding areas (see Supplementary Figs. 6 and 7). For the relevant introduction and discussion\nregarding these two regions, please refer to Lines 202-220.\nHere we briefly introduce the main discovery in these two regions. Our density model\n(Supplementary Fig. 6) shows the northern Flamsteed magmatic intrusion has a thickness of ~11.5\nkm and a diameter of ~150 km. A narrow (~20 km wide), sill-like corridor may connect the\nnorthern and southern Flamsteed intrusions, aligned with surface wrinkle ridges. However, further\ninvestigation is needed to confirm this connection due to the corridor's small width and potential\nerror from using the 160 kg/m³ iso-surface. The southern Aristarchus plateau intrusion\n(Supplementary Fig. 7) has a quasi-circular shape with a ~100 km diameter and ~6.5 km thickness.\nNeither the Flamsteed nor the Aristarchus intrusions appear connected to the western PKT border\nintrusions. We carefully examined the Bouguer gravity anomaly in the Oceanus Procellarum\nregion and did not find more evidence of connections between the western PKT border structures\nand other volcanic systems nor between the different volcanic systems themselves. As a\nconsequence, our discussion focuses on the inversion models (Figs. 2-5) in study areas A and B.\nComments: (2) The Results and Discussion section needs reorganization having separate section\non ‘Study Area A’ and ‘Study Area B’ for clarity and emphasis on the importance of both these\nregions. Further, the importance of the ‘Study Area B’ needs to be discussed in detail.\nResponse: Thank you for your valuable feedback. In the Results and Discussion section, we\nintroduce and discuss PKT border structures as a whole. Study areas A and B contain the western\nborder structures and the northwestern three branches, respectively, which are not suitable for a\nseparate discussion. Moreover, since the important interconnected magmatic network between\nMarius-NH, Marius-SH, and their connection to the PKT border structure lies within Study area\nA, the discussion section of our paper devotes considerable attention to explaining the potential\ntransport mechanism and favorable evidence for this connectivity. If we were to separate the\ndiscussion into distinct sections for areas A and B, the paper would become unbalanced.\nThat said, study area B remains significant. Alongside study area A, it presents a more complete\n3D subsurface model of the magmatic structure along the western PKT border, providing valuable\ninsights for the PKT border plumbing systems. Additionally, intrusions of varying sizes have been\ndiscovered beneath several craters in study area B, such as Briggs B crater (Fig. 3d), Lichtenberg\ncrater (Fig. 3d), and Seleucus crater (Fig. 3e). These findings suggest that magma accumulation\nbeneath craters is likely a common process. Importantly, the source of this magma may not always\nbe from the deep lunar mantle; rather, it could have migrated laterally through fractures in the lunar\ncrust before accumulating beneath impact craters. These discoveries, along with the observed\nconnections between volcanic systems in Oceanus Procellarum, underscore the critical role of\nlateral magma transport in the Moon's volcanic evolution (refer to Lines 433-449).\nComments: (3) Carry out a detailed mineralogical analysis of the region. This is essential to\nsubstantiate the proposed hypothesis of interconnections, as shown in Fig. 6. The smaller craters\nin the area may provide important insights regarding the same.\nResponse: Thank you for your insightful suggestions. In response to your suggestion, we have\nincorporated a brief mineralogical analysis of both the Marius Hills region (Lines 393-415) and\nthe Lichtenberg crater area (Lines 423-432). These analyses are based on previously published\nmineralogical analyses through spectral data and are used to support and validate our 3D inversion\nresults and findings.\nMost of the mare basalts had not been visited by previous missions. Thus, the understanding relied\nmainly on remote sensing techniques, which, however, lack ground truth verification. A detailed\nmineralogical study of the Oceanus Procellarum region would require extensive and rigorous\nresearch, primarily involving remote sensing and lunar rock sampling analysis, which goes beyond\nthe scope of our current study. Our primary focus remains on interpreting the deep magmatic\nstructures of the Moon using GRAIL gravity data. We hope that our findings will contribute to\nfuture research on lunar volcanic evolution as well as mineralogical and compositional analyses.\nComments: (4) The temporal relationships between the formation of the domes, surrounding\nbasalts, and the wrinkle ridges need to be discussed.\nResponse: Thank you for your insightful suggestions. We added a paragraph to address the\ntemporal relationships between the formation of PKT border structures, Marius-SH, surrounding\nbasalts, and wrinkle ridges. Please refer to Lines 302-326.\nAdditionally, the following concerns are needed to be addressed:\nComments: Authors are suggesting that the magma that accumulates beneath the craters need not\nto have deep crustal or mantle origin. Rather, they can migrate through the fractures present in the\nlunar crust. To make such a statement, it is necessary to address the factors that control the magma\ntransport at shallow crustal levels and the local structural features that exert control over the\nmagma migration.\nResponse: Thank you for your suggestions. The factors that control the magma transport at\nshallow crustal levels and the local structural features that exert control over the magma migration\nhave been detailed and discussed in the manuscript.\n1. Fragmented Lunar Megaregolith (Lines 340-346)\n2. Reactivation of Pre-existing Fractures (Lines 358-361).\n3. Systemwide Pressure Weaknesses, which can propagate quickly through the magmatic web,\nfacilitating lateral magma transport in different volcanic systems (Lines 347-357).\n4. The horizontal pressure gradient along the rift, which exerts a local structure control for\nthe magma migration from the PKT border to Marius-SH and Marius-NH (Lines 362-376)\n5. Impact-induced Fractures (Lines 376-378)\n6. High heat flow (Lines 389-392).\nComments: Line 132, What is the level of uncertainties or errors?\nResponse: Thank you for your question. We added the following text to address this point. “When\nusing the iso-surfaces of 150 kg/m3 and 170 kg/m3 as the interpreted magmatic bodies, the total\nvolume increases by 5% and decreases by 6%, respectively, compared to the 160 kg/m3 iso-\nsurfaces” (refer to Lines 134-136).\nComments: Line 194, It is not clear from the description that how it can be inferred from the 3D\ninversion model that initially 3 discrete magma chambers gradually interconnected with the\nprogress in the filling processes. This needs to be elaborated.\nResponse: Thank you for your comment. We added the Ⅰ, Ⅱ, Ⅲ to indicate the locations of the\npotential 3 discrete magma chambers inferred from our density model (see Figs. 3 and 5). The\nspatial architecture revealed by the density model shows that the middle branch is likely to have\nformed from three discrete magma chamber compartments (indicated by Ⅰ, Ⅱ, Ⅲ in Figs. 3c, 5a,\nand 5b). The three inferred magma chamber compartments gradually interconnected with each\nother, possibly throughout the filling process. This procedure is similar to the formation style of\nthe Great Dyke (Wilson, 1996).\nComments: Line 222-224, In order to establish that the buried magma intrusion in the region\nexerted structural controls on the wrinkle ridges, it is important to constrain to timescale of these\ntwo processes. In case of a younger magmatic episode, it is also possible that the wrinkle ridges\ncontrolled the emplacement of the magmatic unit. Consider this, especially since this region\nhas experienced prolonged magmatic activity.\nResponse: Thank you for your insightful comments and suggestions. We added a paragraph to\naddress the temporal relationships between the formation of PKT border structures, Marius-SH,\nsurrounding basalts, and wrinkle ridges. Please refer to Lines 302-326.\nThe faults associated with the wrinkle ridges could facilitate structural adjustments potentially\ntriggering volcanic eruptions in the Marius-SH, as discussed by Zhang et al. (2023). Except for\nthis, we have not found any evidence to support the view that the wrinkle ridges currently control\nthe emplacement of the magmatic unit in the PKT border structures or the Marius SH. The\ninversion models support the view (Andrews-Hanna et al., 2014; Zhang et al., 2023) that buried\nmagma intrusion in the region exerted structural controls on the wrinkle ridges (refer to Lines 302-\n326). However, we will continue to investigate the relationship between the wrinkle ridges and the\nemplacement of the magmatic unit. Thank you again for your insightful comments.\nComments: Lines 266-272, Is megaregolith the only cause for the formation of the\ninterconnections? Do subsurface structural features such as faults associated with the wrinkle\nridges have a role to play?\nResponse: Thank you for your questions. As we addressed above, megaregolith is not the only\ncause for the formation of the interconnections. Other factors include reactivation of pre-existing\nfractures, systemwide pressure weaknesses, horizontal pressure gradient along the rift, impact-\ninduced fractures and high heat flow.\nThe faults associated with the wrinkle ridges could facilitate structural adjustments potentially\ntriggering volcanic eruptions in the Marius-SH. We have added this point in the revised manuscript\n(refer to Lines 318-324).\nComments: Lines 272-273 It is mentioned that higher heat flux in the PKT region could have\nplayed a role in these magmatic activities. How?\nResponse: Thank you for your comments. We added the following text to address this point in\nthe revised manuscript below (Lines 389-392).\n“Higher-than-average heat flux (Wieczorek and Phillips, 2000) in the PKT region could also have\nplayed a role in these magmatic activities. According to geodynamic models, the elevated heat flux,\ndriven by shallow enrichment in heat-producing elements, may have induced mantle upwelling,\nwhich in turn fueled the long-lived volcanism observed in the PKT region (Laneuville et al., 2018;\nLaneuville et al., 2013).”\nComments: Line 322, Deep feeder dykes below the PKT border are not shown in the schematic\nfig. 6.\nResponse: Thank you for your comments. Deep feeder dykes below the PKT border were added\nin the new Fig. 6.\nComments: Line 323-325. It is essential to discuss the implications of the proposed scenario on\nthe mineralogy and the surface age of the basalts in the region to substantiate the findings. I suggest\na brief discussion on the same may be provided.\nResponse: Thank you for your insightful suggestions. In response to your suggestion, we have\nincorporated a brief mineralogical analysis of both the Marius Hills region (Lines 393-415) and\nthe Lichtenberg crater area (Lines 423-432).\nComments: Line 333, “100nT” looks typo error; instead, it should be “100mGal”\nResponse:Thank you for your comment. We corrected the typo accordingly.\nComments: The paper is generally well-written, presenting some interesting and significant\nresults, but the novelty of the study needs to be brought out and highlighted more clearly and\nprominently.\nResponse: Thank you very much for your valuable feedback. We have revised the manuscript to\nmore clearly and prominently emphasize the unique contributions and significance of our findings\n(Lines 437-446). We hope these changes address your concern and enhance the overall impact of\nthe study.\nComments: Reviewer #3 (Remarks on code availability):\nWe were able to install and run the code and reproduce the results presented in the paper.\nResponse: Thank you very much for your efforts to run and reproduce the results."}, {"round": 2, "reviewer_comments": {"reviewer_2": "Review of “Unveiling Magmatic Structures and Connectivity beneath 1 Lunar Oceanus Procellarum Region: Insights from\nGRAIL Gravity Data” by Geng and colleagues.\nThis manuscript uses 3D density inversions to study the gravity anomalies in Oceanus Procellarum associated with the\nMarius Hills volcanic complex and the Aristarchus plateau. The results reveal a number of density anomalies consistent with\nlarge intrusions, consistent with previous studies. The results also reveal an interconnection of these density anomalies and\nconnection to the nearby border structures surrounding the Procellarum region. These results suggest an interconnected\nmagmatic system responsible for both the border anomalies and prominent volcanic regions.\nOverall, I find the authors have satisfactorily replied to my previous comments and revised the manuscript accordingly. The\nresults of this study are noteworthy and will be of interest to many in the planetary science community. The manuscript is\nwell written, the analyses are robust, and the conclusions are well supported. Below I have a number of minor comments for\nconsideration by the authors.\n[72] revise to “detailed modeling”\n[77] “based on data collected by” the spherical harmonic model of topography itself was not collected by LRO\n[FIgure 1] revise “Gmal” to “mGal”, also found in other figures\n[94-95] It may help some readers to indicate the “bridge” and “belt” in the figure (perhaps with arrows)\n[103] “topographical gradients” revise to “topographic gradients leading up to the adjacent highlands” to provide more\ncontext\n[123] “isosurface” It may be helpful to justify your use of isosurfaces here. Given the geological evidence for possible\nintrusive activity, the 3D continuous density inversions are thresholded at discrete cutoffs in order to delineate the extent of\npossible intrusive bodies.\n[127-128] “The bottom depth denotes…” It would be good to mention that bottom depths in these sorts of inversions are very\npoorly constrained because of the attenuation of gravity with distance (depth). I would caution against over-interpreting\nbottom depths. Other inversions of the PKT border anomalies had them at somewhat greater depths than this study (in which\nthe anomalies only extend to ~10 km depth). Either could be correct. Very small-scale anomalies must arise at very shallow\ndepths, but these features are broad enough that they could arise from deeper levels.\n[150] “The density contrasts… lie between” The data could likely also be fit by a discrete body with a single density, which\nmight be more likely for intrusions. The data doesn’t require such a wide range of densities within the intrusions. Some\ndiscussion of this point may be helpful to the reader.\n[202] It isn’t clear what this paragraph about these other regions (Flamsteed and Aristarchus) adds to this discussion. The\nresults of these inversions aren’t obviously inconsistent with previous interpretations of these anomalies as buried impact\ncraters, and no further interpretations are made. I would suggest moving most of this discussion to the supplement, and\nperhaps just stating in the main text that these features were also analyzed but the interpretation as magmatic structures was\ninconclusive. Similarly, much of the discussion of the sensitivity to density contrast could be moved to the supplement, with a\nsimple statement in the main text that different density contrasts have a modest effect on the thicknesses of the bodies but do\nnot change your conclusions or interpretations.\n[231] In the discussion of the relative merits of the buried impact crater vs intrusion interpretation, you might point out that the\napparent interconnection of these anomalies provides further support for the intrusive interpretation.\n[270-275] “likely to have formed from…” This proposed scenario is not unreasonable, but I do think you may be over-\ninterpreting these non-unique inversions. I suggest at least softening from “likely” to “may have formed”, and similarly\nelsewhere.\n[290] I believe you have the discussion of the seismics backwards. Should this be “increase in velocity below the upper 20-\n25 km of crust”?\n[302] “strong correlation with overlying wrinkle ridges” A similar correlation has been documented in Mare Frigoris (Williams\net al., https://doi.org/10.1016/j.icarus.2019.03.002)\n[347] revise “pressure weaknesses” to “structural weaknesses”\n[347-449] Although I appreciate the added discussion of the broader implications, much of this text is rather speculative. I\ndon’t have any particular disagreements, but it may be preferable to focus on the more robust conclusions in the main text\nand expand some of these more speculative interpretations in the supplement.\n[379] “magma plumbing… might not have achieved a state of complete cooling” In this paragraph, you discuss a history of 2-\n3 billion years of volcanism. An individual intrusion of the scale found in this study could not remain molten for billions of\nyears. As pointed out by the first reviewer, the cooling timescales are relatively easily calculated (there is an analytic solution\nto the solidification of a dike or sill presented in Geodynamics by Turcotte and Schubert). Rather, this long period would\nindicate reactivation of the system by multiple batches of magma. However, can all of these features be related to the\nmagma plumbing system? For instance, the lunar maria are found throughout the nearside mare region, not always in\nconnection with these gravity anomalies.", "reviewer_3": "The manuscript is well written, and the authors have adequately addressed the comments and incorporated all the\nsuggestions. Using GRAIL data, the manuscript presents a holistic view of the magmatic structures and interconnections\nunderneath the Oceanus Procellarum region of the Moon.\nThe study efficiently highlights the role of lateral magma transport in lunar magmatism.\nI strongly recommend that the manuscript be published."}, "author_response": "Response: Thank you for your positive feedback and thoughtful comments. We appreciate your recognition of our revisions and the value of our findings. We address the minor comments accordingly.\n[72] revise to “detailed modeling”\nResponse: It was revised accordingly.\n[77] “based on data collected by” the spherical harmonic model of topography itself was not collected by LRO\nResponse: It was revised accordingly.\n[FIgure 1] revise “Gmal” to “mGal”, also found in other figures\nResponse: We have revised “Gmal” to “mGal” in Figure 1 and other relevant figures as requested.\n[94-95] It may help some readers to indicate the “bridge” and “belt” in the figure (perhaps with arrows)\nResponse: Thank you for the suggestion. We decided not to add arrows for the “bridge” and “belt” in the figure. Since the belts generated by the border structures and the bridge connecting the Marius-SH and Marius-NH are clearly observable in Figure 2b. Additionally, some of these belts are intermittent, and adding such annotations could obscure detailed information. Therefore, we opted to keep the figure without these labels.\n[103] “topographical gradients” revise to “topographic gradients leading up to the adjacent highlands” to provide more context\nResponse: It was revised accordingly.\n[123] “isosurface” It may be helpful to justify your use of isosurfaces here. Given the geological evidence for possible intrusive activity, the 3D continuous density inversions are thresholded at discrete cutoffs in order to delineate the extent of possible intrusive bodies.\nResponse: Text below was added in the revised manuscript to explain this point.\n“Given the geological evidence for potential intrusive activity24, 25, 27, the 3D continuous density inversions are volumed-rendered using discrete cutoffs to delineate the extent of possible intrusive bodies. ”\n[127-128] “The bottom depth denotes…” It would be good to mention that bottom depths in these sorts of inversions are very poorly constrained because of the attenuation of gravity with distance (depth). I would caution against over-interpreting bottom depths. Other inversions of the PKT border anomalies had them at somewhat greater depths than this study (in which the anomalies only extend to ~10 km depth). Either could be correct. Very small-scale anomalies must arise at very shallow depths, but these features are broad enough that they could arise from deeper levels.\nResponse: We have discussed the uncertainties regarding the bottom depths of the inversion models and their causes in multiple sections of the manuscript. For instance, we note that the resolution of the gravity signal decreases with depth, and the density contrast between the anomalous body and surrounding rock diminishes at greater depths, as predicted by seismic models. These factors result in the bottom of the anomaly appearing shallower than its actual depth. Additionally, the impact of Bouguer correction density values, filtering, and inversion parameters on the inversion results, particularly on bottom depths, is addressed in the Methods section. Detailed analysis and images are provided in the supplementary materials for further clarification.\n[150] “The density contrasts… lie between” The data could likely also be fit by a discrete body with a single density, which might be more likely for intrusions. The data doesn’t require such a wide range of densities within the intrusions. Some discussion of this point may be helpful to the reader.\nResponse: Text below was added in the revised manuscript to explain this point.\n“We observe that the density contrast of the intrusions from the inversion gradually decreases from the center to the edges. This is mainly due to the smoothing nature of the inversion method used44; the actual density contrast variation within the intrusions may not be as pronounced. ”\n[202] It isn’t clear what this paragraph about these other regions (Flamsteed and Aristarchus) adds to this discussion. The results of these inversions aren’t obviously inconsistent with previous interpretations of these anomalies as buried impact craters, and no further interpretations are made. I would suggest moving most of this discussion to the supplement, and perhaps just stating in the main text that these features were also analyzed but the interpretation as magmatic structures was inconclusive. Similarly, much of the discussion of the sensitivity to density contrast could be moved to the supplement, with a simple statement in the main text that different density contrasts have a modest effect on the thicknesses of the bodies but do not change your conclusions or interpretations.\nResponse: Revised accordingly. We have moved the paragraph about Flamsteed and Aristarchus to Supplementary Note 1. The discussion of the sensitivity was moved to Supplementary Note 2. In the main text, we provide simple statements as suggested.\n[231] In the discussion of the relative merits of the buried impact crater vs intrusion interpretation, you might point out that the apparent interconnection of these anomalies provides further support for the intrusive interpretation.\nResponse: Thank you for the suggestion.It is a good point. We added this point to the manuscript as shown below.\n“Moreover, the observed interconnection of these anomalies, as discussed below, further supports their interpretation as intrusive bodies.”\n[270-275] “likely to have formed from…” This proposed scenario is not unreasonable, but I do think you may be over-interpreting these non-unique inversions. I suggest at least softening from “likely” to “may have formed”, and similarly elsewhere.\nResponse: Revised accordingly.\n[290] I believe you have the discussion of the seismics backwards. Should this be “increase in velocity below the upper 20-25 km of crust”?\nResponse: Our statement is correct. Here we discussed “velocity within the upper 20-25 km of the crust”. The text following this statement in our manuscript are as below, which discuss velocity within the upper 20-25 km of the crust and based on previous research: “Seismic data from the Apollo missions revealed a marked increase in velocity within the upper 20-25 km of the crust, correlating with porosity decrease due to increasing overburden pressure. For instance, it has been reported that lunar interior models M1 and M3 estimate subsurface densities of 2,680 kg/m³ at 10 km and 3,050 kg/m³ at 20 km for M1, and 2,760 kg/m³ at 10 km and 2,890 kg/m³ at 20 km for M3, respectively58. ”\n[302] “strong correlation with overlying wrinkle ridges” A similar correlation has been documented in Mare Frigoris (Williams et al., https://doi.org/10.1016/j.icarus.2019.03.002)\nResponse: Thank you for the information.This study area extends beyond our research area; therefore, it is not included in our paper.\n[347] revise “pressure weaknesses” to “structural weaknesses”\nResponse: Thank you for your suggestion. We reviewed the referenced literature and believe that \"pressure gradient\" is more accurate, so we have changed \"pressure weaknesses\" to \"pressure gradient.\" The text in our manuscript is as below: “Previous studies have suggested that systemwide pressure gradients can propagate quickly through the magmatic web, facilitating lateral magma transport in different volcanic systems, such as beneath Hawaiian70 and Kamchatka volcanoes78.”\n[347-449] Although I appreciate the added discussion of the broader implications, much of this text is rather speculative. I don’t have any particular disagreements, but it may be preferable to focus on the more robust conclusions in the main text and expand some of these more speculative interpretations in the supplement.\nResponse: Thank you for the suggestion. In our previous revision, based on Reviewer 3's comment that \"the arguments presented in the study need to be substantiated by the compositional and temporal relationships to build a complete picture of volcanism in the Oceanus Procellarum,\" we added an analysis of the mineralogical composition and temporal relationships among the formation of the Marius Hills volcanic complex, the western PKT border structures, and the wrinkle ridges. These discussions integrate our inversion model with existing studies on mineralogical and temporal relationships, further supporting the consistency of our results with other research findings. Additionally, our model provides more detailed support for related studies. Therefore, we believe it is reasonable to retain these discussions.\n[379] “magma plumbing… might not have achieved a state of complete cooling” In this paragraph, you discuss a history of 2-3 billion years of volcanism. An individual intrusion of the scale found in this study could not remain molten for billions of years. As pointed out by the first reviewer, the cooling timescales are relatively easily calculated (there is an analytic solution to the solidification of a dike or sill presented in Geodynamics by Turcotte and Schubert). Rather, this long period would indicate reactivation of the system by multiple batches of magma. However, can all of these features be related to the magma plumbing system? For instance, the lunar maria are found throughout the nearside mare region, not always in connection with these gravity anomalies.\nResponse: Thank you for your insightful comments. We agree with your point that \"this long period would indicate reactivation of the system by multiple batches of magma.\" We have added the following text to clarify this: “Consequently, the magma plumbing system beneath the PKT border is likely long-lasting and possibly involved reactivation of the system through injection of multiple batches of magma.”\nRegarding the second question, “can all of these features be related to the magma plumbing system?”, it is unlikely that all these features are related to the magma plumbing system. In regions like Flamsteed and Aristarchus, our model indicates no connection to the western PKT border intrusions. This study focuses on the Oceanus Procellarum region along the western border of the PKT. We believe that examining the magma plumbing system along the entire PKT border will provide deeper insights into its role in influencing volcanic activity within the PKT region."}]

Decision letter