1 March 19, 2019 Shallow seismic activity and young thrust faults on the Moon Thomas R. Watters1, Renee C. Weber2, Geoffrey C. Collins3, Ian J. Howley2, Nicholas C. Schmerr4 and Catherine L. Johnson5,6 1. Center for Earth and Planetary Studies, Smithsonian Institution, Washington, DC 20560, USA. 2. NASA Marshall Space Flight Center, Huntsville, AL 35805, USA. 3. Physics and Astronomy Department, Wheaton College, Norton, MA 02766, USA. 4. University of Maryland, Department of Geology, College Park, MD 20742, USA. 5. Dept. of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, British Columbia, V6T 1Z4, Canada. 6. Planetary Science Institute, Tucson, AZ 85719, USA. Correspondence and requests for materials should be addressed to 1T.R. Watters (e-mail: watterst@si.edu). 2 March 19, 2019 Supplementary Information Supplementary Text, Supplementary Figure S1-S6, Supplementary Tables S1-S3, and Supplementary References. LOCSMITH Seismic Event Location The LOCSMITH event location algorithm1 was developed for sparse seismic networks. The scheme uses an iterative, adaptive grid search and accounts for uncertainty in arrival time by means of windows around the true arrival time. The grid point (possible quake location) is refined for the next iteration when the gridded theoretical arrival time is found to lie within the arrival time widows for all stations and all seismic phases. If the theoretical arrival time is not found, the grid point is considered falsified. The algorithm thus produces a cloud of non-falsified grid points, points non-falsified by any arrival time or back azimuth data, with acceptable solutions. The spatial distribution of points in a relocation cloud is controlled by the number of observed arrivals, with the best clouds having a ball shape similar to a typical quake location error ellipse2,3. The LOCSMITH method only requires arrival uncertainties. The reference arrival uncertainties assigned were 3 seconds for the P wave and 10 seconds for S wave, following the method used for deep moonquakes by Hempel et al. (3). LOCSMITH returned acceptable solutions ranging in depth from the surface to 300 km beyond which the search was terminated. However, only relocated epicenters with surface locations were evaluated. Although surface solutions have been obtained in previous studies, the shallow moonquakes have generally been attributed to deeper crustal sources (see 4). Based on the study of moonquake coda (the total length of the seismic wavetrain), Gillet et al. (4) found that the energy decay curve differs significantly when the source is above or below the depth of 3 March 19, 2019 the high attenuation and high scattering megaregolith. In the study of coda of impacts and shallow moonquakes, Gillet et al. (4) find a persistent separation between shallow moonquakes and impacts (their Fig. 7), whereas Blanchette-Guertin et al. (5) using a different metric find very similar distributions of decay times for shallow moonquakes and impacts (their Fig. 5). It is important to note that both studies find overlap in the populations of surface impact events and shallow moonquakes. These studies suggest that the source mechanism and how the energy is partitioned into the wavefield does affect the resulting coda and that source depth is not the only source of differences in the coda. Moreover, the overlap in decay times of the populations of surface impact events and shallow moonquakes supports near-surface sources for the shallow moonquakes. Seismic Shaking Map In an effort to constrain the distance over which ground motion from a shallow seismic event is significant, seismic shaking proximal to an active scarp is estimated using an open- source 3-D finite-difference code, Serpentine Wave Propagation6 (WPP), adapted for the Moon to produce time-acceleration histories (i.e., synthetic seismograms) for individual rupture events. WPP solves the equation of motion on a 3-D finite difference, non-uniform Cartesian mesh, and is designed to run on parallel multi-CPU clusters. The user may specify elastic structures in 3-D, including topography, layering, lateral heterogeneity, attenuation, and multiple sources with varying source time functions. The code is primarily used for regional simulations of wave propagation in complex media. The elastic response of the medium is calculated for dominant frequencies up to 10 Hz or greater, though computational resources required grow for higher frequencies. The source is implemented as a double-couple earthquake, with fault parameters 4 March 19, 2019 obtained from the elastic dislocation model software COULOMB available from USGS Earthquake Hazards Program. For our simulations, we produced seismic wavefields with WPP up to dominant frequencies < 10 Hz. Adaption of the code for the Moon involves specifying the subsurface structure using results from the Apollo Passive Seismic Experiment, including minimal attenuation (i.e., a high Q value) (7), and reasonable velocity values for the 1D layered structure (e.g., regolith shear wave velocity vp of 200 m/s) from Toksoz et al. (8). We produced synthetic seismograms in a station grid surrounding the event, allowing us to the map the predicted peak ground acceleration proximal to an active scarp (Fig. 4). The WPP modeling incorporates the surface topography of the Mandel’shtam region obtained from LOLA9 on LRO and implemented in the non-Cartesian mesh portion at the top of the WPP grid. The subsurface model used is modified from Toksoz el al., (8), and is presented in Table S1. No attenuation structure is specified in the model, but random seismic scattering in the uppermost 1 km of lunar megaregolith is simulated by adding 25% root mean square heterogeneity in Vs and Vp with a correlation scale length of 50 m to the uppermost 1 km of the seismic model10. Displacements associated with fault slip events using COULOMB11, 12 are constrained by topographic profiles across the fault scarp. Elastic dislocation modeling allows the fault geometry and magnitude and sense of slip to be specified and the stresses and material displacements in the vicinity of the fault are then completely determined by using the stress functions for an elastic half-space13 using boundary element methods14. Assumed model parameters for near-surface lunar crust are a coefficient of static friction µs of 0.85, Poisson’s ratio ν of 0.25, and values for the elastic modulus E of 40 GPa, appropriate for a highly fractured, weakened upper lithosphere (e.g., 15, 16). We have chosen a scarp in the Mandel’shtam cluster (6.9°N, 161°E) as representative of lunar lobate thrust fault scarps17. The topographic profile 5 March 19, 2019 obtained from a digital elevation model derived from LROC NAC stereo images across the northern segment (Fig. S1, S2) indicates the fault scarp has ~70 m of relief at this location. Figure S1. Digital elevation model of the Mandel’shtam scarp generated using LROC NAC stereo images (frames M191895630 and M191909925). The DEM has a horizontal spatial scale of 5 m/pixel (NAC stereo images have a resolution of ~1 m/pixel) and a vertical precision of ~0.5 m. Elevations are referenced to a sphere of 1,737,400 m. Note that the line is the location of the topographic profile used to constrain the COULOMB modeling. The profile is shown in Fig. S2. 6 March 19, 2019 The best fit model parameters for Mandel’shtam are a fault dip of 35°, maximum fault depth of 700 m, and a cumulative slip of 110 m (Fig. S2). These are similar to model parameter values for five other lunar scarps that includes another scarp in the Mandel’shtam cluster16. Coulomb calculates the total seismic moment for the source from the fault area, slip, Young’s modulus, and Poisson’s ratio. The total length of the scarp is ~12.3 km. The length of the fault segment characterized by the topographic profile is ~1.98 km. The corresponding total seismic moment = 3.89*1018 Nm and a Mw = 6.36. Figure S2. Mandel’shtam scarp topography and elastic dislocation model of displacement. The topographic profile (black) was extracted from an LROC NAC stereo derived digital elevation model (DTM) with horizontal spatial scale of 5 m/pixel. The best-fit model displacement (red) is for a linear fault with a dip of 35°, a maximum fault depth of 700 m, and a cumulative slip of 110 m. The location of the topographic profile is shown in Fig. S1. It is expected that the greatest surface disruption proximal to a scarp will occur during the largest event the fault experienced (as shown for Mandel’shtam in Fig. 4). The shakemaps, based on a modeled fault plane determined from the fault parameterization values obtained from elastic 7 March 19, 2019 dislocation modeling (strike=0º, dip=35º, and rake=90º), give a sense of the area effected by “significant” seismic shaking, and its distance from the Mandel’shtam scarp. The peak ground surface motion is measured in every 1x1 km block in the model at a reference frequency of 1 Hz. The spread in points in the plots is due to the effect of topography and radiation pattern of the fault. In an effort to objectify what is meant by “significant” seismic shaking we adopt the intensity scale based on Worden et al. (18). On this scale, shaking intensity is classified as follows: Weak: 0.3% g; -51 dB Light: 2.8% g; -31 dB Moderate: 6.2% g; -24 dB Strong: 12% g; -18 dB Very Strong: 22% g; -13 dB Severe: 40% g; -8 dB Violent: 75% g; -3 dB Extreme: 139% g; +3 dB Where g is the acceleration due to lunar gravity (1.62 m/s2) and 100% g is equal to 0 dB. The shakemap for vertical ground motion (Fig. 4A, B) indicates strong to moderate shaking out to a distance of ~12-30 km and the shakemap for horizontal motion (Fig. 4C, D) indicates strong to moderate shaking out to ~7-15 km. The magnitude obtained for Mandel’shtam (Mw = 6.36) from elastic dislocation modeling is a cumulative maximum event magnitude. A series of smaller magnitude events is expected. Table S2 shows the peak vertical and horizontal ground shaking that would occur at the source for varying moment magnitudes of shallow moonquakes. Similar to the shake maps (Fig. 4), these are reported relative to lunar gravity, where a 40 dB change in shaking corresponds to a 100 times change in acceleration. The peak is evaluated by finding the maximum ground acceleration occurring within the simulation for each component of motion. 8 March 19, 2019 Distances between Mapped Lobate Scarps and Relocated Epicenters The spatial relation of relocated epicenters and lobate scarps were determined by measuring the distance between mapped scarps and epicenters with surface solutions in the relocation clouds. Relocated epicenters with the minimum distance to a nearby scarp were selected. Distances were measured between the relocated epicentral locations and points near the center of each mapped scarp. Results show the minimum distance between 8 of the 17 relocated epicenters and mapped lobate scarps is less than 30 km (Fig. S3A). In an effort to determine if the proximity of the 8 relocated epicenters is simply a random effect, 10,000 sets of 17 epicentral locations with a globally uniform random spatial distribution were generated. The analysis indicates that less than 4% of these 170,000 random quake locations are within 30 km of a fault scarp (Fig. S3B) and that none of the 10,000 sets have more than five random quake locations within 30 km of a scarp (Fig. S3C). Extending the minimum distance by a factor of 2 results in only three of the 10,000 sets with eight or more random quake locations within 60 km of a scarp (Fig. S3D). 9 March 19, 2019 Figure S3. Histograms of distribution of the distances between relocated and random epicenters and lobate scarps. A) Distribution of minimum distances between 17 relocated epicentral locations and mapped fault scarps. B) Distribution of the minimum distance between 170,000 randomly generated epicenters and mapped fault scarps. C) Distribution of the minimum distance of random epicentral locations within 30 km of mapped nearside scarps in 10,000 sets of 17 randomly generated epicentral locations. D) Distribution of the minimum distance of random epicentral locations within 60 km of mapped nearside scarps in 10,000 sets of 17 randomly generated epicentral locations. Demodulating the Time Series The time series is demodulated by first removing its mean, computing the envelope function and then dividing the de-meaned series by the difference between the upper and lower 10 March 19, 2019 envelopes evaluated at the times of the input series values. The resulting plot is a time series with a flat envelope consisting of values between ±0.5 (see Fig. S4). Demodulation removes long-period variation from a time series, but does not remove the effect of the non-uniform lunar orbital motion. This is illustrated in Fig. S4A which shows increased bin values in the first and last bins (times near perigee and apogee) - this is the same as the distribution of a sine wave. The background EMD distribution is made uniform by multiplying it by 2 and taking the inverse sine (Fig. S4B). This ensures that the histogram computed at moonquake times is not biased by the uneven sampling of the Moon’s orbital position with respect to the Earth. The histogram of the demodulated, uniform EMD at shallow quake times is shown in Fig. S4D. The ratio of the number of occurrences of a particular EMD at quake times to the number of occurrences of this EMD over the Apollo experiment duration, and the ratio of the demodulated EMD at quake times to the background demodulated EMD are shown in Fig. S4E and S4F respectively. 11 March 19, 2019 Figure S4. Histograms of the EMD and demodulated EMD at the times of the shallow moonquakes. The background EMD (S4A, B) is a binning of the value of the EMD as discretized in 1-day time steps. The left column figures (S4A, C) are before demodulation and right column figures (S4B, D) are after demodulation. Demodulating removes all but the monthly periodicity (anomalistic month, or time from perigee to perigee), accounting for the bias introduced by the fact that the Moon spends more time near apogee than near perigee. Stresses due to True Polar Wander True polar wander (TPW) occurs when a body reorients itself relative to its axis of rotation driven by mass redistribution19. It involves motion of the surface through the body’s 12 March 19, 2019 polar axis19, 20. The reorientation direction is generally roughly perpendicular to the tidal axis. True polar wander induces extensional stresses in the region leading the reorientation direction because the spin axis there must lengthen and compressional stresses in the trailing region because the spin axis there must shorten19, 21, 22. Equations for the stresses can be found in Matsuyama and Nimmo (19) and Kean and Matsuyama (23). Comparison of Predicted and Observed Fault Orientations The net non-isotropic compressional stresses sn from sr+sw+st result in thrust faults with preferred orientations24. Although the contribution of diurnal tidal stresses to sn is small (≤5 kPa), the addition of st results in peak stresses generally occurring when the Moon is near apogee or perigee. The magnitude of the orbital recession stresses sr is very sensitive to h2, the Love number defining the vertical displacement in response to a gravitational potential, and the change in orbital radius a - ao determined by the period of recession. The value of h2 used to estimate a maximum sr is the long-term (time scale more than ~104 yr) h2 @ 1.8. This value assumes the presence of a 25 km thick elastic lithosphere and is a more appropriate value for tectonic time scales (I. Matsuyama, Personal Communication; see 25). The maximum stress for a period of recession of 1.1 Byr (likely duration of the Copernican period) is ~220 kPa and assumes an average rate of recession over the Copernican of 3 cm/yr (see equations in 19). The short-term, tidal forcing period value of h2 = 0.0371 (26) was used to estimate st (see 24). Consequently, the stress inferred from the population of young lobate thrust fault scarps and modeling indicates a present-day low level of compressional stress in the lunar lithosphere. The orientation of each mapped fault segment was compared to the local orientation of the principal axis of most compressive stress calculated from general models for orbital 13 March 19, 2019 recession, despinning, orbital recession and diurnal tides, and orbital recession and true polar wander (Fig. S5). The expected faults should be close to perpendicular to the local axis of most compressive stress if the tidal and polar wander stresses are influencing orientations of the fault population21. The predicted fault orientations from the four models are generally in good agreement with the observed orientations of the lobate scarps, and are distinctly different from a simulated population of faults with random orientations from isotropic contraction alone. Figure S5. Comparison of observed and expected local thrust fault orientations for different stress models. An angle of 0° indicates fault segments that are orthogonal to local direction of most compressive stress, and 90° indicates fault segments that are parallel to most compressive stress. Analysis of stresses around times of the shallow moonquakes To determine which of the near apogee shallow moonquakes occur when compressional stresses are dominant, the stresses from orbital recession plua diurnal tidal stresses were evaluated at the 14 March 19, 2019 coordinates of the relocated epicenters for 15 days on either side of the recorded date of the shallow moonquake. Diurnal stresses were calculated using the lunar ephemeris from 1971 to 1976, which was downloaded from the Jet Propulsion Laboratory (JPL) Horizons website. Plots of compressional stress around the moonquake dates (see examples in Fig. S6) show that compressional stress remains high for five to seven days on either side of the peak stress. Of the set of 12 relocated epicenters that occur near apogee, three events occur at or within one day of local peak compression, another five occur within seven days of local peak compression, and the remaining four events occur more than seven days from local peak compression (Table S3). Of the 8 events that occur near apogee and are within 30 km of a fault scarp, 5 occur within 7 days of peak stresses (Table S3) and the other 2 events greater than 30 km but within 60 km occur less than 7 days before or after peak stresses. This does not exclude the possibility that events more than seven days from local peak compression are not triggered by tidal stresses. Interestingly, the N2 event, which occurs outside of near-apogee, also occurs when the most compressive stress is near the minimum for its location. Figure S6. Plots showing stress accumulation due to recession plus diurnal stresses over the period example shallow moonquakes occur near apogee. Isotropic compressional stresses due to 15 March 19, 2019 global contraction are not included. The addition of the stress component due to global contraction will increase the magnitude of the total stress by a constant value (i.e., 2 MPa). The plots show the increase and decrease of compressional stress (tension is negative) at the latitude and longitude coordinates for the relocated epicenters of the shallow moonquakes A) N14(1), and B) N25. The N14(1) event occurs when compressional stress peak. The N25 event occurs before peak compression is reached. Day 0 is the day of the event. Note that the total local most and least compressive stresses will never become tensile because isotropic compressional stresses due to global contract far exceed the maximum tidal stresses. References 1. Knapmeyer, M. Location of seismic events using inaccurate data from very sparse networks. Geophys. J. Int., 175, 975-991, (2008). 2. Weber, R. C., Knapmeyer, M., Panning, M. & Schmerr, N. Modeling approaches in planetary seismology. in Extraterrestrial Seismology, V. C. H. Tong, R. A. Garcia, Eds. (Cambridge University Press, New York, NY, 2015), pp. 140-158. 3. Hempel, S., Knapmeyer, M., Jonkers, A. R. T., Orberst, J. Uncertainty of Apollo deep moonquake locations and implications for future network designs, Icarus, 220, 971-980, (2012). 4. Gillet, K., Magerin, L., Calvet, M. & Monnereau, M. Scattering attenuation profile of the Moon: Implications for shallow moonquakes and the structure of the megaregolith. Phys. Earth and Planet. Inter., 262, 28-40 (2017). 5. Blanchette-Guertin, J. -F., Johnson C. L. & Lawrence, J. F. Investigation of scattering in lunar seismic coda. J. Geophys. Res., 117, E06003, doi:10.1029/2011JE004042 (2012). 6. Petersson, N. A., & Sjögreen, B. Reference guide to WPP version 2.0, LLNL Tech. Rep., LLNL-TR-422928 (2010). 7. Dainty, A. M. & Toksoz, M. N. Elastic wave propagation in a highly scattering medium – A diffusion approach. J. Geophys. Res., 43, 375-388 (1977). 8. Toksoz, M. N. et al. Seismic scattering and shallow structure of the moon in Oceanus 16 March 19, 2019 Procellarum, Moon, 9, 11-29 (1974). 9. Smith, D. E. et al. Initial observations from the Lunar Orbiter Laser Altimeter (LOLA). Geophs. Res. Lett. 37, L18204, doi:10.1029/2010GL043751 (2010). 10. Frankel, A. & Clayton, R. W. Finite-difference simulations of seismic scattering - Implications for the propagation of short-period seismic-waves in the crust and models of crustal heterogeneity, J. Geophys. Res., 91, 6465-6489, doi:10.1029/JB091iB06p06465, 1986. 11. Toda, S. et al. Forecasting the evolution of seismicity in southern California: Animations built on earthquake stress transfer. J. Geophys. Res., 110, B05S16, doi:10.1029/2004JB003415 (2005). 12. Lin, J. & Stein, R. S. Stress triggering in thrust and subduction earthquakes, and stress interaction between the southern San Andreas and nearby thrust and strike-slip faults. J. Geophys. Res., 109, B02303, doi:10.1029/2003JB002607 (2005). 13. Okada, Y. Internal deformation due to shear and tensile faults in a half-space. Bull. Seism. Soc. Am., 82, 1018-1040 (1992). 14. Crouch, S. L. & Starfield, A. M. Boundary Element Methods in Solid Mechanics, 322 pp., Allen and Unwin, Winchester, Mass., 1983. 15. Watters, T.R., Thomas, P.C. & Robinson, M.S. Thrust faults and the near-surface strength of Asteroid 433 Eros, Geophys. Res. Lett., 38, L02202, doi:10.1029/2010GL045302 (2011). 16. Williams, N. R. et al. Fault Dislocation Modeled Structure of Lobate Scarps from Lunar Reconnaissance Orbiter Camera Digital Terrain Models. J. Geophys. Res., doi: 10.1002/jgre.20051 (2013). 17. Watters, T. R. et al. Evidence of recent thrust faulting on the Moon revealed by the Lunar Reconnaissance Orbiter Camera. Science, 329, 936-940 (2010). 18. Worden, C. B., Gerstenberger, M. C., Rhoades, D. A. & Wald, D. J. Probabilistic relationships between ground-motion parameters and modified Mercalli intensity in California. Bull. Seism. Soc. Am., 102, 204-221 (2012). 19. Matsuyama, I. & Nimmo, F., Tectonic patterns on reoriented and despun planetary bodies, Icarus, 195, 459-473 (2008). 20. Melosh, H. J. Tectonic patterns on a reoriented planet—Mars. Icarus, 44, 745–751 (1980). 17 March 19, 2019 21. Collins, G. C. et al. Tectonics of the outer planet satellites. in Planetary Tectonics, T. R. Watters, R. A. Schultz, Eds. (Cambridge Univ. Press, New York, NY, 2010), pp. 264– 350. 22. Leith, A. C., & McKinnon, W. B. Is there evidence for polar wander on Europa? Icarus, 120, 387–398 (1996). 23. Keane, J. T., & Matsuyama, I. Evidence of lunar true polar wander, and a past low- eccentricity, synchronous lunar orbit. Geophys. Res. Lett., 41, 6610-6619 (2014). 24. Watters, T. R. et al. Global thrust faulting on the Moon and the influence of tidal stresses. Geology, 43, 851–854 (2015). 25. Matsuyama, I. et al. GRAIL, LLR, and LOLA constraints on the interior structure of the Moon, Geophys. Res. Lett., 43, 8365-8375, http://doi.org/10.1002/2016GL069952, (2016). 26. Mazarico, E. M. et al. Detection of the lunar body tide by the Lunar Orbiter Laser Altimeter: Geophys. Res. Lett., 41, 2282–2288, doi:10.1002/2013GL059085 (2014). 27. Nakamura, Y. et al. Shallow moonquakes: Depth, distribution and implications as to the present state of the lunar interior. Proc. Lunar Sci. Conf. 10th, 2299-2309 (1979). 28. Oberst, J. Unusually high stress drops associated with shallow moonquakes. J. Geophys. Res., 92(B2), 1397–1405, doi:10.1029/JB092iB02p01397 (1987). Table S1. 1-D Background seismic velocity model Depth (km) Vp (m/s) Vs (m/s) Density (kg/m3) 0 500 285 1200 1 1000 560 2500 15 3000 1715 2800 38 5000 2850 3300 This subsurface velocity model is modified from Toksoz el al. (6). 18 March 19, 2019 Table S2. Peak Ground Motion for Thrust Fault Sources Mw Peak Vertical Peak Horizontal (ms-2) (dB @ 1.62 ms-2) (ms-2) (dB @ 1.62 ms2) 6.36 3.7 x100 7.2 2.9 x100 5.1 6.00 1.1 x100 -3.4 8.4 x10-1 -5.7 5.00 3.4 x10-2 -33 2.7 x10-2 -35 4.00 1.1 x10-3 -63 8.4 x10-4 -65 3.00 3.4 x10-5 -94 2.7 x10-5 -95 The peak ground motion is determined by finding the maximum ground acceleration occurring within the simulation for each component of motion. Table S3. Relocated Epicenters that Occur Near Apogee N# lon lat Distance Timing of Event with respect to time of local diurnal peak compression (including recession stress) Stress Drop (MPa)* N6 42.96 37.28 57.81 2 days after peak 2.5 N7 -74.42 -12.17 15.05 11 days after peak 1 N10 -74.38 -8.07 129.73 13 days after peak 30 N11 -33.57 -16.86 82.44 4 day after peak 2 N14 (1) 30.99 20.09 13.23 At peak compression 0.5 N14 (2) 25.56 -42.89 16.37 1 day before peak 0.5 N16 82.23 19.39 4.23 At peak compression 55 N20 -23.49 -15.96 131.34 11 days after peak 4 N21 -29.81 -51.77 57.76 5 days after peak 4 N23 59.25 8.42 62.98 14 days before peak 10 N25 39.8 33.58 28.59 5 days before peak 15 N26 -17.18 47.75 23.62 7 days before peak 40 Distance is the minimum distance between the relocated epicenter and a mapped lobate scarp. The N# is event number of the 28 recorded shallow moonquakes27. Hypocenters within 30 km of a fault scarp and the quake occurred within 7 days of peak stresses are shown in red and those beyond 30 km but within 60 km are shown in orange. Events that occur more than 11 days before or after peak stresses are shown in green. One events within 7 days of peak stresses, but beyond 60 km from a fault scarp is shown in black. *Stress drop estimates are from Oberst (28).