Velocity Measurements of Powdered Rock at Low Confining Pressures and Comparison to Lunar Shallow Seismic Velocity

Seismic methods will be useful for future lunar near‐surface characterization, and high‐fidelity elastic models will be required to aid interpretation of seismic observations. To develop an elastic lunar near‐surface model, we performed ultrasonic velocity measurements of lunar regolith simulant at low confining pressure and developed a rock physics model calibrated to these measurements. Grain contact models based on Hertz‐Mindlin theory produce accurate results at high confining pressure (i.e., several hundred meters or more burial depth) but historically fail to predict observed velocities in unconsolidated media at low pressure. Therefore, we heuristically modified existing models to fit our measured data over a range of porosities and confining pressures. To compare with Apollo 14 and 16 active seismic experiments, we used our new heuristic rock physics model to produce lunar subsurface velocity profiles. We performed ray tracing through our velocity profiles to calculate seismic traveltime, which results in good agreement with first arrivals interpreted from the Apollo experiments. Our model suggests a slightly higher velocity‐pressure dependence than inferred from in situ measurements, which may be due to porosity reduction in the lunar regolith from impact‐induced and natural vibrations.


Introduction
Active source seismic experiments have been used to investigate Earth's subsurface since the early 1900s (Dragoset, 2005;Hübscher & Gohl, 2014), and much of what is known about how the Earth has structurally evolved over time has come from seismic studies.Similarly, seismic studies using controlled active sources as well as natural passive seismic sources (e.g., earthquakes and meteorite impacts) have revealed details of the interior structures of the Moon and Mars (e.g., Cooper et al., 1974;Lognonné et al., 2020;Toksőz et al., 1972).While historically the seismic studies of extraterrestrial bodies have focused on body-wide or large features (kilometer depth scale), less work has focused on the near-surface seismic properties of these bodies (meter depth scale).Near-surface seismic studies are regularly employed on Earth for geohazard detection, site characterization, and resource exploration (e.g., Abd El-Aal & Mohamed, 2010;Gálfi & Pálos, 1970;Hack, 2000;McCann & Forster, 1990;Pegah & Liu, 2016;Steeples, 1990) and may be a viable technique for similar studies on extraterrestrial bodies.With the increase in lunar exploration activity associated with the Artemis program, which may place humans on the Moon as early as 2026, increasing our understanding of the Moon's near-surface in terms of geotechnical safety and resource potential will be paramount for sustainable space development (Metzger, 2016;Spudis & Lavoie, 2011).
Seismic velocity is the speed at which a mechanical seismic wave travels through a material, and is controlled by the elastic properties (i.e., compressive stiffness, shear stiffness, and bulk density) of the material.A highly dense and incompressible material will display high seismic velocity, and a compressible porous material will display low seismic velocity.Measuring seismic velocity and its changes with depth provides information to help identify subsurface materials and structures.For example, a sharp change in seismic velocity observed from the Apollo 14 active seismic data is interpreted to occur at a depth of ∼8.5 m where unconsolidated regolith overlies a much more rigid rock layer at the experiment location (e.g., Kovach & Watkins, 1973).To offer plausible interpretations of observed seismic velocities, we can create elastic models of various subsurface scenarios to see if these models can reproduce our velocity measurements.To create these models for the lunar near-surface, we must understand the relationships between the elastic properties of lunar regolith and properties such as porosity and confining pressure with depth.The results we present provide calibration for these relationships to enable lunar near-surface seismic modeling.
This study focuses on laboratory ultrasonic velocity analysis of a lunar regolith simulant.Our goal is to develop a predictive rock physics model that is calibrated to laboratory velocity measurements and is consistent with seismic velocity observed in the shallow (<10 m) lunar regolith from Apollo seismic experiments.We achieve this by (a) conducting compressional and shear velocity measurements of lunar regolith simulant under variable confining pressure at constant low temperature, (b) modifying existing rock physics models to predict our experimental velocities, and (c) forward-modeling seismic traveltimes to compare with the Apollo 14 and 16 active seismic experiments.Our data and rock physics model may be used in future studies for predicting the nearsurface seismic response of various geotechnical and resource scenarios to evaluate the potential effectiveness of lunar seismic exploration campaigns.

