Numerical simulations of impact crater formation with dilatancy

Impact‐induced fracturing creates porosity that is responsible for many aspects of the geophysical signature of an impact crater. This paper describes a simple model of dilatancy—the creation of porosity in a shearing geological material—and its implementation in the iSALE shock physics code. The model is used to investigate impact‐induced dilatancy during simple and complex crater formation on Earth. Simulations of simple crater formation produce porosity distributions consistent with observations. Dilatancy model parameters appropriate for low‐quality rock masses give the best agreement with observation; more strongly dilatant behavior would require substantial postimpact porosity reduction. The tendency for rock to dilate less when shearing under high pressure is an important property of the model. Pressure suppresses impact‐induced dilatancy: in the shock wave, at depth beneath the crater floor, and in the convergent subcrater flow that forms the central uplift. Consequently, subsurface porosity distribution is a strong function of crater size, which is reflected in the inferred gravity anomaly. The Bouguer gravity anomaly for simulated craters smaller than 25 km is a broad low with a magnitude proportional to the crater radius; larger craters exhibit a central gravity high within a suppressed gravity low. Lower crustal pressures on the Moon relative to Earth imply that impact‐induced dilatancy is more effective on the Moon than Earth for the same size impact in an initially nonporous target. This difference may be mitigated by the presence of porosity in the lunar crust.


Introduction
One of the most important collateral effects of impact cratering on planetary surfaces is fracturing and fragmentation of the target rocks surrounding the crater. Impact-induced fracturing increases the porosity and permeability of the cratered target, which has important implications for postimpact hydrothermal activity [Kirsimäe and Osinski, 2012], fluid (e.g., hydrocarbon) migration [Grieve, 2005], and possible microbial colonization [Cockell et al., 2012]. In addition, fracturing and brecciation are responsible for many aspects of the geophysical signature of an impact crater, including the most characteristic feature: a circular negative gravity anomaly centered over the crater [Pilkington and Grieve, 1992]. Impact-induced fracturing and brecciation is caused by both the passage of the shock wave and the subsequent shear-dominated cratering flow [Collins et al., 2004]. The fracturing and comminution processes create pore space between fragments and within fractures, reducing the bulk density of the subcrater material. The tendency of rocks and granular materials to change (typically increase) volume during shear failure is known as dilatancy, after Reynolds [1885] who first introduced the term to describe shearing granular materials. Calculation of fracturing-i.e., damage accumulation-is routine in modern numerical impact simulations [e.g., Ivanov et al., 1997;Collins et al., 2004]; accounting for dilatancy is not. As a result, most impact simulations do not predict correctly density changes beneath an impact crater, which has limited the scope for comparison of model results with geophysical data.

