Stochastic Model of Regolith Gardening

• The paper demonstrates that breccia-lens infilling dominates regolith production over ejecta blanketing in small lunar craters.
• It employs a probabilistic framework to quantify regolith thickness as a function of random, discrete impact events using analytical methods.
• The model reconciles Apollo 15 seismic data with predicted spatial variability in regolith accumulation and thickness.

Searching arXiv for relevant papers on stochastic models of regolith gardening and adjacent observational/constitutive constraints. A stochastic model of regolith gardening treats regolith evolution on airless bodies as the cumulative outcome of discrete, spatially random impact events that excavate, bury, and overlap one another through time. In the formulation developed for small, simple lunar craters by Hirabayashi et al., the central quantity is not a single bulk regolith thickness but the probability distribution of thickness at a point, expressed as the fraction of area whose regolith thickness exceeds a specified depth $h$ (Hirabayashi et al., 2018). Within this framework, “gardening” is modeled primarily as regolith production and spatial accumulation by impact cratering, with ejecta blanketing outside crater rims and breccia-lens infilling inside crater cavities as the two explicit regolith-generating processes (Hirabayashi et al., 2018). The resulting picture is intrinsically heterogeneous: regolith thickness varies laterally because impacts occur randomly in space, and continued cratering can thicken regolith even after crater equilibrium is reached (Hirabayashi et al., 2018).

1. Conceptual definition and scope

In the stochastic formulation of regolith gardening, repeated impacts into an initially intact, flat surface generate a thickness field that is probabilistic rather than spatially uniform. The modeled processes are excavation and deposition of ejecta blankets, regolith infilling of the transient crater cavity represented as a breccia lens beneath the final crater floor, overlap and superposition of multiple impacts, lateral variability in regolith thickness, and temporal evolution of regolith thickness as crater number increases through time (Hirabayashi et al., 2018).

This scope is deliberately narrower than many broader uses of the term “gardening.” The model does not explicitly include topographic diffusion of crater forms, regolith mixing or turnover rates in the sense of vertical overturn statistics, particle size evolution, compaction, melting, shock-damage halos outside the adopted crater-regolith units, or spatially heterogeneous distal rays (Hirabayashi et al., 2018). A common misconception is that a regolith-gardening model necessarily constitutes a complete maturity model. In this case it does not: the formulation is directed at impact-generated regolith production and accumulation by small, simple craters rather than grain comminution, agglutinate production, or optical maturity (Hirabayashi et al., 2018).

The stochastic aspect is central because deterministic average-thickness models cannot represent the fact that thickness at a given point depends on the random history of crater overlaps of different sizes and positions (Hirabayashi et al., 2018). This is closely related, but not identical, to later continuum treatments in which discrete impact statistics are coarse-grained into effective vertical transport. A 2026 lunar model recasts gardening as a competition between impact-driven advection and stochastic diffusion in depth, arguing that pure diffusion cannot reproduce sharp gradients and displaced centroids in maturity and radionuclide profiles (Costello et al., 10 Apr 2026). This suggests that “stochastic regolith gardening” has at least two complementary meanings in current literature: lateral stochastic coverage by crater-generated regolith units (Hirabayashi et al., 2018), and stochastic-impact-informed continuum transport in depth (Costello et al., 10 Apr 2026).

2. Crater geometry, assumptions, and single-event regolith production

The analytical model of Hirabayashi et al. is restricted to populations of small, simple craters (Hirabayashi et al., 2018). Each crater contributes two idealized regolith units: an ejecta blanket outside the final rim and a breccia lens inside the crater cavity. For the Apollo 15 application, craters with radius $r \le 381\ \mathrm{m}$ contribute via ejecta blanketing, whereas craters with radius $r \le 100\ \mathrm{m}$ contribute via breccia-lens infilling (Hirabayashi et al., 2018). The asymmetry reflects the distinction between local cavity infilling and ejecta contributions that can come from larger or more distant craters.