Background
Seismic exploration of the Moon was attempted with the early Ranger missions; however, all of these attempts were unsuccessful due to various component failures (Lognonné & Pike, 2015).The first successful lunar seismometer was deployed as part of the Early Apollo Surface Experiment Package during the Apollo 11 mission and was followed by a passive seismic experiment as well as multiple active seismic experiments (Nakamura, 2015;Nunn et al., 2020).For the active source experiments, astronauts deployed various types and sizes of explosives over a limited study area to create seismic waves.The seismic data collected during these studies were crucial to develop an early understanding of the structure of the lunar crust (Cooper et al., 1974;Gangi & Yen, 1979;Horvath et al., 1980;Kovach & Watkins, 1973;Latham et al., 1970;Mark & Sutton, 1975;Nakamura et al., 1982;Toksőz et al., 1972;Watkins & Kovach, 1973).Modern technology and techniques are allowing for further analysis of these seismic data over 40 years later (Dal Moro, 2015;Sollberger et al., 2016).
Collating existing analyses of various lunar seismic experiments (Table 1) yields a near-surface (<10 m) compressional seismic velocity (V P ) of roughly 100-140 m/s and shear velocity (V S ) of 35-100 m/s.This range of interpreted velocities reflects the difficulty in identifying the first arrival time of seismic waves in the recorded Apollo data due to a mixture of experimental uncertainty and data quality (Brzostowski & Brzostowski, 2009) and is lower than velocities typically observed in near-surface terrestrial settings (e.g., Bachrach et al., 1998;Shen et al., 2016).Seismic velocities depend on materials' elastic moduli and density, and are given by where K eff and μ eff are the effective bulk and shear moduli and ρ is the bulk density.As the moduli (i.e., stiffness) of the material increase, the velocity increases.The stiffness of a granular material is controlled by porosity, confining pressure (i.e., grain contact force), sorting, grain size, and grain shape (Prasad & Meissner, 1992;Zimmer et al., 2007aZimmer et al., , 2007b)), whereas grain mineralogy appears to have a negligible effect (Blangy et al., 1993).When comparing seismic velocity observations from lunar regolith to terrestrial examples, we must also consider the reduced lunar gravity and lack of moisture in lunar conditions.Although reduced gravity would weaken grain contacts and lower velocity, the lack of moisture may result in higher elastic moduli and velocity (Clark et al., 1980;Murphy et al., 1984).
The low near-surface seismic velocities interpreted from Apollo data are caused by a layer of fine-grained regolith.This ubiquitous layer of pulverized, angular, and glassy material with a global depth up to ∼20 m was formed by innumerable impact events over billions of years (McKay et al., 1991).Samples of the lunar regolith as well as intact rocks were returned to Earth for detailed analysis including ultrasonic velocity measurements (Anderson et al., 1970;Chung, 1972;Johnson et al., 1982;Kanamori et al., 1970Kanamori et al., , 1971;;Mizutani & Newbigging, 1973;Talwani et al., 1974;Todd et al., 1972), although only a few of these studies report measured velocities at low confining pressures (Table 1).Talwani et al. (1974) measured velocities at 0.1-250 MPa confining pressure (equal to 25-1,400 m depth), with their low-pressure measurements yielding V P between 170 and 250 m/s from Apollo 17 samples.Johnson et al. (1982) measured sample velocity at several pressure steps between 0.004 and 0.2 MPa, representing the most detailed dataset of lunar regolith velocity at near-surface confining pressures and yielding low-pressure V P measurements of 125-220 m/s from Apollo 15 regolith.The lowest pressure measurements reported by Johnson et al. correspond to a lunar depth >5 m, and their lowest reported velocity of 125 m/s is the only velocity measurement from all studies which explicitly agrees with the low (100-140 m/s) velocity interpreted from the Apollo active seismic experiments.Increasing our comprehension of the elastic properties of the lunar near-surface requires a better understanding of the elastic behavior of granular materials at low confining pressure, a topic which is relatively underexplored compared to higherpressure settings (Prasad et al., 2005).
Seismic velocity in granular media is sensitive to changes in confining pressure and is often expressed by a pressure exponent.The pressure exponent describes the shape of the polynomial relationship between seismic velocity and confining pressure.Table 1 illustrates pressure exponents published for lunar and terrestrial granular media.Laboratory measurements by Hardin and Richart (1963), Johnson et al. (1982), andZimmer et al. (2007a) showed pressure exponents between one-fifth and one-third.These observations do not agree with several existing contact radius models (CRM) (e.g., Mindlin, 1949;Walton, 1987), which assume grains are equal-sized spheres and predict that velocity should vary with effective pressure to the one-sixth power.Duffy and Mindlin (1957) experimentally showed that the velocity-pressure dependence of metal spheres transitioned from onefourth at low pressure to one-sixth at high pressure, suggesting that idealized contact models may not accurately describe the low-pressure behavior of granular materials.Additionally, Goddard (1990) described how altering the idealized contact geometry (i.e., a cone contacting a plane rather than two spheres in contact) can produce a one-fourth power pressure dependence, suggesting that complex contact geometries resulting from angular grains may contribute to the pressure dependencies observed from experimental data.These observations point to the inadequacy of idealized contact models for accurately predicting how velocity varies with pressure from increasing subsurface depth, especially at low-pressure near-surface conditions.
An important consequence of contact models not accurately representing granular media at low effective pressure is an overprediction of velocity and an underprediction of V P /V S ratio.Prasad (2002) reported anomalously high V P /V S ratios in sediments at low pressure, values which are hard to reconcile with contact models.Makse et al. (1999Makse et al. ( , 2004) ) discussed how the shear modulus of high-porosity granular systems violates assumptions inherent in many contact models, which could contribute to the discrepancy between modeled and observed velocities.While some existing contact models may be parameterized to provide an acceptable match to the observed V P in some studies, a recurring observation is an overprediction of V S .To account for these differences, several authors have proposed models and model modifications that reduce the tangential stiffness of the grain pack to better predict experimental observations.Bachrach et al. (2000) proposed to extend standard Hertz-Mindlin contact theory (Mindlin, 1949) with two additional parameters that allow for fractional slip at grain contacts and correct for grain contact area of non-spherical grains.The reduced contact area due to grain angularity is illustrated in Bachrach et al. (2000, their Figure 2).This modeling approach has recently been used by Caicedo et al. (2023) to describe the elastic behavior of Martian regolith simulant at low pressure.Zimmer et al. (2007a) used an empirical fit to laboratory data that relies on a void-ratio relation and corrects for compaction.Dutta et al. (2010) focused on empirical relations between coordination number (the average number of contacts per grain), porosity, and pressure to correct for observed V P /V S ratios, although one shortcoming is the need for differing coordination numbers in the calculation of V P and V S .It should also be noted that these rock physics modeling studies focused on terrestrial sands and glass beads, which may not share the same geotechnical properties as the angular and fine-grained lunar regolith.However, these studies provide a starting point to model the elastic behavior of granular material at low effective pressures.