Dilatancy Model
Dilatancy is a well-known property of many geologic materials and is often quantified from triaxial compression experiments using the observed postfailure ratio of volumetric strain to axial strain (−d v ∕d 1 ) [Vermeer and De Borst, 1984]. Nonporous rocks and dense sands typically exhibit a maximum value of this ratio of order 1; the ratio is much lower (and may be negative) for loose sands, highly porous materials, clay-rich soils, and poor-quality rock masses [Vermeer and De Borst, 1984;Bolton, 1986;Hoek and Brown, 1997].
Dilatancy is also often quantified by a dilatancy angle that is defined as wherė( ⋅) is the appropriate plastic strain rate [e.g., Vermeer and De Borst, 1984]. This particular definition (using the sine function) is most appropriate for analysis of triaxial compression tests, where d v ∕d 1 is the quantity measured, under the assumption of a Mohr-Coulomb strength model, where shear stress m =( 1 − 3 )∕2 is linearly proportional to the mean stress m = ( 1 + 3 )∕2. As the plasticity model in iSALE is defined in terms of different measures of stress and strain, dilatancy will be quantified here by a "dilatancy coefficient" , given as where p is the accumulated plastic shear strain given by p = ∫ √ 2̇p ij̇p ij dt anḋp ij is the deviatoric plastic strain rate tensor. For triaxial deformation (̇2 =̇3), and are related by for simple shear deformation (̇2 = 0), ≈ sin . can also be expressed in terms of An alternative way to conceptualize dilatancy, described in detail by Reiner [1945], is to consider the change in internal pressure in a closely packed granular material undergoing isochoric deformation. If the deforming volume of granular material is not allowed to dilate to occupy a greater volume, then deforming grains will exert a pressure on the boundary of the volume in an effort to move over one another: permanent shear strain results in an increase in pressure. It is this concept that will be exploited here to account for dilatancy in the iSALE shock physics code.
A simple approach to account for dilatancy during shear failure in impact simulations is to supplement the pressure computed by the equation of state with a "dilatancy pressure, " representing the outward force of grains moving past one another, in cells where shear failure has occurred [Johnson and Holmquist, 1994]. This additional pressure effectively shifts the pressure-density relationship for the dilatant material up (to a higher pressure) so that if the material unloads to the reference pressure the density drops to a (dilated) bulk density that is below the reference density of the pristine material. Johnson and Holmquist [1994] defined the dilatancy pressure from energy considerations, such that a user-specified fraction of the elastic strain energy released during permanent deformation is converted into potential hydrostatic energy COLLINS ©2014. The Authors.  (1)) as a function of relative distension for a variety of sands, as measured in triaxial compression tests in which the mean effective pressure was 0.15-0.6 MPa. Data are taken from Bolton [1986], and references therein. Relative distension is a measure of the state of compaction of a granular material. It is the normalized distension between min , the distension at maximum closest packing of the granular material, and max , the maximum distension, which is often defined by the porosity that results from quickly inverting a container of the granular material. It is equivalent to the more commonly defined relative density. Note that the trend is linear and the dilatancy angle is zero at a relative density of approximately 0.25.
through an increase in pressure. O' Keefe et al. [2001] adopted a similar approach, but instead modified the equation-of-state surface at low pressure such that the reference density of the material decreased linearly as a function of plastic strain from 2800 kg/m 3 at the point of failure to a minimum reference density of 2000 kg/m 3 at a plastic strain of 0.1.
In the approach proposed here the supplemental dilatancy pressure results from an addition of distension (related to porosity by = [1 − ] −1 ) to cells undergoing shear failure (plastic shear deformation). Distension is the ratio of grain density to bulk density ( = s ∕ ); hence, if the bulk density remains constant, an increase in distension produces an increase in the grain density s and an associated increase in pressure. This approach is favored for three reasons: First, it is computationally expedient as it does not require a new field (dilatancy pressure) to be advected and stored. Second, in anticipation of later application to impact crater formation, this approach allows density changes owing to dilatancy to be separated from density changes owing to heating, as the former is recorded as a change in porosity whereas the latter is recorded as a change in material density. Third, it is readily combined with the existing algorithm to account for compaction of pore space in iSALE (the -model) . The recent discovery that the Moon's crust has substantial (∼10%) porosity suggests that both dilatancy and compaction were important processes in the evolution of the lunar crust, which provides strong motivation for an approach that can account for both processes at the same time.
The equation used to update the distension in cells undergoing shear deformation is The first term in this equation, , relates the dilation to a change in distension. It can be defined by considering the idealized case of a dilatant material that experiences simple shear deformation at constant internal pressure; that is, for which the grain density s remains constant. From the definitions of distension and volume strain, v = ln( 0 ∕ ), where 0 is the initial bulk density, it follows that d The second and most important term in equation (5) is the dilatancy coefficient, = d v d p , as described previously. A semiempirical approach is adopted to define this term, which rock and soil mechanics experiments indicate is dependent on porosity, strain, and confining pressure [Vermeer and De Borst, 1984;Bolton, 1986;Hoek and Brown, 1997]. Figures 1 and 2 show maximum dilatancy angle (defined using equation (1)) as a function of relative distension and confining pressure, respectively, from a range of sources [Bolton, 1986;Alejano and Alonso, 2005, and references therein]. Trends derived from these empirical data form the basis of the dilatancy model described here: is defined as a function of distension (porosity), pressure p, and temperature T according to   (1)) as a function of pressure for a variety of sands and rocks as measured in triaxial compression tests. Data for sands are taken from Bolton [1986], and references therein; data for rocks are taken from Alejano and Alonso [2005], and references therein. The abscissa represents the mean effective stress at the point of failure p for the sand data and the confining pressure 3 for the rock data. The data are well fit by a logarithmic function.
where max , c , min , p lim , and are constant model parameters described in detail below and T m is the melt temperature of the material at pressure p.
Hence, the dilatancy coefficient has a maximum value ( max ) at zero porosity ( = 1), pressure, and temperature and decreases as any of these three variables increase. Equation (6) is valid for T ≤ T m ; min < < max ; and 10 5 Pa < p ≤ p lim . The dilatancy coefficient tends to zero as (i) the distension approaches a critical distension c ; (ii) the pressure approaches the limiting pressure p lim ; or (iii) the temperature reaches the melt temperature T m . Note that T m is the melt temperature at pressure p and not a constant. The form of the distension term is supported by dilatancy measurements for sand (see Figure 1) [Bolton, 1986, and references therein] and permits negative dilatancy coefficients for c < < max . The pressure term is supported by dilatancy measurements for rocks and sand (see Figure 2) [Bolton, 1986, and references therein]. These functional forms are similar to the model of dilatancy proposed by Bolton [1986]. The efficacy of both these terms was also tested by comparing numerical triaxial compression simulations against experimental data (see Appendix A). The form of the temperature term has no direct empirical basis; it was chosen for consistency with the thermal softening function in the strength model used in iSALE [Ohnaka, 1995;Collins et al., 2004]. Equation (6) contains four material constants ( max , c , p lim , and ). The last term, which is also a constant in the strength model, can be found from experimental measurements of strength degradation as a function of temperature [e.g., Ohnaka, 1995], which suggest ≈ 0.7-1.2 [Collins et al., 2004;Ivanov et al., 2010]. Measurements of dilatancy in sand suggest that c ≈ 0.25 min + 0.75 max (see Figure 1), where max (typically 1.8-2.1) is the distension of the sand in its most distended state and min (typically 1.45-1.6) is the distension of the sand under maximum closest packing (cf. maximum and minimum voids ratio). Typical distension values for a variety of sands imply c ≈ 1.7-2. However, if distension is redefined relative to min (i.e., the sand is considered to have a distension of one at maximum closest packing), the same data imply c ≈ 1.15-1.25, which corresponds to a density decrease of 15-20%. Although no empirical support for a value of c in this range for nonporous rock is available, validation simulations conducted in this work are consistent with c ≈ 1.2. We also note that bulking observed in the ejecta and the breccia lens at terrestrial craters is typically 15-20% [e.g., Pilkington and Grieve, 1992], consistent with this value representing the minimum bulk density caused by impact-related dilatancy.
Measurements of dilatancy suggest that the confining pressure 3 required to suppress dilatancy is ∼10 MPa for sand and ∼100 MPa for a range of sedimentary rocks (see Figure 2); the mean stress (pressure) associated with this suppression will be somewhat higher. Simulations conducted to validate the dilatancy model as part of this work suggest p lim ≈ 50 MPa for sand and p lim ≈ 800 MPa for rock (see Appendix A). It is also noteworthy that the range of confining pressures p lim = 100-800 MPa is consistent with the depths at which fractures in the crust are expected to close under lithostatic pressure [e.g., Pilkington and Grieve, 1992].
An estimate of the maximum dilatancy coefficient max as defined here can be derived from measurements of the maximum dilatancy angle obtained from shearbox and triaxial compression tests (see equation (3)). The estimate is not perfect, however, as pressure is never zero in such experiments. Vermeer and De Borst [1984] suggest that typical dilatancy angles for a range of geologic materials including sand, concrete, and rocks are between 0 and 20 • . Alejano and Alonso [2005] show peak dilatancy angles at the onset of failure COLLINS ©2014. The Authors. as large as 50 • for a range of sedimentary rocks under low-confining pressures (e.g., see Figure 2); however, in such cases the dilatancy angle drops rapidly to below 20 • after less than 1% strain. Hence, a realistic expected range of max ≈ 0 − 0.5, although values of order 1 may be possible for some rocks. This range is consistent with results of simulations of trixial compression tests described in Appendix A.
It is worth noting that when applying the present dilatancy model to the deformation of large rock masses, it may not be appropriate to derive model parameters from measurements of centimeter-scale samples. Although triaxial tests exploring the sensitivity of peak dilatancy angle to sample size suggest only a modest (∼10%) reduction in peak dilatancy angle over an order of magnitude range in sample diameter [Medhurst, 1997], if extrapolated to the kilometer scale, this would result in a factor of 2 reduction in dilatancy angle. Moreover, experience from engineering assessments of rock mass strength and failure suggests that dilatancy angle varies between 0 and 12 • and is a function of rock mass quality [Hoek and Brown, 1997]. A dilatancy angle one quarter of the friction angle (∼ 12 • ) is expected for a very good quality hard rock mass with a Geological Strength Index (GSI) >75; an angle of 4 • (1/8 of the friction angle) is expected for an average quality hard rock mass with a GSI ≈50; and a zero degree dilatancy angle is expected for a very poor quality rock mass with a GSI <30 [Hoek and Brown, 1997]. As will be demonstrated subsequently, impact simulations that assume a maximum dilatancy angle of a few degrees generate density anomalies that are most consistent with gravity data collected at terrestrial impact structures, suggesting that these representative values for rock masses are more appropriate for planetary cratering studies than those derived from laboratory measurements.