The model assumes that impacts occur randomly on an initially intact, flat surface; newly formed craters are circular; each crater produces regolith only in the ejecta blanket and breccia lens; no distinction is made between regolith particle sizes in ejecta and breccia lens; and no melting, compaction, topographic evolution, explicit mixing law, thermal fatigue, or volcanism is included in the analytical formulation (Hirabayashi et al., 2018). Crater production follows a cumulative size-frequency distribution of power-law form, crater overlap is treated statistically rather than geometrically, the minimum crater radius is taken effectively to zero for analytical simplicity, and maximum contributing crater radii are imposed separately for infilling and ejecta blanketing (Hirabayashi et al., 2018). The model also does not assume that crater equilibrium halts regolith growth (Hirabayashi et al., 2018).

Outside the crater, ejecta thickness is prescribed by a power law,

$$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$

where $r_c$ is the final crater radius, $\sigma$ is the ejecta thickness ratio at the rim, and $\kappa$ is the ejecta thickness power-law slope (Hirabayashi et al., 2018). For most of the paper, the adopted values are $\sigma=0.014$ and $\kappa=3$, following Sharpton (2014) (Hirabayashi et al., 2018).

Inside the crater, breccia-lens thickness is approximated by

$h_{\mathrm{in}}= \begin{cases} \delta r_c\left(1-\dfrac{r^2}{r_c^2}\right), & r\le r_c,\[6pt] 0, & r>r_c, \end{cases}$

where $r \le 381\ \mathrm{m}$0 is the ratio of maximum breccia-lens thickness at crater center to crater radius (Hirabayashi et al., 2018). For lunar fresh small simple craters, the inferred range is $r \le 381\ \mathrm{m}$1, and for the Apollo 15 application the adopted value is $r \le 381\ \mathrm{m}$2 (Hirabayashi et al., 2018). The physical significance of this geometry is that the breccia lens is much thicker than the ejecta blanket near crater center (Hirabayashi et al., 2018).

3. Probabilistic formulation of regolith thickness

The stochastic model asks what fraction of area has regolith thickness greater than a threshold $r \le 381\ \mathrm{m}$3 (Hirabayashi et al., 2018). For craters of a single size $r \le 381\ \mathrm{m}$4, the area fraction at depth $r \le 381\ \mathrm{m}$5 affected by one crater is

$r \le 381\ \mathrm{m}$6

where $r \le 381\ \mathrm{m}$7 is the area at depth $r \le 381\ \mathrm{m}$8 affected by a crater of radius $r \le 381\ \mathrm{m}$9, and $r \le 100\ \mathrm{m}$0 is the total target area (Hirabayashi et al., 2018). If $r \le 100\ \mathrm{m}$1 craters of that size are randomly emplaced, the empty fraction approaches the Poisson limit

$r \le 100\ \mathrm{m}$2

and extension to all crater sizes yields

$r \le 100\ \mathrm{m}$3

where $r \le 100\ \mathrm{m}$4 is the fraction of surface area where regolith thickness exceeds $r \le 100\ \mathrm{m}$5 (Hirabayashi et al., 2018).

This quantity has several equivalent interpretations. It is an areal thickness exceedance function, the survival function of regolith thickness at a randomly chosen point, and a probabilistic description that naturally incorporates crater overlap (Hirabayashi et al., 2018). Because the formulation returns a full thickness distribution rather than only a mean, it predicts that some locations are buried beneath thick breccia lenses and multiple ejecta blankets while others are affected only by small thin deposits (Hirabayashi et al., 2018).

Crater production is represented through a cumulative size-frequency distribution

$r \le 100\ \mathrm{m}$6

where $r \le 100\ \mathrm{m}$7 is the normalization coefficient, $r \le 100\ \mathrm{m}$8 is a dimensionless normalized crater-production number proportional to exposure time for constant impact flux, and $r \le 100\ \mathrm{m}$9 is the cumulative power-law slope (Hirabayashi et al., 2018). In the Apollo 15 application, the inferred parameters are $$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$0, $$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$1, and $$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$2 for the present produced crater CSFD (Hirabayashi et al., 2018).

For a given threshold $$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$3, the area affected by a single crater is computed separately for ejecta and breccia-lens infill. For ejecta,

$$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$4

and for the breccia lens,

$$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$5

for $$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$6 (Hirabayashi et al., 2018). These relations are then integrated over crater sizes to obtain the continuous stochastic coverage law (Hirabayashi et al., 2018).

Because the model yields a thickness distribution over area, the expected regolith thickness for a randomly selected point is defined by

$$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$7

This expected value is the quantity compared with empirical constraints such as Apollo seismic regolith thickness (Hirabayashi et al., 2018).