Lunar Regolith Simulant
Colorado School of Mines Lunar Highlands Type 1 (CSM-LHT-1) regolith simulant was designed and manufactured on campus at Colorado School of Mines.CSM-LHT-1G, used in this study, is a variant which includes synthetic anorthosite-composition glass fragments.This simulant is designed to mimic the elemental and geotechnical properties of the lunar highlands regolith, which was encountered at the Apollo 16 landing site.
Particle size analysis was performed for CSM-LHT-1G and two Apollo 16 regolith samples, as shown in Figure 1a.CSM-LHT-1G and sample 64501 display very similar grain size distributions, with sample 67461 displaying a slightly coarser distribution.Overall distribution shapes for Apollo 16 samples and CSM-LHT-1G are similar, with mode grain sizes between 350 and 590 μm.In a following section, we compare our collected velocity data to laboratory measurements by Johnson et al. (1982) from sample 15301 of Apollo 15.The grain size distribution for sample 15301 reported by Graf (1993) is presented in Figure 1b along with our CSM-LHT-1G data binned to match the reported sieve sizes, showing an overall finer distribution from sample 15301.
Figure 2 shows the variety of grain sizes and angularities of CSM-LHT-1G.
Mineral phase analysis identified the volume distribution of individual minerals contained within CSM-LHT-1G.
Plagioclase feldspar is the dominant mineral phase present, and together with augite and anorthosite glass represents 88% by volume.Table 2 lists the primary mineral constituents along with their elastic properties reported from various sources.

Average Grain Density
Average grain density of the regolith simulant is required for accurate porosity estimation of samples.We determined the average grain density using a helium pycnometer with a reported accuracy of 0.03% of the measured reading.Individual sample mass was relatively low (10-14 g); therefore, grain density measurements were taken on 18 separate samples to capture any potential sample variability.Average grain density measurements ranged between 2.83 and 3.04 g/cm 3 with an average of 2.98 g/cm 3 , which is in excellent agreement with pycnometer measurements of Apollo fines (Anderson et al., 1970).

Experimental Configuration
We designed and constructed an experimental apparatus to determine the compressional and shear velocity of unconsolidated granular media at low confining pressure by measuring the ultrasonic time-of-flight through the material (Birch, 1960).Figure 3 shows the experimental configuration.The three primary components of the apparatus are the sample holder, contact transducers, and the sample chamber.The sample holder consists of two precision-milled aluminum end plates connected by steel dowels that serve to support and align the transducers via recesses milled into each plate.The upper plate slides freely along the steel dowels, allowing for easy loading and unloading of transducers and samples within the holder.Additionally, the upper plate was designed with minimal mass, resulting in very low confining pressure (0.005 MPa) being applied to the sample under vertically oriented static conditions.
Two pairs of ultrasonic transducers were used in this analysis: Olympus V601 compressional transducers and V151 shear transducers.Both sets of transducers were constructed with 500 kHz piezoelectric crystals with a nominal element size of 2.5 cm.Vinyl tubing with an inner diameter of 3.3 cm was used for the sample chamber, which fit snugly around the 3.2 cm diameter face plates.The upper edge of the vinyl sample chamber was beveled to ensure that the upper transmitting transducer was only in contact with the granular sample.
Friction between the sliding upper plate and steel dowels was assumed to be negligible since the guide holes in the upper plate were milled 1 mm larger than the diameter of the dowels.The confining pressure being applied to a sample within the sample chamber was calculated based on the surface area of the transducer element in contact with the sample as well as the weight being applied above the transducer element.Combining the weight of the transducer and the upper plate, the calculated baseline confining pressure is 0.005 MPa.By adding calibration masses on top of the upper plate, experimental pressure steps up to 0.08 MPa were achieved.The precise pressure applied to samples was not calibrated due to equipment and environment limitations; future improvements to this work should include precise calibration for confining pressure.
The entire experimental apparatus was kept inside a large walk-in freezer with a nominal temperature of 26°C.Freezing conditions were not strictly required for the current study but were driven by additional experimental goals not presently covered.Coaxial cables connected the transducers to a pulse generator and Tektronix TDS 3014B oscilloscope located outside of the freezer.High-quality waveforms recorded by the oscilloscope meant that a pre-amplifier was unnecessary for this experiment.