Cratering Simulations
The dilatancy model described in section 2 was implemented in the iSALE shock physics code Collins et al., 2004], a multirheology, multimaterial extension to the finite-difference SALE hydrocode [Amsden et al., 1980]. iSALE is well tested against laboratory experiments at low-and high-strain rates , as well as other impact simulation codes [Pierazzo et al., 2008], and has been widely used to numerically simulate terrestrial impact crater formation [e.g., Collins and Wünnemann, 2005;Collins et al., 2008;Wünnemann et al., 2005;Goldin et al., 2006]. Dilatancy model implementation involved modifying the "radial return" plasticity algorithm [Wilkins, 1964], which defines how deviatoric stresses that exceed the yield envelope are mapped back to the yield surface during plastic deformation. A description of these modifications is given in supporting information, and model verification and testing are described in Appendix A.
To explore the role of dilatancy in impact crater formation, iSALE was used to simulate impacts on Earth at a range of scales. Simulations considered a range of impactor diameters between 100 m and 20 km, which produced a range in crater size that spans the known terrestrial impact crater record (for those impacts COLLINS ©2014. The Authors. Figure 3. Impact-induced dilatancy during simple crater formation. Porosity evolution from a simulation of a 250 m diameter asteroid impact at 15 km/s into a crystalline target with a maximum dilatancy coefficient max = 0.045, which results in a ∼4 km diameter simple crater. Thin solid lines show porosity contours of 0.1, 1, 5, 10, and 15%; the dotted lines indicate target rock deformation. not significantly affected by atmospheric entry). A spatial resolution of 10-20 cells per projectile radius was used for all simulations. For simplicity and computational expediency, all impacts were treated as vertical and assumed a constant impact velocity of 15 km/s. The influence of impact angle and velocity on impact-induced dilatancy will be explored in future work. As the principal aim of this work was to explore trends in impact-related dilatancy with crater size, rather than replicate specific impact events, most simulations used a simple, uniform, nonporous crystalline rock material model (granite) to represent the target and a constant gravitational acceleration g = 9.81 m/s 2 . The exceptions were the two largest impact scenarios where a mantle layer was included below 35 km depth. The granite material model was also used to represent the impactor. To describe the thermodynamic behavior of the material, equation-of-state (EoS) tables generated using the analytic EoS (ANEOS [Thompson and Lauson, 1972]) software package were used, with input parameters appropriate for the granite crust [Pierazzo et al., 1997] and dunite mantle [Benz et al., 1989].
To describe the resistance of the material to shear deformation, the strength and damage models described by Collins et al. [2004] and Ivanov et al. [2010] were used, in conjunction with the modifications for dilatancy. Note that the dilatancy modifications did not include changes to Poisson's ratio , which is assumed to be constant (see Table 1), implying a fixed ratio between the shear modulus G and the bulk modulus K (determined by the equation of state). As Poisson's ratio of most rocks depends on porosity [Walsh, 1965], and can increase or decrease with porosity depending on the crack geometry [Shearer, 1988], it may be important to account for this variation in future refinements of the dilatancy model (and the related porous compaction model). For now, the assumption of constant Poisson's ratio, which is widely employed in numerical impact models, serves as a useful first approximation. The strength model modifications for dilatancy include a new dependence of the coefficient of friction (for the damaged material) on distension km diameter simple craters and the resultant Bouguer gravity anomaly at the level of the preimpact surface. Simulation results are shown for three different values of the maximum dilatancy coefficient max = 0.045, 0.09 and 0.18, which are appropriate for low-GSI, average-GSI, and high-GSI rock masses, respectively. Thin solid lines represent porosity contours of 0.1, 1, 5, 10, and 15%; the dotted lines indicate target rock deformation.
(supporting information), which was assumed to drop linearly from d at = 1 to c at c ( c is the coefficient of friction in the critically dilatant state and was taken to be 0.6 d based on the behavior of dense sand [Hettler and Vardoulakis, 1984]. While this modification has little effect on complex crater collapse, it was found to be important to replicate the observed thickness of breccia in terrestrial simple craters. For large impact crater formation the strength model must include some form of transient target weakening model that facilitates deep-seated gravitational collapse of the transient crater Melosh and Ivanov, 1999]. The physical rationale for this apparent target weakening is still a matter of debate [e.g., Senft and Stewart, 2009]; in the models presented here it is assumed that acoustic fluidization of the target rocks is the primary weakening mechanism [Melosh, 1979]. In the present version of iSALE, the effects of acoustic fluidization are incorporated using the "block (oscillation) model" [Melosh and Ivanov, 1999;Ivanov and Artemieva, 2002;Wünnemann and Ivanov, 2003]. The choice of block model and other strength model input parameters was based on previous successful models of large crater formation on Earth [e.g., Collins and Wünnemann, 2005;Collins et al., 2008;Wünnemann et al., 2005;Goldin et al., 2006] and the Moon [Wünnemann and Ivanov, 2003]; however, it is noted that the choice of block-model parameters is nonunique and the present model parameters are not regarded as definitive. All the important strength model parameters used in the simulations presented here are included in Table 1; the interested reader is referred to Collins et al. [2004] and Ivanov et al. [2010] for more detailed parameter definitions.
To explore the influence of dilatancy on impact-generated porosity, simulations were performed with three different values of the most important model parameter, the maximum dilatancy coefficient max = 0.045, 0.09, and 0.18. These represent approximate values appropriate for rock masses with a low, average, and high Geological Strength Index (GSI), respectively, following guidelines relating GSI to the ratio of dilatancy angle and friction angle [Hoek and Brown, 1997]. These and the remaining dilatancy model parameters are presented in Table 1. COLLINS ©2014. The Authors.