4. Breccia lenses as the dominant regolith source

The central conclusion of the analytical model is that breccia-lens infilling dominates impact-generated regolith production (Hirabayashi et al., 2018). In the adopted geometry, the breccia lens is a paraboloid whose maximum thickness at crater center is $$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$8, tapering quadratically to zero at the rim (Hirabayashi et al., 2018). The parameter $$h_{\mathrm{out}} \sim \sigma r_c \left(\frac{r}{r_c}\right)^{-\kappa}, \qquad r>r_c,$$9 is inferred from crater geometry using a transient crater depth/diameter ratio of $r_c$0, a final crater diameter $r_c$1 times transient diameter, and observed final crater depth-to-diameter ratios of $r_c$2 to $r_c$3, leading to $r_c$4 for final $r_c$5 and $r_c$6 for final $r_c$7 (Hirabayashi et al., 2018).

At Apollo 15, the quantitative contrast between ejecta blanketing and breccia-lens production is strong. Using local crater counts, the ejecta-only expected thickness is $r_c$8, whereas the ejecta-plus-breccia-lens expected thickness is $r_c$9 locally (Hirabayashi et al., 2018). Sensitivity tests reinforce this conclusion: varying the ejecta thickness parameter $\sigma$0 has little effect, whereas sensitivity to $\sigma$1 is strong (Hirabayashi et al., 2018). The reported $\sigma$2-dependence is $\sigma$3 for $\sigma$4, $\sigma$5 for $\sigma$6, and $\sigma$7 for $\sigma$8 (Hirabayashi et al., 2018).

Volume arguments point in the same direction. For the nominal Apollo 15 parameters, total regolith volume generated by one crater is partitioned approximately as 15% ejecta and 85% breccia lens (Hirabayashi et al., 2018). Even under an approximate volume-conserving adjustment of the rim and ejecta assumptions, the partition becomes 42% ejecta and 58% breccia lens, while predicted total regolith thickness remains nearly unchanged: $\sigma$9 for the volume-conserving case versus $\kappa$0 for the original case (Hirabayashi et al., 2018).

A frequent misconception in earlier crater-gardening discussions is that ejecta blanketing alone explains observed lunar regolith thicknesses. The Hirabayashi model explicitly argues the opposite: ejecta blanketing alone is far too weak, and breccia-lens formation controls total regolith production (Hirabayashi et al., 2018). This conclusion is specific to the adopted crater class and regolith units, but within that domain it is the decisive result.

5. Apollo 15 calibration, regolith growth, and lunar thermophysical context

The Apollo 15 application computes regolith thickness using the observed crater size-frequency distribution of small, simple lunar craters, with $\kappa$1 for ejecta blanketing and $\kappa$2 for regolith infilling (Hirabayashi et al., 2018). Allowing for some amount of regolith coming from outside the area, the modeled result is consistent with the empirical result from the Apollo 15 seismic experiment (Hirabayashi et al., 2018). The paper also concludes that the timescale of regolith growth is longer than that of crater equilibrium, implying that crater equilibrium on a surface does not imply a fixed regolith thickness (Hirabayashi et al., 2018).

This distinction matters for interpretation of lunar surface states. Crater equilibrium is a statement about crater populations, whereas the stochastic regolith-thickness distribution continues to evolve through additional impacts (Hirabayashi et al., 2018). A plausible implication is that surfaces classified as saturated in crater counts can still be evolving in subsurface physical structure.

Independent lunar thermophysical observations provide a complementary near-surface constraint. Diviner mapping shows that the upper regolith fines layer is globally remarkably uniform, with globally averaged thermal inertia at $\kappa$3 of $\kappa$4, global mean $\kappa$5, and global-scale variations over $\kappa$6 scales of $\kappa$7 (Hayne et al., 2017). Hayne et al. interpret this uniformity as evidence that the upper $\kappa$8 has been rapidly homogenized by impact gardening and lateral mixing on timescales $\kappa$9 (Hayne et al., 2017). The diurnally active near-surface layer is $\sigma=0.014$0, and density increases from $\sigma=0.014$1 at the surface to $\sigma=0.014$2 at about $\sigma=0.014$3 depth (Hayne et al., 2017).