Sample Preparation
CSM-LHT-1G simulant was vacuum dried for 24 hr and then stored in an airtight container within the walk-in freezer to minimize any moisture within the samples.Simulant was loaded into the sample chamber using three methods to create porosity variation between samples: (a) gentle air pluviation for high relative porosity, (b) lightly tamping the sample for medium relative porosity, and (c) heavily tamping and tapping the sample for low relative porosity.This process resulted in sample porosities of ∼36%-49%.
The mass of each sample was recorded using a digital balance with 0.01 g readability.Measurement drift of >1 g over 30 min was observed due to the low temperature environment affecting the balance, so efforts were made to warm the equipment between sample measurements.The transducers and sample chamber were loaded into the sample holder and sample length was recorded using a digital caliper with <0.2 mm error, which was determined using a calibration block.

Porosity Determination
We calculated the average porosity (Φ) for each measurement sample as follows: where m is the sample mass, V is the sample volume, and ρ gr is the average grain density.Porosity values presented throughout this study are calculated using an average grain density value of 2.98 g/cm 3 .In Section 3.2, we  (Graf, 1993).Data for sample 15301 were converted to volume percentage assuming uniform particle density.
reported that the average grain density measurements of CSM-LHT-1G ranged between 2.83 and 3.04 g/cm 3 , which would result in ∼±2% porosity uncertainty.
While the maximum calculated sample porosity was 49% after minor compaction at 0.005 MPa, a value of 66% porosity was reported using the helium pycnometer under completely uncompacted conditions.The higher value from the pycnometer also includes a small amount of sample headspace, which would artificially increase the porosity measurement.These values provide bounds for the critical porosity, which is the crossover point above which a material behaves as a suspension (Nur et al., 1995).For this study, a critical porosity value of 60% was chosen, which is above the average estimated porosity of 50%-54% for the upper 15 cm of the lunar subsurface (Carrier et al., 1991).Both the critical porosity and current porosity are inputs to the rock physics modeling presented in Section 5.

Velocity Determination
Compressional and shear wave first arrivals were manually determined for each measurement.Compressional first arrivals were determined by selecting the time at which the recorded signal began to deviate from zero amplitude.The first arrivals from shear measurements are more difficult to identify due to contamination from preceding compressional waves and were determined by locating a slope change or start of a perceived interference pattern that denotes the onset of the shear wave.Due to this increased ambiguity of the shear first arrival, the error for shear wave picking is estimated to be up to 10% as opposed to compressional pick error up to 5%. Figure 4 illustrates representative examples of compressional and shear waveforms.The system delay time was measured by placing the transducer elements in direct contact and identifying the first arrivals.The system delay time was then subtracted from the total measured delay time from the manual picks, which combined with the measured sample length yielded the velocity of the sample.et al. (2023) provide compressional and shear velocity data as well as the recorded waveforms used in this analysis.It should be noted that compressional and shear velocity measurements were not taken on identical samples due to experimental constraints; therefore, the exact porosity values vary between the two measurement sets.Notably, porosity values down to 36% were achieved during compressional measurements, whereas the lowest porosity value for shear measurements was 40%.a Effective mineral was calculated using a Voigt-Reuss-Hill average (Mavko et al., 2020).

Amos
Both the compressional and shear velocity show clear porosity and pressure trends (Figure 5); velocity increases with decreasing porosity or increasing pressure.Some measurement scatter exists due to irreducible heterogeneity between samples, and may be due to the unique packing structure of each sample (Zhai et al., 2020).Although each pressure step represents a consistent pressure increase of 0.025 MPa, it can be observed that velocity does not linearly increase with confining pressure and instead shows a diminishing pressure dependence as pressure is increased.In other words, a significant velocity gradient should be observed over our depth of interest (<10 m).As the confining pressure was increased from 0.005 to 0.08 MPa, porosity was reduced by 1.5%-4.5%,and porosity reduction was directly proportional to the initial porosity of the sample (i.e., samples with higher initial porosity tended to show a larger porosity reduction with pressure).