Figure 5.
Porosity as a function of depth beneath the preimpact target surface as derived from numerical simulations of a ∼4 km diameter crater on Earth. Simulation results are shown for three different values of the maximum dilatancy coefficient max = 0.045, 0.09, and 0.18, which are appropriate for low-GSI, average-GSI, and high-GSI rock masses, respectively. Shown for comparison is the inferred porosity beneath the floor of the Brent impact crater, Canada, based on density measurements from drill-core sampling [Innes, 1961]. Solid symbols show the mean porosity over each sampling interval; the boxes delimit the full range of porosity (minimum to maximum) over the interval. A depth of 200 m was added to the sampling depths to correct for postimpact erosion.

Impact-Induced Dilatancy in Simple Crater Formation
Impact simulations that used an impactor diameter L of 100-250 m produced simple craters. The crater evolution and generation of porosity by dilatancy for the L = 250 m impact is depicted in Figure 3. In this simulation max = 0.045, which is the assumed value for a low-GSI rock mass; the effect of maximum dilatancy coefficient on the porosity distribution beneath the final crater is shown in Figure 4. As illustrated by the 0.1% porosity contour at 1 s after impact, propagation of the shockwave through the target results in only a minor porosity increase < 1% that extends to a radius of ∼1 km, apart from very near the surface where dilation is enhanced by the interaction between the detached shock and the free surface (Figure 3). Although the plastic strain rate is high during shock propagation, the elevated pressure in the shock wave suppresses dilatancy away from the surface. Dilation of the rocks surrounding the crater increases during the subsequent excavation flow, particularly at shallow depths and in the ejecta curtain, where confining pressure is low. This results in a porosity of a few percent in the wall and floor of the transient crater, which forms at around 10 s. Finally, further dilation occurs during debris sliding as the unstable transient crater rim collapses back into the crater. The final porosity of the breccia lens ranges from 5% at its base to the maximum porosity of 17% defined by the critical dilatancy distension c = 1.2. Simple crater formation in a more strongly dilatant target ( max = 0.09-0.18) results in higher-absolute porosities, but the spatial distribution of porosity is qualitatively similar (Figure 4).
The final porosity distribution beneath the crater can be compared with observations and geophysical models of terrestrial simple impact craters. Figure 5 shows porosity as a function of depth (relative to the preimpact surface) for three impact simulations of a L = 250 m impactor that each assume a different maximum dilatancy coefficient max corresponding to nominal values for a low-GSI rock mass, an average-GSI rock mass, and a high-GSI rock mass [Hoek and Brown, 1997]. For comparison, porosity as a function of depth beneath the Brent impact crater was inferred from density measurements of drill-core samples [Innes, 1961]. The well-studied Brent impact crater, Canada, is the type-example of a simple crater formed in a crystalline target [Dence et al., 1977;Grieve et al., 1989]. It has an apparent crater diameter of 3 km; the estimated rim diameter of the pristine crater is 3.8 km [Dence et al., 1977]. The present level of erosion is estimated to be about 200 m below the preimpact target surface, the apparent crater floor is 260 m below the surface, and the crater is filled with postimpact sediments [Grieve et al., 1989]. Density measurements (mean, min, max) from ∼30 m intervals of a ∼1 km drill core near the crater center [Innes, 1961] were converted to porosity assuming a nominal constant grain density of 2.7 g/cc. To adjust for postimpact erosion, sampling depths were shifted down by 200 m. The porosity-depth profiles from all three impact simulations are broadly consistent with the observed decrease in porosity with depth beneath the crater, particularly at the base of the breccia lens (z ≈ 1 km) and below. Within the breccia lens the model that assumed the lowest maximum dilatancy angle (low GSI) is most consistent with observations. More strongly dilatant behavior is not supported by this comparison but might be accommodated by postimpact closure or filling of pore space, in particular by water or other pore fluids not accounted for in the model. Interestingly, the numerical simulations predict enhanced dilatancy at the base of the breccia lens caused by shear localization; this is  [Innes, 1961;Grieve et al., 1989]. For direct comparison with the Brent gravity anomaly, the Bouguer gravity anomalies for the simulated craters were calculated at the level z = −0.2 km (mass differences above this level were neglected), to account for postimpact erosion, and assuming a relative density of −0.17 g/cc for the sedimentary crater fill. Simulation results are shown for three different values of the maximum dilatancy coefficient max = 0.045, 0.09, and 0.18, which are appropriate for low-GSI, average-GSI, and high-GSI rock masses, respectively. consistent with the observed local maximum in porosity at 1 km depth. It is possible that higher-resolution simulations would better capture shear localization at the base of the breccia lens, resulting in a more prominent porosity high that is a better match to the observations.
An alternative method to test the veracity of the numerical simulation results is to compare the Bouguer gravity anomaly of the final simulated impact structure with Bouguer anomalies over terrestrial craters. Figure 4 shows the relationship between the simulated porosity anomaly beneath a ∼4 km diameter simple crater (for L = 250 m; low GSI, max = 0.045; Figure 3) and the resulting Bouguer gravity anomaly. The method for calculating the Bouguer anomaly from the porosity distribution is described in the supporting information: the anomaly is reduced to the preimpact surface (i.e., the mass distribution above the preimpact surface is removed, and the crater is filled to the preimpact surface with material with a density equivalent to the preimpact target). The magnitude of the Bouguer gravity anomaly is −3.6 mgals; the equivalent anomaly magnitudes for impacts on the average-GSI and high-GSI targets are −5.1 mgals and −6.8 mgals, respectively. This range is within the spread of negative Bouguer anomaly magnitudes for large terrestrial simple craters compiled by Pilkington and Grieve [1992]; for example, Brent has a Bouguer gravity anomaly of −5 mgals. However, for many of these craters the Bouguer gravity anomaly is affected by postimpact erosion and/or infill by low-density sediments. To more accurately compare numerical model predictions with the observed gravity anomaly at Brent, the Bouguer gravity anomalies were recalculated to account for postimpact erosion and sedimentary infill, by reduction to a datum 200 m below the preimpact surface and by assuming a relative density of −170 kg/m 3 for the crater fill, consistent with drill-core measurements [Innes, 1961]. Figure 6 compares the observed Bouguer gravity anomaly at Brent to the anomalies produced by impact simulations that assumed a low-, average-and high-GSI target responses. While the diameter of the anomaly is consistent with all three numerical results, the −5 mgal Bouguer anomaly magnitude is most consistent with impact on a weakly dilatant, low-GSI target. However, postimpact filling of pore space by water or compaction by burial are not accounted for in the model and would reduce the synthetic gravity anomaly magnitude.
Although the maximum dilatancy angle has a considerable effect on the gravity anomaly magnitude, it has only a minor influence on transient crater dimensions and final crater topography ( Figure 6). Transient (precollapse) crater volume (below the preimpact surface) is less than 3% larger for max = 0.045 (low GSI) compared with max = 0.18 (high GSI), owing to a small (∼1%) decrease in transient crater diameter with max . The effect of dilatancy on final crater dimensions is more pronounced, emphasizing the importance of the collapse stage in porosity generation. Final crater depth is 20% larger and final diameter is 5% smaller for max = 0.045 compared with max = 0.18 (Figure 6).