The stochastic thickness model and the Diviner thermophysical results refer to different aspects of regolith structure. The former addresses cumulative production and lateral variability at meter scales and greater (Hirabayashi et al., 2018), whereas the latter constrains the thermally active fines layer in the upper few centimeters (Hayne et al., 2017). Taken together, they indicate that laterally heterogeneous regolith production at depth is compatible with large-scale homogenization of the very shallow fines layer.

6. Extensions, analogues, and limitations

A later lunar development generalizes stochastic gardening into an effective one-dimensional transport model in which impact statistics produce a net gardening depth $\sigma=0.014$4, an effective advection velocity

$\sigma=0.014$5

and a diffusion coefficient

$\sigma=0.014$6

which are then inserted into an advection-diffusion equation for tracer concentration (Costello et al., 10 Apr 2026). In that formulation, gardening is not pure diffusion but the competition between downward burial and stochastic local mixing (Costello et al., 10 Apr 2026). The model is validated against Apollo maturity profiles over $\sigma=0.014$7 to $\sigma=0.014$8 years and applied to $\sigma=0.014$9Fe and $\kappa=3$0Pu depth profiles (Costello et al., 10 Apr 2026). Relative to the lateral overlap model of Hirabayashi et al., this is a vertical coarse-grained transport closure rather than a pointwise areal-thickness distribution (Hirabayashi et al., 2018, Costello et al., 10 Apr 2026).

Analogous gardening concepts also appear outside the lunar case. On Vesta, photometric mapping with Dawn reveals a phase-curve slope parameter $\kappa=3$1 that is interpreted as a proxy for regolith physical roughness (Schröder et al., 2017). High $\kappa=3$2 occurs in ejecta of young craters such as Antonia and Cornelia, while low $\kappa=3$3 correlates strongly with steep slopes and likely mass wasting (Schröder et al., 2017). Older crater ejecta are not photometrically extreme, leading to the inference that rough ejecta are smoothed over several tens of Myr, most likely by impact gardening (Schröder et al., 2017). This suggests that on Vesta gardening is used less as a thickness-production model than as a roughness-relaxation process.

For S-complex asteroids, a two-process weathering-plus-gardening model represents gardening as the conversion of weathered surface back into fresh surface on a timescale $\kappa=3$4, with best-fit values $\kappa=3$5 and $\kappa=3$6 (Willman et al., 2010). The corresponding mean-field equation,

$\kappa=3$7

has a natural stochastic reading as a two-state patch process, although that interpretation is explicitly identified in the source synthesis as an inference rather than a statement of the paper itself (Willman et al., 2010). The same work notes a discrepancy between the fitted gardening time and an impact-derived estimate of $\kappa=3$8, proposing a “honeycomb” or armoring mechanism in which small impactors are inefficient resurfacers (Willman et al., 2010). This controversy illustrates that effective gardening rates may depend strongly on unresolved assumptions about crater formation, ejecta efficiency, and surface structure.

A further limitation arises from constitutive state. The static stratification model of Schräpler et al. does not describe impact flux or overturn rates, but it provides a pressure-controlled compaction law,

$\kappa=3$9

with $h_{\mathrm{in}}= \begin{cases} \delta r_c\left(1-\dfrac{r^2}{r_c^2}\right), & r\le r_c,\[6pt] 0, & r>r_c, \end{cases}$0 for solid grains (Schräpler et al., 2015). In a gardening context, this functions as a closure for post-disturbance porosity structure rather than a gardening law itself (Schräpler et al., 2015). A plausible implication is that stochastic overturn and burial models require a separate constitutive description of how redeposited material compacts under self-weight.

Across these formulations, the phrase “stochastic model of regolith gardening” therefore denotes a family of related but nonidentical approaches. In the narrow sense established by Hirabayashi et al., it is an analytical description of the probability distribution of impact-generated regolith thickness produced by random crater emplacement, with breccia-lens infilling as the dominant source term (Hirabayashi et al., 2018). In broader current usage, it can also refer to effective stochastic transport in depth (Costello et al., 10 Apr 2026), roughness evolution under impact reworking (Schröder et al., 2017), or mean-field resurfacing of surface-state fractions (Willman et al., 2010). The common element is that discrete impacts are treated as random events whose ensemble behavior governs the structure, thickness, or physical state of regolith on airless bodies.