Rock Physics Modeling
As previously discussed, established CRM struggle to predict the velocity of granular materials at low confining pressure (Bachrach et al., 2000;Zimmer et al., 2007a).Figure 6 illustrates a select few of these models and how they compare to our lowest pressure (0.005 MPa) velocity measurements.The Hertz-Mindlin model (Equation A1) with no grain slippage agrees with only a few of the measured bulk moduli values around 40% porosity; however, most of the data as well as the data trend is not captured.This model also severely overpredicts the shear modulus, leading to overpredicting both V P and V S .The modified CRM proposed by Bachrach et al. (2000) (Equation A2; hereafter referred to as the Bachrach model) agrees with the bulk moduli measurements above 40% porosity but overpredicts shear modulus, V P , and V S , even with a high reduction in shear stiffness of the material compared to the Hertz-Mindlin model.This mismatch was not observed by Caicedo et al. (2023), who used the Bachrach model to describe Martian regolith simulant, perhaps due to the larger rounded grains used in their study.The uncemented sand model (Dvorkin & Nur, 1996, Equation A4) shown in Figure 6 has been modified to use the effective moduli from the Bachrach model at critical porosity, which agrees with measured bulk moduli above 40% porosity but again overpredicting shear modulus, V P , and V S .Although a moderate model fit can be achieved to our V P measurements or bulk modulus by adjusting input parameters such as grain contact radius and coordination number, the combination of bulk modulus, shear modulus, V P , and V S is not matched by any existing CRM.
While the uncemented model does desirably trend to the mineral moduli at zero porosity, the measured bulk modulus shows an increasing (i.e., stiffening) trend below ∼40% porosity, which is not captured by the model.Preparing samples below 40% porosity required heavy tamping and tapping, possibly resulting in the emergence of force chains which would stiffen the grain pack (Peters et al., 2005).A similar compacting effect may be caused in the lunar regolith by impact-induced vibrations, evidenced by the bulk density in the upper 3 m of the lunar subsurface being higher than anticipated based on self-compaction (Carrier et al., 1991).We expect porosity values above 25% in the lunar near-surface (Carrier et al., 1991), so capturing this stiffening behavior is critical for our model.
Observing that existing CRMs do not effectively capture our measured data (namely shear modulus), a new method is needed to model the elastic properties of lunar regolith simulants.To this end, we have developed a new workflow that draws from existing models and makes heuristic modifications to fit our data.Our approach addresses model inaccuracies due to uncertainty in CRM input parameters and model assumptions, and specifically corrects for our observed low shear modulus values.Additionally, we introduce a model modification to account for the sharp increase in velocity observed below 40% porosity.We accomplish this by adopting the modified CRM model from Bachrach et al. (2000) to calculate the bulk modulus at critical porosity and use a lower Hashin-Shtrikman (HS;Hashin & Shtrikman, 1963) bound similar to the uncemented sand model (Dvorkin & Nur, 1996) to calculate the bulk modulus below critical porosity.Additionally, we adopt a gradual isoframe model approach  (Fabricius, 2003) to capture the stiffening trend below 40% porosity.We then calculate shear modulus based on our modeled bulk modulus and Poisson's ratio observations from our measured data as where K is bulk modulus, μ is shear modulus, and ν is Poisson's ratio.We can then calculate V P and V S from our modeled effective moduli and bulk density (Equation 1).
Specifically, our rock physics modeling approach is as follows: 1. Use the Bachrach model to calculate the effective moduli at critical porosity (Equation A2). 2. Fitting a lower HS bound to this critical porosity point, calculate the effective bulk modulus as porosity is reduced to a point we term the transition porosity (Equation A4). 3. Below the transition porosity, calculate the effective bulk modulus based on a gradual trend from the HS lower bound at the transition porosity to the HS modified upper bound at zero porosity (Equation A3). 4. Calculate a fitting function between the grain Poisson's ratio at zero porosity and the measured Poisson's ratio at experimental porosity.5. Use the Poisson's ratio fit function with the modeled bulk modulus to calculate the shear modulus (Equation 3).6. Repeat steps 1-5 for each experimental pressure step.7. Based on the results from Step 6, calculate a fitting function to describe how Poisson's ratio varies with pressure.
This procedure produces a model which is calibrated to experimental observations and describes elastic behavior over a range of porosities and confining pressures.Figure 7 shows the results of this modeling approach against our measured data, which adequately captures the observed moduli and velocities, including the observed increase in compressive stiffness below the 40% transition porosity.Due to experimental limitations, the shear modulus is not calibrated below this transition porosity.
The Poisson's ratio fitting functions are described in more detail in Appendix B. By relying on observed Poisson's ratio values to calibrate our fitting function, we are able to match the low measured shear modulus values while keeping our modeled shear modulus coupled with our modeled bulk modulus.This results in rational V P /V S ratio values over a wide porosity range (Figure 7e).
Implications of this model may be more readily understood by considering the V P /V S ratio (Figure 7e).At low pressure and high porosity, the V P /V S ratio is high (up to 4.5) and reduces to the mineral V P /V S ratio at zero porosity.As the confining pressure is increased, there is a non-linear decrease in the V P /V S ratio.This implies that tangential stiffness increases more rapidly with confining pressure than normal stiffness.

Discussion
We now compare our model that was calibrated with lunar regolith simulant to reported velocity from both laboratory and in situ experiments from the Apollo program.