Impact-Induced Dilatancy in Complex Crater Formation
Impact simulations that used an impactor diameter L ≥ 500 m produced complex craters. As shown in several previous numerical studies [e.g., Melosh and Ivanov, 1999;Wünnemann and Ivanov, 2003], dynamic weakening of the target material by acoustic fluidization facilitates deep-seated uplift of the crater floor and collapse of the rim, resulting in a crater evolution qualitatively different to simple crater formation. With increasing impactor diameter the extent and longevity of weakening is increased, resulting in more extensive collapse. In small complex craters central uplift produces a central peak; for larger impactor sizes the central uplift overshoots the preimpact surface and then collapses downward and outward to form an annular topographic high (i.e., peak ring) [Collins et al., 2002;Ivanov and Artemieva, 2002]. This size-morphology progression is well documented in previous numerical impact modeling studies [e.g., O' Keefe and Ahrens, 1999;Melosh and Ivanov, 1999;Wünnemann and Ivanov, 2003]. The new information afforded by the present study is the predicted generation of porosity by dilatancy during complex crater formation. Figure 7 depicts the generation of porosity during a simulation of a L = 3.2 km impact, which results in a crater approximately 40 km in diameter (at an impact velocity of 15 km/s and for max = 0.045). The effect of maximum dilatancy coefficient on the final porosity distribution and Bouguer gravity anomaly is shown in Figure 8.
As for simple crater formation, porosity created during the passage of the shock wave is limited to <<1% because of the high pressure in the shock. In addition, the higher overburden pressure suppresses shock-induced dilatancy below a few km depth, as shown by the extent of the 0.1% porosity contour at T = 5s (Figure 7). Dilation of the rocks surrounding the crater increases during the subsequent excavation flow, particularly at shallow depths and in the ejecta curtain, where confining pressure is low, but high pressure limits dilation beneath the crater floor. This results in a porosity of a few percent in the rim of the transient crater and <1% beneath the crater floor, which forms at around 40 s. Further dilation occurs during COLLINS ©2014. The Authors. Figure 8. Porosity distribution beneath the final simulated ∼40 km diameter complex craters and the resultant Bouguer gravity anomaly at the level of the preimpact surface. Simulation results are shown for three different values of the maximum dilatancy coefficient max = 0.045, 0.09, and 0.18, which are appropriate for low-GSI, average-GSI, and high-GSI rock masses, respectively. Thin solid lines represent porosity contours of 0.1, 1, 5, 10, and 15%; the dotted lines indicate target rock deformation. transient crater collapse, particularly in the subsiding shallow, near-rim area. More modest dilation occurs in the central uplift, again owing to the high pressures associated with the convergent flow-only during outward (divergent) collapse of the central uplift (T > 100 s) are conditions in the central region conducive to dilation. The consequence of dilatancy suppression in the central uplift is that at a given depth, porosity first increases with radial distance until it reaches a maximum beneath the annular trough, before decreasing again at greater distances. The higher overburden pressure (compared to the simple crater scenario) also implies that the maximum porosity is <15%; the critical dilatant state is not achieved in this simulation, although it is possible that higher-resolution simulations might produce higher porosity in the highly strained near-surface region.
The attenuation of porosity with depth beneath the final craters in Figure 8 can be compared with results from the deep drill core beneath the center of the ∼40 km diameter Puchezh-Katunkski impact structure, Russia [Masaitis and Pevzner, 1999]. It is worth noting that the assumption of a uniform crystalline target is not a good approximation for the Puchezh-Katunkski preimpact target, which comprised a thick sedimentary sequence above a heterogeneous crystalline basement. However, as the sedimentary layer would have been almost entirely excavated from the center of the crater, it is reasonable to compare the present model results with data from the drill core through the central uplift of crystalline basement. At Puchezh-Katunkski, porosity measurements of crystalline rock samples to depths exceeding 5 km [Masaitis and Pevzner, 1999] show a decrease in mean porosity from ∼12% at the top of the impact-affected sequence to 1-2% at 5 km depth. Figure 9 compares the range of measured porosities (mean, standard deviation, minimum, and maximum) to porosity-depth profiles beneath the center of the final crater for simulations that assumed three different maximum dilatancy coefficients, appropriate for low-, average-, and high-GSI rock masses. The model results are broadly consistent with observations for the range of dilatancy model parameters used. At shallow depths the model results for the strongly dilatant, high-GSI target model are most consistent with the measured porosities; nearer the base of the drill core, the more weakly dilatant target models (low-GSI and average-GSI) give results that are more Figure 9. Porosity as a function of depth beneath the preimpact target surface as derived from numerical simulations of a ∼40 km diameter crater on Earth (Figure 7). Simulation results are shown for three different values of the maximum dilatancy coefficient max = 0.045, 0.09, and 0.18, which are appropriate for low-GSI, average-GSI, and high-GSI rock masses, respectively. Shown for comparison are the ranges of observed porosity beneath the Puchezh-Katunkski impact crater in Russia, based on drill-core sampling [Masaitis and Pevzner, 1999]. The grey boxes span 1 standard deviation about the mean sample porosity; the white boxes span the minimum-maximum range of porosity. consistent with observations. A potential source of discrepancy between the models and observations at Puchezh-Katunkski and Brent is that the inferred porosity of the core samples only accounts for sample-scale porosity, such as small fractures and intergranular pore space, and does not include macroscale porosity, such as large open fractures or interblock shear zones. The presence of large-scale porosity at the boundaries between displaced megablocks beneath the floor of the Puchezh-Katunkski crater has been inferred based on core output from the deep drilling [Ivanov et al., 1996] but is not accounted for in Figure 9. Hence, the inferred porosity from drill-core samples may underestimate the total subcrater porosity.
Also shown in Figure 9 is the porosity-depth profile in the average-GSI model at a radius of 10 km (dot-dashed line), through the annular trough where dilatancy is most effective (see Figure 7). Porosities beneath the annular trough are slightly above the maximum porosities observed in the deep drill core at the center. A comparison between this profile and the profile beneath the center of the crater in the same model (dashed line) is a measure of the effect of the convergent flow in suppressing dilatancy in the central uplift.
The Bouguer gravity anomaly of the final simulated impact structure shown in Figure 7 is depicted in Figure 8, together with the corresponding anomalies for simulations that assumed max = 0.09 and 0.18. As with the simple crater scenario, the magnitude of the porosity and Bouguer gravity anomalies increases with maximum dilatancy coefficient, but the qualitative distribution of porosity and associated character of the Bouguer anomaly is the same. Figure 10 shows the final porosity distribution beneath other simulated complex craters and their associated Bouguer gravity anomalies. Results are shown for impactor diameters of L = 1 km, 2 km, 5 km, and 10 km, which produced craters with rim diameters D = 14 km, 24 km, 60 km, and 110 km, respectively. The simplified target representation (uniform crystalline rock) precludes a detailed comparison of numerical model results with observed anomalies at specific craters. However, the change in anomaly character and magnitude with crater size can be compared with trends in terrestrial observations [Pilkington and Grieve, 1992]. With increasing crater diameter the shape of the Bouguer anomaly changes from a broad gravity low with no central anomaly (D < 24 km) to a broad low with a small, central (relative) positive anomaly at D = 40 km. This is broadly consistent with the gravity signature of terrestrial craters: most terrestrial craters larger than 30 km diameter exhibit a central gravity high, whereas most smaller craters do not. The amplitude of this relative central high in the Bouguer gravity increases with increasing crater diameter; it is almost equivalent in magnitude to the broad negative anomaly at a crater diameter of D = 110 km. In the numerical simulations, which almost all assume an approximately uniform density preimpact target, the cause of this central positive anomaly is the suppression of dilatancy in the central uplift: less porosity is created in the central uplift than the annular trough, owing to elevated pressure in the convergent uplift and higher temperatures beneath the crater floor. Only in the 110 km diameter crater simulation does uplift of high-density rocks (in this case the mantle) contributes to the gravity anomaly, and even in this case, the contribution is very small.
For each simulated impact crater, the magnitude of the broad negative Bouguer gravity anomaly Δg is shown as a function of crater diameter D in Figure 11. Results are shown for the three suites of simulations using different values for the maximum dilatancy coefficient, appropriate for low-, average-, and high-GSI rock masses. Results from all three simulation suites show a consistent trend. For craters smaller COLLINS ©2014. The Authors. Figure 10. Porosity distribution beneath the final simulated complex craters formed by impactor diameters of (a) L = 1 km, (b) L = 2 km, (c) L = 5 km, (d) L = 10 km, and the resultant Bouguer gravity anomaly at the level of the preimpact surface. All simulations used a low-GSI rock mass maximum dilatancy coefficient max = 0.045. Thin solid lines represent porosity contours of 0.1, 1, 5, 10, and 15 %; the dotted lines indicate target rock deformation. than D ≈ 24 km, negative gravity anomaly magnitude increases approximately linearly with crater diameter; for larger craters, anomaly magnitude increases with diameter at a smaller rate. Intuitively, gravity anomaly magnitudes are consistently larger (more negative) for models that assumed a more strongly dilatant material response; negative anomalies are 2-3 times larger for the high-GSI target models than for the low-GSI target models.
The model results can be compared with the observed magnitude of the negative anomaly at terrestrial impact structures in crystalline and mixed crystalline/sedimentary targets compiled by Pilkington and Grieve [1992]. There is considerable scatter in the terrestrial gravity anomaly data, particularly at low crater diameters. This is in part because many gravity anomalies include the contribution from low-density postimpact crater fill, which can be the dominant cause of the anomaly in simple craters. For example, the large negative anomalies at Lonar (−3.6 mgals), Pretoria Salt Pan (−5 mgals), Tenoumer (−10 mgals), and Roter Kamm (−9.3 mgals) are all dominated by sedimentary infill (grey symbols in Figure 11). The residual contribution Figure 11. Negative Bouguer gravity anomaly magnitude as a function of crater diameter for various simulation suites. Shown for comparison are Bouguer gravity anomaly magnitudes for terrestrial craters in crystalline and mixed sedimentary and crystalline targets after Pilkington and Grieve [1992]. The line Δg = D is shown to illustrate proportionality between gravity anomaly magnitude and crater diameter. The grey symbols denote the gravity anomalies discussed in the text that include the effect of postimpact sedimentary infill; the corrected values are shown in white.
Accounting for these corrections (which are not exhaustive), several aspects of the terrestrial data are consistent with the model results. The spread in gravity anomaly magnitude at a fixed crater size is approximately a factor of 2-3; at small crater diameters the negative anomaly magnitude is approximately proportional to crater diameter; and there is an apparent reduction in the slope of log(Δg)-log(D) at a crater diameter of 10-20 km diameter. At small crater sizes D < 24 km, simulated anomaly magnitudes produced in low-to high-GSI targets span the range of observed negative gravity anomaly magnitudes, although a greater number of observed anomalies are consistent with simulated impacts in low-GSI targets. At larger crater sizes, only the gravity anomalies produced by impacts in low-GSI targets are consistent with observed anomaly magnitudes.

Discussion
Numerical simulations of terrestrial impacts in uniform crystalline targets, using the dilatancy model with parameters appropriate for rock masses of varying degrees of quality, result in porosity and gravity anomalies that are broadly consistent with observation. The impact simulations that assumed a low maximum dilatancy coefficient, max = 0.045, appropriate for low-quality, low-GSI rock masses, appear to generate porosity distributions and gravity anomalies most consistent with observation. More strongly dilatant behavior (i.e., 0.045 < max < 0.18) cannot be ruled out (and indeed may be appropriate in some cases) but would require some of the impact-induced porosity to be closed or filled after the impact to be generally applicable. At least partial pore filling is expected in environments where liquid water is stable, which would reduce the density and gravity anomaly of the crater. A lack of correlation between impact crater age and gravity anomaly magnitude [Pilkington and Grieve, 1992] suggests that postimpact pore closure is not universal; however, a more detailed comparison between model predictions of postimpact porosity and geophysical observations is required to better constrain the dilatancy model parameters.
Another factor that will influence gravity anomaly magnitude is the initial porosity of the target. In the simulations presented here, for simplicity the preimpact target was assumed to be nonporous. The presence of initial pore space will have two effects. First, the impact will permanently compact pore space in a proximal zone where pressures exceed the crushing strength , increasing the relative density of the target near the crater center. Second, as the dilatancy coefficient is a decreasing function of porosity, target rocks with a high initial porosity will be more weakly dilatant than those with little-to-no porosity. Hence, it is likely that initial porosity will act to reduce the overall mass deficit caused by impact-induced damage and may even result in a relative mass excess if pore space compaction dominates pore space creation by dilatancy. The role of preimpact porosity is of particular importance for understanding impact-induced density changes on the Moon, where the present average crustal porosity is ∼12% [Wieczorek et al., 2012].
The decrease in dilatancy coefficient with increasing pressure is an important property of the model and has three effects on impact-induced dilatancy. First, while a small amount of dilation occurs during shock wave propagation, in general the high pressures in the shock wave suppress the generation of porosity as it propagates through the target rocks. The exception is the shallow near-surface interference zone where shock and release wave interactions shield the material from high pressure. Second, at depths exceeding about 10 km on Earth the confining pressure is sufficient to suppress porosity generation at any stage during crater formation. As a result, the majority of the impact-generated pore space is created by shear deformation near the surface and late in the crater formation process-during excavation and collapse. It is this effect that results in the change in the log slope of negative gravity anomaly magnitude versus crater diameter ( Figure 11). Third, in complex craters, where deep rocks are raised to the surface in a central uplift, dilatancy is also suppressed in the central region because of the high pressure in the convergent flow. As a consequence, the central gravity high in many terrestrial complex craters may be caused, at least in part, by suppression of dilatancy, rather than (or in addition to) uplift of dense rocks [Pilkington and Grieve, 1992].
The effect of pressure on dilatancy also suggests that impact-induced dilatancy may be more effective on planetary surfaces with lower gravity. To test this hypothesis, a numerical impact simulation of a 5 km diameter impactor striking the lunar surface at 15 km/s was performed for comparison with the terrestrial impact simulations. For simplicity and to aid direct comparison, the lunar impact simulation employed the same material models for the impactor, crust, and mantle as the terrestrial simulations (with max = 0.09); only the crustal thickness (40 km) and gravitational acceleration (1.62 m/s 2 ) were modified to represent the Moon. The porosity distribution beneath the final ∼70 km diameter crater is compared with the corresponding terrestrial impact simulation result in Figure 12.
Impact-induced dilatancy is dramatically more effective on the Moon than on Earth; equivalent porosity contours are 2-3 times deeper beneath the lunar crater compared to the terrestrial crater. At this approximate crater size, that negative Bouguer gravity anomaly magnitude is a factor of 2 larger on the Moon than on Earth (Figure 12). Moreover, while dilatancy is suppressed in the central uplifts in both craters, the character of the Bouguer anomaly over each crater is very different. The terrestrial crater gravity anomaly exhibits a broad gravity low with a relative high over the central region, whereas the lunar crater gravity anomaly is a bowl-shaped low. The efficiency of impact-induced dilatancy on the Moon supports the idea that the high-inferred porosity of the Moon's crust is the result of impact-induced fracturing and fragmentation [Wieczorek et al., 2012]. However, as the lunar impact simulation assumed a nonporous preimpact target, the resultant porosity distribution probably represents an upper limit to what can be produced by a single impact of this size. If impact fracturing progressively increased the porosity of the lunar crust over time, the effectiveness of impact-induced dilatancy is likely to have diminished concomitantly.
The numerical simulations of terrestrial crater formation presented here provide new insight into impact-induced density changes that accompany fracturing and fragmentation. However, several limitations of the simulations remain. Principal among these is the nonunique choice of acoustic fluidization model parameters ( and ) [Wünnemann and Ivanov, 2003] which dictate the extent of target weakening in the simulations and hence the resultant size-morphology progression. The parameters used in this work produce a range of final simulated craters that are broadly consistent with observational constraints, such as the simple-to-complex transition, particularly given the simple assumption of a uniform crystalline target. However, alternative parameter choices may produce a suite of simulated (different) craters that are more consistent with observational constraints than those presented here. In particular, despite dilatancy shallowing the craters, the simulations produce complex craters that are systematically deeper than is observed, which may suggest that the magnitude of target weakening in the present models is insufficient. However, tests using alternative acoustic fluidization model parameters resulted in very similar final porosity distributions for the same final crater size, suggesting that the porosity predictions presented here are robust. It is also noted that the effect of impact velocity, impact angle, initial target porosity, and the remaining dilatancy model parameters ( c and p lim ) on impact-induced dilatancy has not been investigated. Impact angle and initial porosity, in particular, are expected to influence the spatial distribution of impactinduced porosity.

Conclusions
A simple model to account for dilatancy during impact crater formation has been developed and implemented in the iSALE shock physics code. The model accounts for the effects of pressure, temperature, and distension on the tendency for geological materials to dilate, and has been tested and calibrated against triaxial experiments on rocks and dense sand. Numerical impact simulations of simple crater formation using the new model produce porosity distributions consistent with drill-core sampling, as well as geophysical modeling and observed gravity anomalies. Dilatancy model parameters appropriate for low-quality rock masses (i.e., with a low Geological Strength Index; GSI) produce a porosity distribution most consistent with observation; high-GSI behavior would require postimpact pore closure, healing or filling by water to explain observed anomaly magnitudes. While dilatancy has a profound influence on the postimpact density distribution beneath an impact structure, it has a relatively minor influence on topography, reducing crater depth by ∼20%.
The tendency for rock to dilate less when shearing under high pressure is an important property of the model. Pressure suppresses impact-induced dilatancy: (i) in the shock wave, (ii) at depth beneath the crater floor, and (iii) in the convergent flow that forms the central uplift. As a result, numerical impact simulations predict differences in subsurface porosity distribution between simple and complex craters. In simple craters, bulking in the debris lens is quite uniform and porosity decays exponentially both beneath the crater floor and radially from the crater rim; whereas, in complex craters, porosity is greatest in the annular trough and decays exponentially with depth and radial distance.
Simulation results predict that negative gravity anomaly magnitude is proportional to crater size for crater radii less than the depth at which dilatancy is suppressed, broadly consistent with observed gravity anomalies at terrestrial craters [Pilkington and Grieve, 1992]. At larger crater sizes, the increase in negative gravity anomaly magnitude with diameter is less pronounced and dilatancy suppression in the central uplift results in a relative central gravity high that increases in magnitude with increasing diameter.
The lower confining pressure at an equivalent depth on the Moon relative to Earth implies that impactinduced dilatancy is dramatically more effective on the Moon than on Earth for the same size impact in an initially nonporous target. This may be mitigated by the presence of substantial initial pore space in the lunar crust, as is the case for the modern Moon [Wieczorek et al., 2012].
Future application of the dilatancy model in more refined numerical simulations of specific craters that more faithfully represent preimpact target conditions will allow numerical predictions to be directly compared with geophysical observations, such as gravity and seismic velocity anomalies. This will provide much Figure A2. (top) Differential stress and (bottom) volumetric strain as a function of axial strain in a series of numerical triaxial stress tests for dense sand. Experimental data were sourced from Hettler and Vardoulakis [1984].
This description defines the critically dilatant state as = c and p∕K 0 = 1−1∕ c . Figure A1 demonstrates that the dilatancy algorithm in iSALE performs as expected under isochoric simple shear. Observe how distension and dilatancy pressure increase indefinitely under the assumption of constant dilatancy angle; whereas, definition of a critically dilatant state implies that distension and dilatancy pressure converge to the expected values with increasing strain. Note also that the relationship between dilatancy coefficient and distension in the present model implies that the critically dilatant state is only approached at an equivalent plastic shear strain ≈2. This contrasts with the model of O' Keefe et al. [2001], which assumed maximum dilatancy pressure was reached at a plastic strain of only 0.1.

A2. Triaxial Compression
Simulations replicating triaxial strength tests were performed using a single computational cell in two-dimensional, cylindrical geometry. Pressure boundary conditions were applied to the top and right boundaries to replicate conditions in typical triaxial strength tests: a constant radial confining pressure ( 3 ; right boundary) and increasing axial stress ( 1 ; top boundary). Free-slip boundary conditions were enforced on the left (symmetry axis) and bottom boundaries. The axial stress 1 was increased at a constant rate from an initial stress equivalent to the confining stress until an equivalent plastic strain of 2 was achieved. The initial state (density and internal energy) of the material in the cell was set consistent with the isotropic initial confining pressure. Suites of numerical experiments at various confining pressures were performed for both a rock-like material (marble) and dense sand; the influence of initial distension was also examined for sand. Elastic and strength model parameters for marble and dense sand were inferred from experimental results. A simple polynomial equation of state relating pressure to density (p = A + B 2 , where = ∕ 0 − 1.) was combined with a Drucker-Prager strength model (Y = min{ p, Y lim }). Dilatancy model parameters ( c , p lim , max ) were varied to achieve the best overall fit to the appropriate triaxial experiment data. Figure A2 shows an example of how the model fits experimental data for sand. Material parameters for these simulations are given in Table A1. Note that although the dense sand considered in these experiments contains some porosity at the reference density 0 , for simplicity it was assumed that the minimum distension min = 1.0 (i.e., s0 = 0 ) and that only pore space added during shear resulted in an increase in distension (i.e., = ∕ 0 ). In other words, permanent porosity change owing to compaction was not considered and, instead, a simple polynomial equation of state was used to relate pressure to density. Shown are stress ratio ( 1 ∕ 3 ) and volumetric strain versus axial strain for dense sand at three different confining pressures (0.05 MPa, 0.6 MPa, and 1 MPa). The experimental data exhibit three regimes with increasing strain: an elastic regime (< 1% strain) where stress rises linearly with strain and volume change is negative; a so-called hardening regime (between < 1% and ∼3% axial strain) where permanent failure begins and volume change transitions from negative to positive, resulting in a strength increase; and a plastic regime (>3% strain) where stress and positive volume change approach approximately constant values. The absence of a hardening model implies that the dilatancy model is unable to fit the data in the intermediate, hardening regime. Accounting for the dependence of porosity on the Poisson's ratio would improve the model in this regime. However, the initial elastic response and the large-deformation (plastic) response is adequately represented at the confining pressures considered.