Comparison With Apollo Regolith Samples
Figure 8 shows V P data that we estimated from published figures in Johnson et al. (1982) along with our modeled V P .While porosity was not reported by Johnson et al. (1982), using an average grain density of 3.0 g/cm 3 along with their reported initial and final bulk densities suggests that their experimental porosities ranged from 37% to 57%.Our V P model and measurements are slightly higher than most V P values reported by Johnson et al. (1982); however, there is some agreement at both low and high pressures.No clear initial bulk density dependence is observed in the Johnson et al. (1982) measurements, but some hysteresis is observed based on the compaction history of the sample; the points representing the pressure loading cycle (unfilled symbols) are generally lower than measurements from the pressure unloading cycle (solid symbols).This hysteresis leads to these data displaying a variable pressure dependence based on compaction history.Although there is considerable variability in velocity versus pressure, the data are generally explained with pressure exponents between 0.25 and 0.33.Our estimated pressure exponent for V P is 0.2, which is slightly lower than that observed from the Johnson et al. (1982) data.Prasad and Meissner (1992) observed that effective elastic moduli decreased with decreasing grain size and increasing grain angularity.The grain size distribution of Apollo sample 15301 (Graf, 1993) used in the Johnson et al. study is shown in Figure 1b alongside the CSM-LHT-1G grain size distribution.Sample 15301 appears more heavily weighted to finer grain sizes compared to CSM-LHT-1G, which likely contributes to the lower V P observed by Johnson et al. (1982) compared to our measurements.Detailed grain angularity analysis is not available to determine if grain shape could be a contributing factor to the velocity discrepancy observed between CSM-LHT-1G and lunar regolith but should be considered in future studies.
Velocity data measured by Talwani et al. (1974) at atmospheric pressure is also presented in Figure 8.The reported V P generally agrees with measurements reported by Johnson et al. (1982) but is slightly higher than our V P model.V S is slightly higher than our model but appears to agree with modeled V S from Zimmer et al. (2007a).
One potential contributing factor to the observed velocity discrepancy between our measurements and those of Johnson et al. (1982) is experimental error.It was observed that the vinyl sample chamber became very stiff under our cold experimental conditions, potentially providing a low-loss medium with a higher velocity than the CSM-LHT-1G sample.This could result in an inaccurately high measured velocity if the ultrasonic waves refracted along the sample chamber.Additional experimental configurations and materials should be included in future work to test this hypothesis.
We compare these velocity results of actual and simulated lunar regolith to measurements of dry terrestrial sands in Figure 8.Our V P measurements and resulting model agree with several low pressure measurements by Prasad (2002), although there is a group of measurements at roughly 0.8 MPa that falls below our modeled velocity.This reduced velocity is attributed to higher grain angularity in these samples, suggesting that our model provides a closer match to sands with rounded grains at pressures higher than our experimental range.Our modeled V S is slightly higher than data from Prasad (2002) above 0.6 MPa.We also compare our model with V P and V S calculated with equations from Zimmer et al. (2007a), which show close agreement to our modeled V P at low pressures but higher V S than our model.CSM-LHT-1G contains finer and more angular grains than sands measured by Zimmer et al. (2007a), which would reduce the shear stiffness of the grain pack and produce our lower observed shear velocities.

Depth Profile Modeling
Our model can be used to calculate 1D velocity versus depth profiles for estimated lunar near-surface conditions.Key inputs to create these velocity profiles are density, porosity, and pressure profiles with depth.
Bulk density has been measured from several core samples that were retrieved during the Apollo program (Carrier et al., 1991), but no density calibration exists deeper than 3 m.In situ bulk density is relatively poorly constrained due to sample disturbance from the coring and transportation process however Carrier et al. (1991) offer two equations to relate bulk density (ρ) to depth (z) given the range of available core measurements: Additionally, the studies compiled by Carrier et al. (1991) provide a range of regolith grain density between 2.9 and 3.24 g/cm 3 .Combining these two bulk density equations with the grain density range, we can calculate a range of porosity values with depth (Figure 9a).Our measured bulk density values are slightly lower than those estimated by these depth-density relationships, suggesting that in situ lunar regolith may be slightly more compacted at a given depth/pressure than our measured samples.Lastly, we calculate confining pressure (P) with depth with where g is lunar acceleration due to gravity (1.625 m/s 2 ) and ρ is the depth-dependent bulk density.
Using these density, porosity, and pressure profile ranges, we use our rock physics model to calculate 1D instantaneous V P and V S profiles with depth, as shown in Figure 9b to a depth of 25 m.To compare with published average velocity values interpreted from the Apollo active seismic experiments, we compute average velocity profiles (Figure 9c) from these instantaneous velocity profiles.Much of our average V P and V S ranges agree with Journal of Geophysical Research: Planets 10.1029/2024JE008287 the 100-140 m/s and 35-100 m/s ranges interpreted from Apollo experiments.It is important to note that our experimental conditions only produced one group of measurements at a confining pressure equivalent to being within the upper 10 m of the lunar surface, and our model is uncalibrated at pressures shallower than 7 m depth equivalent.Future work should attempt to calibrate velocity under these extremely low confining pressures to better predict very near-surface seismic velocity.

Comparison With Apollo Seismic Experiments
To further evaluate our model, we used the open-source pyGIMLi package (Rücker et al., 2017) to build a 2D V P profile based on our velocity-depth relations and simulate the Apollo 14 and 16 active seismic experiments.Our model provides a good fit to measurements for bulk and shear moduli (a), (b) as well as V P and V S (c), (d).V P /V S ratio (e) and Poisson's ratio (f) are calculated using V P data and a least squares polynomial fit to V S data since both V P and V S measurements could not be taken from identical samples.
Journal of Geophysical Research: Planets 10.1029/2024JE008287 Calculating the shortest traveltime path through this V P model (Moser, 1991) allows us to compare simulated first arrivals to those interpreted from Apollo seismic records.Although the regolith mineralogic composition was different at the Apollo 14 and 16 landing sites, Blangy et al. (1993) notes that mineralogy is not a dominant factor in elastic moduli of poorly consolidated material.This observation coupled with relatively similar seismic velocity interpreted at both landing sites gives us confidence that mineralogy differences within the regolith are not a main contributing factor to any seismic velocity variation in the near-surface.
To build our V P profile, we assume a uniform layer of regolith with depthdependent velocity overlying a layer with a constant higher velocity (∼300 m/s).The thickness of the lunar regolith layer estimated from in situ seismic measurements at Apollo 14 and 16 landing sites ranges from roughly 9-12 m (e.g., Cooper et al., 1974;Horvath et al., 1980) and is treated as an uncertainty in our model.For the seismic receivers, three geophones are simulated with 45.7 m spacing as in the Apollo experiments.Rather than simulate ∼20 seismic shots to mimic the Apollo experiments, we simulate a very dense shot spacing to produce a continuous estimation of first arrivals versus distance between shots and receivers (offset).Figure 10 illustrates the shortest raypaths through our V P model between a single seismic shot and the three receivers.Note the curved raypaths through the shallow regolith layer as well as the refracted raypath at the top of the deeper constant-velocity layer.Martella et al. (2022) presented an analysis of the Apollo 14 and 16 active seismic data, which is compared against our simulated first arrivals (Figure 11).The best match between our modeled range of V P and the first arrivals from Martella et al. comes from our slowest V P model.This velocity model results in slightly larger traveltime (lower average velocity) at short offsets and slightly lower traveltime (higher average velocity) at further offsets; however, our model appears to be in overall good agreement with the interpreted Apollo data.This same trend can be observed in Figure 9c.
A ∼0.17 depth exponent provides the best fit to first arrival data from Martella et al. (2022), whereas our calibrated rock physics model yields depth exponents of 0.42-0.45under our estimated lunar subsurface conditions.This is different from our previously mentioned pressure exponent of 0.2 since pressure does not increase linearly with depth in our lunar subsurface model due to porosity reduction.The difference between V P depth exponents from our model and interpreted Apollo data suggests that we are overestimating the pressure dependence of velocity in the near-surface.Previous studies have shown a reduced pressure dependence in unconsolidated sands that had been previously loaded under pressure (Vega et al., 2006;Zimmer et al., 2007a), suggesting that irrecoverable porosity reduction may lower the pressure dependence of a grain pack.Additionally, Figure 9a shows that our experimental porosity range at 0.005 MPa (∼7 m depth) is mostly higher than the estimated lunar subsurface conditions, suggesting that our measurements taken during sample loading represent porosities higher than anticipated for the lunar near-surface.While there is no feasible mechanism to bury lunar regolith at high pressure/depth and then exhume this regolith to the surface in an intact state, it has been observed that the shallow lunar regolith displays higher relative density than can be explained by compaction from self-weight alone and may be caused by repeated impact events and natural seismic quakes (Carrier et al., 1991).Such a densification would result in a similar reduced-porosity grain pack as pressure cycling and could reduce V P pressure dependence in the near-surface regolith.Future studies should include additional porosity reduction techniques and velocity measurements during pressure unloading to further calibrate the velocity-pressure dependence of simulated lunar materials.
Another potential source of discrepancy between our model and Apollo data is moisture.Although the velocity measurements used to calibrate our model were made on vacuum-dried samples, any remaining moisture on grain surfaces may impact the recorded velocity (e.g., Lisabeth et al., 2023).Previous studies have shown that extensive outgassing procedures to remove adsorbed volatiles from rock samples result in seismic quality factor measurements in agreement with Apollo seismic data (Tittmann, 1977;Tittmann et al., 1974); such extensive procedures were not used in this study.Future work should focus on the effect of low moisture content in regolith simulant to quantify any moisture-induced discrepancy between laboratory and lunar in situ data.

Future Model Applications
With the encouraging match between our modeled near-surface velocity and Apollo seismic observations, we can use our model to predict the seismic response of untested lunar subsurface scenarios.For instance, with water ice being present near the lunar poles (Colaprete et al., 2010), we can begin to model how seismic velocity may change with varying ice concentration or ice form (i.e., granular vs. cementing ice).Such modeling would help describe the effectiveness of seismic methods for lunar ice resource exploration and will be the focus of a future study.

Conclusions
We developed a rock physics model to describe the elastic behavior of a fine-grained lunar regolith simulant at variable porosity and confining pressure and calibrated this model with laboratory ultrasonic velocity measurements.Forward-modeling seismic traveltime using our velocity-depth relations yields a good agreement with observed Apollo near-surface seismic velocity.While our observed velocity-pressure dependence agrees with other published studies on unconsolidated sands, we note a slight difference in pressure dependence compared with interpreted Apollo first arrivals.This difference in pressure dependence may be due to porosity reduction in the lunar regolith caused by vibrations from impact events.
Future work may focus on improving our experimental technique and conditions, such as testing the effect of apparatus materials on velocity measurements at low temperatures, exploring additional porosity reduction procedures, and performing measurements at even lower confining pressures to calibrate elastic properties of regolith at <7 m depth equivalent.Additional analysis and comparison of CSM-LHT-1G characteristics (e.g., grain angularity) versus Apollo regolith samples may yield further insight into any observed velocity discrepancies.Experimental errors should also be analyzed and propagated through the rock physics modeling results to understand the velocity uncertainty associated with this model.
Finally, as our current calibrated rock physics model provides a level of confidence by producing a good agreement with observed Apollo seismic velocities, future work should focus on perturbing this model to predict seismic velocity under varying conditions; for example, predicting how ice may impact near-surface lunar seismic velocity.π 2 (1 ν gr ) 2 ) where K eff and μ eff are the effective bulk and shear moduli, respectively.R/R relates the effective contact radius to grain radius and f t denotes the fraction of non-slipping contacts.Adjusting these parameters can yield a lower effective shear stiffness and increase the predicted V P /V S ratio.
In addition to models which describe the effective moduli of a granular material, some model modifications incorporate the HS bounds (1963), which describe the maximum and minimum possible effective moduli.For a two-phase composite, these bounds are given by Mavko et al. (2020): where f 1 and f 2 are the volume fractions of the two phases, K 1 and K 2 are the bulk moduli of the two phases, and μ 1 and μ 2 are the shear moduli of the two phases.K m and μ m are the maximum or minimum phase moduli; the upper bound (K HS+ , μ HS+ ) is computed using the maximum moduli, and the lower bound (K HS , μ HS ) is computed using the minimum moduli.
The uncemented sand model proposed by Dvorkin and Nur (1996) is based on the observation that velocities in unconsolidated sandstones can be described by fitting a lower HS bound between the mineral moduli at zero porosity and the Hertz-Mindlin effective moduli at the critical porosity.This model is given by where Φ c is the critical porosity, and K HM and μ HM are the effective moduli given by the Hertz-Mindlin equations.
A modified upper HS bound can also be calculated in a similar fashion as the uncemented model by interchanging the grain moduli in place of the Hertz-Mindlin moduli.Fabricius (2003) discussed the isoframe model, where moduli trend from the Reuss bound (the "softest" possible bound describing a suspension of grains) toward the modified upper HS bound as sediments become more compacted and cemented.Between these bounds, it is assumed that some fraction of grains are in contact and the remaining fraction are in suspension, that is, not contributing to the stiffness of the material.

Figure 2 .
Figure 2. Optical microscope image of CSM-LHT-1G Illustrating an assortment of grain shapes and sizes.Each black scale division is 100 μm.Dark green grains are anorthosite glass particles.

Figure 3 .
Figure 3. Experimental setup inside freezer.Two transducers (a) are held in a vertical orientation separated by the sample chamber (b) and secured by the sliding end plate (c).Calibration masses (d) were placed on top of the sliding end plate to increase the confining pressure applied to the sample.

Figure 4 .
Figure 4. (a) Compressional and (b) shear waveforms under variable confining pressure.Note the very low amplitude values at 0.005 MPa make these waveforms hard to discern when compared to higher pressure measurements; therefore, the dotted black lines show these waveforms at 10× amplitude for increased visibility.

Figure 5 .
Figure 5. (a) V P and (b) V S measurement results over variable confining pressure and porosity.

Figure 6 .
Figure 6.Elastic moduli (a), (b) and velocity (c), (d) of select rock physics models versus our model.Black circles are from our lowest pressure (0.005 MPa) measurements.Variable agreement is observed between models and measurements for bulk modulus (a); however, our model provides a significantly improved match to the observed shear modulus (b) and velocities (c), (d).

Figure 7 .
Figure7.Elastic moduli and velocities of our rock physics model (lines) versus measurements (points) under variable confining pressure.Our model provides a good fit to measurements for bulk and shear moduli (a), (b) as well as V P and V S (c), (d).V P /V S ratio (e) and Poisson's ratio (f) are calculated using V P data and a least squares polynomial fit to V S data since both V P and V S measurements could not be taken from identical samples.

Figure 8 .
Figure8.Modeled V P and V S versus pressure.V P measurements of lunar regolith fromJohnson et al. (1982) are shown by black symbols corresponding to their initial bulk density (square = 1.29 g/cm 3 , triangle = 1.39 g/cm 3 , circle = 1.53 g/cm 3 ).Velocity measurements of lunar regolith fromTalwani et al. (1974) are shown by colored diamonds, and velocity measurements of terrestrial sands fromPrasad (2002) are shown by colored squares.V P and V S calculated from equations given inZimmer et al. (2007a) for dry terrestrial sands are shown by dashed lines.

Figure 9 .
Figure 9. (a) Calculated bulk density, porosity, and confining pressure profiles.Experimental bulk density and porosity are shown by horizontal bars at their respective confining pressure (depth), (b) Vertical V P and V S velocity ranges from our rock physics model.A one-sixth power velocity-depth fit to Apollo V P first arrivals from Martella et al. (2022) is illustrated by the black dashed line, (c) Computed average velocity profiles.

Figure 10 .
Figure 10.V P profile of the lunar near-surface produced using our rock physics model.Black lines illustrate the shortest seismic raypath between the seismic source and the three receivers.The displayed configuration would be representative of shot 5 from the Apollo 16 seismic survey.Note the horizontal refracted raypath along the top of the second layer at 12 m depth.

Figure 11 .
Figure 11.V P first arrivals from our velocity model (blue) versus interpreted first arrivals from Apollo 14 and 16 (gold) fromMartella et al. (2022).Black dotted lines represent average velocity bounds of 100-140 m/s for the regolith layer, and the vertical black line represents the maximum horizontal offset distance associated with the first arrivals of the regolith based on our velocity model.Our model is in good agreement with the interpreted Apollo shallow velocity but displays a lower average V P over the first several m horizontal offset and underpredicts refracted first arrivals from the deeper layer.

Table 1
In Situ and Laboratory Seismic Velocity Measurements