Impact cratering is the dominant surface process on most solid bodies of the solar system. The process of hypervelocity impact causes irreversible changes to the nature and physical properties of rocks (Melosh, 1989). These changes produce anomalies in geophysical data that make craters conspicuous in comparison to the surrounding rocks. The identification of potential craters, on Earth and, increasingly, across the 10.1029/2019JE005929 Solar System (e.g., Zuber et al., 2013) is primarily achieved by remote geophysical surveying. The most commonly observed geophysical signature of an impact structure is a circular region of negative residual gravity (Pilkington & Grieve, 1992). Negative residual gravity anomalies at craters are caused by the replacement of high-density intact rocks in the preimpact target with fractured para-autochthonous rocks, porous impact breccias, and/or the relative loss of mass due to excavation and basin infill (whether it is filled by atmosphere, water, post-impact sedimentary rocks, or vacuum). Within complex craters, which are shallow compared to their diameters, the main contributor to the gravity signature derives from fractured para-autochthonous target rocks beneath the floor of the crater (Pilkington & Grieve, 1992). Additionally, large complex craters (>40 km) may have relative gravity highs in their center due to the uplift of high-density material from depth during cratering. The presence of low-density, porous rocks in impact structures has important implications for the transport of fluids within craters: affecting the viability and longevity of hydrothermal activity (e.g., Abramov & Kring, 2007), the distribution of potential hydrocarbon reservoirs (Grieve, 2005), and the habitability of impact structures for microbial life (Cockell, 2006). Consequently, understanding how small-scale properties of impactites link to the large-scale geophysical characteristics of craters is of critical importance.
Impact simulations are now capable of modeling the tendency of rocks and granular materials to change volume during shear failure, that is, dilatancy. Results of impact simulations using this model are consistent with porosity measurements and the gravity anomalies of simple craters (Collins, 2014). However, limited porosity data have hindered comparisons with complex craters. Furthermore, while the importance of fracturing in reducing the density of para-autochthonous impactites has been widely acknowledged, detailed analysis of the microporosity and microfracturing of shocked para-autochthonous rocks (e.g., Winkler et al., 2018) is limited, particularly for naturally shocked rocks.
Here, the Chicxulub impact structure is used as a case study to understand the production of porosity during complex crater formation. We analyze the petrophysical and microstructural character of uplifted basement rocks in the Chicxulub peak ring and, through numerical impact simulations, link our results to the gravity signature across the crater. Our study provides insight into the formation of complex craters and, more specifically, the processes that produce the remarkable physical properties of impact-deformed rocks.
The Chicxulub structure, located in the Yucatán peninsula of Mexico, is a 190-to 210-km diameter impact structure (Figure S1b in the supporting information). Its identification was first made on the basis of large-scale geophysical signatures: a large, circular, negative Bouguer gravity anomaly and magnetic anomalies (Hildebrand et al., 1991;Penfield & Camargo, 1981). The only surface expression of the crater's presence is a semicircular ring of water-filled sinkholes, cenotes, in the Cenozoic limestones that overlie the onshore portion of the crater.
Following the discovery of the crater, several attempts were made to constrain the size and structure of the Chicxulub crater based on forward modeling of potential field data (Espindola et al., 1995;Hildebrand et al., 1991Hildebrand et al., , 1995;;Sharpton et al., 1993Sharpton et al., , 1996)). The variations between the models, in terms of structure and the density contrasts between lithologies, attest to the nonuniqueness of gravity modeling.
Seismic reflection and refraction surveys were carried out at Chicxulub to resolve the crater dimensions and subcrater structure (Gulick et al., 2008;Morgan et al., 1997). The results of those surveys indicated that the topographic (inner) rim of the crater is 140-170 km in diameter (beyond which are outer ring faults), with a floor approximately 1 km beneath the surface, and an 80-to 90-km diameter peak ring, which rises up to 400 m above the crater floor (Gulick et al., 2013). Beneath the crater center, the mantle has been uplifted by 1.5-2 km (Christeson et al., 2009). The basement, both outside and deep beneath the crater, possesses seismic velocities >5.8 km/ (Christeson et al., 2001), indicating that the preimpact crust, beneath the sedimentary cover, consisted of high-density, nonporous rocks. A ring-shaped Bouguer gravity minimum within the crater corresponds to a radial distance slightly interior to the peak ring. Additionally, beneath the peak ring a strong inward dipping reflector can be observed. Together, these observations indicate that the rocks of the peak ring, and its roots, which dip inwardly, are among the lowest density para-autochtonous rocks within the crater. A detailed review of the geophysical characteristics of the Chicxulub structure can be found in Gulick et al. (2013).
In 2016, a joint expedition of the International Ocean Discovery Program and International Continental scientific Drilling Program (IODP-ICDP Expedition 364) drilled into the peak ring of the Chicxulub struc-10.1029/2019JE005929 ture (Figure S1c). The expedition recovered core between 505.42 and 1334.69 mbsf (meters below sea floor) including, from 750.25-1334.69 mbsf, a sequence of uplifted granitic target rocks, with occasional sheet intrusions of preimpact igneous rocks and impact-related breccias and melt rocks. A schematic lithological column and location map is shown in Figures S1a and S1b, respectively. Physical property data from the expedition (Christeson et al., 2018;Morgan et al., 2017) have shown that the shocked target rocks within the Chicxulub peak ring have porosities of approximately 10%, which is a likely cause of the gravity minimum associated with the peak ring. However, the microstructural cause of this porosity is uncertain and it is not known how porosity is more widely distributed across the crater.
Here, we assess the cause of the unusually low densities and high porosities of the uplifted granitic rocks of the Chicxulub peak ring. These results help constrain cratering-related deformation mechanisms and the timing of porosity generation. Observational results are obtained by petrophysical and petrographic techniques, with additional data from remote geophysical surveys. The observational data are then compared to the results of numerical impact simulations, which can predict porosity distribution and microfracture orientations.
2.1. Geophysical and Petrophysical Data 2.1.1. Gravity Data Numerous gravity surveys have been acquired at the site of the Chicxulub structure. Data used here (Figure S1b) were provided by Alan Hildebrand and Mark Pilkington (see Pilkington et al., 1994), supplemented by data acquired during the 2005 seismic surveying of the structure (see Gulick et al., 2008). Comparison between the results of numerical impact simulations and gravity data requires the determination of an average gravity profile across the structure, without any regional gravity effects. Following Vermeesch et al. (2009), data were initially transformed to a cartesian coordinate system with an origin at the crater center (21.31
• N, 270.46 • E): X = (101.71 × Longitude) + 9, 107.56, (1) Y = (112.02 × Latitude) -2, 386.67.
(
The gravity data were then corrected to account for the regional variation of the Bouguer gravity field across the site. This correction was achieved by excluding all data within 100 km of the crater center, and fitting a best fit plane through the remaining regional data set. The best fit plane, representing the regional Bouguer gravity field, Z, takes the following form:
After correction of the Bouguer gravity data to account for the regional field, radial profiles were acquired in 1 • intervals. These radial profiles avoided anomalous regions of the gravity field to the south and northwest of the crater center (Figure 1a). From the 232 remaining radial profiles, the median and quartile values of the radial residual gravity field were determined with increasing radial distance (Figure 1b). The effect of removing profiles from the anomalous regions can be seen in the supporting information (Figure S2).
IODP-ICDP Expedition 364 collected density data from cores using two different methods. First, cores were analyzed with a Geotek Multi-Sensor Core Logger (MSCL), which, alongside other physical property measurements, measured bulk density by gamma ray attenuation. These measurements were acquired at 2-cm intervals where, due to the machine setup, the individual resolution was 0.5 cm (see Morgan et al., 2017, for details). The second set of density data derives from He-pycnometry measurements (Quantachrome pentapycnometer) of discrete samples taken from the core. Samples had a diameter of 2 cm and an approximate volume of 6 cm 3 . Samples were taken approximately once per core section (approximately every 1.2 m), and a total of 719 samples were measured. He-pycnometry, in addition to allowing determination of bulk density, also allows the determination of grain density and porosity (see Morgan et al., 2017, for details).
The raw expedition data are shown within the supporting information (Figure S3). Here, we have processed the data acquired during Expedition 364 in order to facilitate direct comparison between the results 10.1029/2019JE005929 of numerical simulations and the observations, and between the two observational methods. First, all data from lithologies that were not granitic target rocks were excluded; numerical simulations are not capable of modeling the injection of allochthonous impactites into the target rocks nor can they model fine-scale lithological variations between subvolcanic intrusive rocks and basement rocks. Subsequent data processing involved the calculation of a running average and standard deviation within ±6 m of each He-pycnometry sample. Averages and standard deviations were only calculated where three or more measurements were acquired within the 12-m-wide interval. Comparable measurements from the MSCL data were acquired by taking an average and standard deviation corresponding to the same depth interval for each He-pycnometry sample.
To understand the microstructural cause of the porosity within the granitic basement rocks, 10 thin sections were prepared from samples of the recovered granitic target rocks at different depths within Hole M0077A. Backscattered electron images of the thin sections were obtained using a Zeiss Ultra Plus field emission scanning electron microscope (SEM) with the large-area imaging software and hardware package, ZEISS Atlas 5, at the Natural History Museum, London. Individual images were obtained with an accelerating voltage of 20 kV using the backscattered electron detector, a working distance of 8.0 mm, and a resolution of 748.9 nm per pixel. The Atlas 5 software was then used to stitch the individual images from a section to produce a single composite image. The total areas of each stitched image are shown in Table S1. As a consequence of the pixel size, and due to Nyquist's Theorem, the minimum width of a pore that can be accurately resolved in these images is ∼1.5 μm.
Image processing was carried out in order to analyze and quantify porosity within the thin sections. Processing was achieved using the ImageJ software package (Schindelin et al., 2012). The first stage required to analyze porosity is image thresholding. The images were deliberately acquired such that there was a significant contrast of grayscale values between pores (black), and mineral phases (gray to white). ImageJ provides several algorithms that automatically thresholds images into objects (pores) and background (solid material), resulting in an image consisting of only binary values, black for pores and white for background. Here, the default ImageJ thresholding algorithm was used, and the result then manually checked and, if required, adjusted.
Once pores were thresholded into a binary image file, porosity estimates could be made by assuming the application of the Delesse principle that the volume fraction of pores can be estimated by calculating the area fraction of pores in a 2-D image. Additionally, the orientation and aspect ratio of each pore were determined by fitting ellipses to each pore, an available function in ImageJ.
We simulated the formation of the Chicxulub impact structure using the iSALE shock physics code. iSALE is a multirheology, multimaterial code based on the SALE hydrocode (Amsden et al., 1980). The original code has been modified to include an elastoplastic constitutive model, fragmentation models (Melosh et al., 1992), various equations of state (Ivanov et al., 1997), a porous compaction model (Wünnemann et al., 2006), and a dilatancy model (Collins, 2014). iSALE and its precursor codes, SALES-2 (Collins et al., 2002) and SALE-B (Ivanov, 2005), have been used to simulate the formation of the Chicxulub impact structure in previous studies (Collins et al., 2002(Collins et al., , 2008;;Ivanov, 2005;Morgan et al., 2016;Rae et al., 2019). Here, we use the model parameters presented by Rae et al. (2019), where a 12-km diameter spherical impactor, using an equation of state for granite with a density of 2,630 kg/m 3 strikes a three-layer target at 15 km/s. The three-layer target is composed of a 3-km layer of sedimentary rocks, represented with an equation of state for calcite, overlying a 30-km layer of crystaline basement rocks, represented with an equation of state for granite, overlying mantle, represented with an equation of state for dunite. This simulation was run at a resolution of 60 cells per projectile radius, that is, cell widths of 100 m. A complete list of the parameters used in our simulation is shown in the supporting information (Tables S2 and S3).
Predicting the gravity signal of an impact structure in a numerical model requires a dilatancy model. Dilatancy is the tendency of granular and rocky materials to undergo volume change when subjected to shear deformation. Collins (2014) describes the dilatancy model used in iSALE. The equation used to update the distension in iSALE computational cells undergoing shear deformation is
where the first term, 𝛼, is the distension, which can be defined as a function of porosity, 1 1-𝜙 . The final term of equation ( 4) is the plastic shear strain rate, which is determined based on velocity gradients in the cell. The second, and most important term in equation ( 4), is known as the dilatancy coefficient, 𝛽, which describes the material's tendency to gain volume upon plastic shear strain, and is a function of pressure (p), temperature (T), and the preexisting distension (𝛼) of the material (see Table S3):
𝛽 max is a material-specific value, where large values indicate a large tendency to dilate upon shear deformation. For impact cratering simulations, and large-scale numerical simulations of rock masses, Collins (2014) suggests that 𝛽 max values range between 0.045 and 0.180, where 𝛽 max = 0.045, corresponding to a Low Geological Strength Index (GSI) material, provides the best fit to the gravity and porosity data of simple craters and small complex craters. Here, that value is used to test its applicability to larger complex craters.
From a simulation with this dilatancy model, it is possible to calculate the resultant gravity anomaly associated with the simulated crater including the effects of density reduction due to shear-induced dilatancy and the filling of the impact crater with postimpact sedimentary rocks. Gravity profiles are calculated by using the modeled distension field generated by iSALE; distension values are converted to a mass difference by assigning a reference density to each of the materials in the simulation, this assigned reference density can be varied independently from the reference density in the equation of state of the material. The additional effect of low-density crater fill can then be included by filling, to whatever depth required, the modeled crater with a material of an assigned density. Finally, the effect of erosion can be considered by removing any material above a specified depth. High-density impact melt rocks are accounted for in this calculation because materials approaching or above their melting temperature have suppressed dilatancy. However, due to the simplifications of the model, preimpact vertical variations in crustal composition, and therefore density, are not accounted for. Crustal density is generally expected to increase with depth; therefore, the gravity calculation is likely to underestimate the mass excess of the central uplift and overestimate the excess outside of the central uplift.
The processed residual 2-D gravity profile (Figure 1b) possesses a local gravity high in the crater center of -15 mGal. At 30-to 40-km radial distance, a gravity low of approximately -25 mGal can be observed, beyond which the gravity signal increases to match the regional gravity field. Beyond 80 km, the interquartile range of the radial gravity signal increases due to thickness variations of the Mesozoic and Cenozoic sedimentary rocks at various azimuths. The region to the northeast of the crater contains a thicker sequence of sedimentary rocks and therefore has a more negative residual gravity signal. The features of the average gravity profile are broadly consistent with the features seen in previously published individual profiles and average profiles (e.g., Sharpton et al., 1993). Additionally, the total mass deficit associated with the Chicxulub impact structure can be calculated from the gravity anomaly by applying Gauss's theorem (see Campos-Enriquez et al., 1998). The average gravity profile obtained here suggests a total mass deficit of 9.15 × 10 15 kg. Upper and lower bounds on the deficit can be made using the upper and lower quartile profiles, which suggest mass deficits between 4.90 × 10 15 and 1.25 × 10 16 kg, respectively. These values are consistent with previous estimates of the Chicxulub mass deficit (Campos-Enriquez et al., 1998).
The low seismic velocities of the Chicxulub peak ring, combined with the large gravity low slightly within the peak ring, had suggested, prior to drilling, that the Chicxulub peak ring is composed of unusually low density rocks (Morgan et al., 2000). Results from Expedition 364 have shown that the peak ring is composed of uplifted crystalline target rocks that possess unusually low densities and high porosities compared to typical crystalline basement rocks (Christeson et al., 2018). Processed expedition data are shown in Figure 2.
The target rocks possess consistent and typical grain densities for granite, averaging 2,628 ± 39 kg/m 3 (1𝜎), but have remarkably high porosities. The average porosity of the granitic target rocks is 11.5 ± 4.7%, resulting in average bulk densities of 2,444 ± 75 kg/m 3 ; additionally, the raw porosity measurements show that the granitic target rocks have a pervasive baseline porosity of approximately 8%, where few samples possess less porosity, but where local excursions to higher porosities occur. There is some heterogeneity of porosity and bulk density with depth in the recovered core: in the upper 100 m, average porosities increase beyond 20% and the variability increases up to 𝜎 = 8.1%. Additionally, porosities within granitic rocks beneath 1,225-mbsf increase from <10% to >10%. Despite some discrepancy between bulk density measurements by He-pycnometry and the independently measured bulk density by gamma ray attenuation using the MSCL, particularly from depths below 1,075 mbsf, the 1𝜎 envelopes of both measurements always overlap, indicating full consistency between the data sets.
Inspection of the SEM images shows that the largest contributor to the total porosity of the granitic target rocks is by intragranular fracturing, with additional contribution from cataclasites and large intergranular pores (Figure 3). Cataclasites are easily distinguished by their finer grain sizes, increased grain angularity, and increased porosity relative to their host material, while large intergranular pores are defined by their anomalous sizes compared to the other sources of porosity (see supporting information Figures S4-S13). Across all the thin sections, the average fraction of the porosity contributed by intragranular fractures is 0.55, while the average fraction of porosity contributed by intergranular pores and cataclasites are 0.24 and 0.21, respectively.
Excluding the large intergranular pores and cataclasites, the average porosity of the intragranularly fractured granites is 6.4%. This intragranular fracturing is distributed heterogeneously throughout the thin sections and is primarily dependent upon the mineral that is fractured (Figure 4). The three major rock-forming minerals within the granites are alkali feldspar, quartz, and plagioclase feldspar. Quartz grains are pervaded by planar fracture (PF) sets, which rarely span the entire grain and which occasionally have feather feature lamellae (Figures 4a and 4e; Poelchau & Kenkmann, 2011). Additionally, the quartz grains possess occasional individual shear fractures on planar surfaces with observable submillimeter displacements, and rare nonplanar fractures, which typically contain a gouge of angular to subangular quartz fragments. Compared to the feldspars, the quartz grains are more porous. Plagioclase feldspar, while also pervasively microfractured, possesses a distinctly different fracture pattern to the quartz grains. The microfractures are typically closer spaced than in quartz from the same sample, are rarely planar, and typically branch and converge into other fractures (Figures 4b and 4d). While the plagioclase crystals often have higher fracture densities than the quartz grains, the total dilation on each intragranular plagioclase fracture is significantly smaller than the intragranular fractures in quartz. Compared to the quartz and plagioclase feldspar within the granitic rocks from Hole M0077A, alkali feldspar is remarkably unfractured (Figure 4c), possessing only widely spaced nonplanar fractures. Additional porosity in the granitic rocks is produced by cataclasites. These regions of the thin sections have extremely high porosities. Across the four thin sections that contain cataclasites, the average porosity in the cataclasite regions is 15.8%. Nevertheless, the cataclasites are only a small area of the total rock mass (20.0% of the total imaged area) and thus contribute only a small fraction of the total porosity of the thin sections.
The SEM images show examples of cross-cutting relationships between the cataclasites and grains with intragranular fractures (Figure 4f) that show the cataclasites formed after pervasive fracturing of the target rocks, also described by Riller et al. (2018) . Additionally, Riller et al. (2018) have shown that the cataclasites are truncated by later deformation events associated with crater formation.
The large intergranular pores observed in the samples are almost certainly not representative of the in situ porosity of the granitic target rocks. Their extreme sizes within the pore size distributions (see supporting information Figures S4-S13) of each sample indicate that they are mostly nonnatural pores, that is, produced during sample preparation. Nevertheless, the possibility remains that thin sections are too small to encompass the full pore size distribution of these rocks. Regardless of their origin, the total contribution of these large intergranular pores to the total porosity is equal to the average contribution from cataclasites and significantly less than the contribution from intragranular fractures. These pores can only be found in 4 of the 10 thin sections, all deriving from near to the top of the hole. Total porosity measurements obtained from SEM image analysis, subdivided by the source of porosity, are shown in Figure 5. In general, the estimations from SEM image analysis are consistent with the measurements made by He-pycnometry. All results from SEM analysis, excluding the uppermost thin section, fall within the 1𝜎 envelope of the He-pycnometry porosity measurements. However, most of the SEM-derived porosities are offset to lower values compared to the mean measurement from He-pycnometry, and this would be exaggerated by excluding the large intergranular pores, many of which are due to sample preparation. Potential errors associated with calculating porosity from SEM images include: lack of resolution to quantify pores smaller than ∼1.5 μm, large pores distributed on a scale larger than a thin section (2-3 cm), porosity anisotropy, and sample preparation artifacts. Quartz develops a number of unique deformation features over a range of shock wave pressures (French & Koeberl, 2010). One of these features is PFs, which can occur in one or multiple orientations. The large abundance of PF sets in quartz from Hole M0077A, together with the interconnectedness of the intragranular fractures, suggests that most of the intragranular fracturing is related to deformation during shock (Figures 4a, 4c, and 4e). The partial disaggregation of some grains of quartz and plagioclase feldspar along intragranular fractures into cataclasites and the matrix of monomict granitic breccias indicates that the formation of cataclasites and brecciation occurred after the formation of intragranular fractures (Figure 4f).
The orientation of pores within the thin sections is generally anisotropic (Figure 6). Six thin sections have unimodal pore orientation-frequency distributions, while two more thin sections have a bimodal distribution (Table 1). Comparing these distributions to the images they derive from (Figure S14) indicates that the cause of the high frequency orientations is the orientation of PFs in quartz. 10.1029/2019JE005929 a The thickness of postimpact sedimentary rocks above the peak ring was maintained at 650 m.
Our numerical simulations of the Chicxulub impact event follow a model of the formation of peak rings, sometimes called the Dynamic Collapse Model, that has been developed over a number of decades (Collins et al., 2002;Grieve et al., 1981;Kring et al., 2016;Morgan et al., 2016). Within this model, and our simulation, the formation of a peak ring results from the overheightening of a central uplift, and its subsequent lateral emplacement over the collapsed transient cavity rim (Figure S15). The simulation presented here, from Rae et al. (2019), allows predictions of the distribution and timing of porosity to be made. Porosity generation in rocks during complex crater collapse would be expected throughout the crater due to shear failure. The final distribution of modeled porosity caused by this dilatancy is shown in Figure 7a.
From the final simulated distribution of porosity in the structure, it is possible to produce a vertical profile of porosity in the peak ring that can be compared to observational results (Figure 7c). Using 𝛽 max = 0.045, these results predict that the porosity of the peak ring material should be approximately 7 ± 2%. Uncertainty in the exact comparative location between the model and the drilled core results in the ±2% error. Overall, the modeled porosity closely matches the baseline values of porosity in the recovered core, with notable excursions within the top 200 m, from 300-350 m, and at 700 m within the peak ring. These excursions from the model may be caused by postimpact hydrothermal dissolution, and/or localization of strain that is not captured by the model.
In addition to the simulated porosity in the peak ring, a set of simulated Bouguer gravity profiles of the crater were calculated from the iSALE simulation and are compared to the observed local Bouguer gravity anomaly (Figure 7b). The range in the displayed profiles results from propagating the uncertainties in the assigned reference densities for the materials in the model (Table 2).
One important problem with the calculation of the overall gravity profile comes from the filling of the basin with postimpact sedimentary rocks. In the model, the top of the peak ring is at 1,050-m depth, while in the real crater, the top of the peak ring is at 618 mbsf. The cause of this discrepancy could be due to postimpact relaxation, or inaccuracies of the mechanism of transient weakening. To ensure that the mass deficit from the postimpact sedimentary rocks is consistent with observations, one of two changes could be made: (1)
The thickness of the sediments in the model could be maintained, but the reference density of the sedimentary infill increased such that the mass deficit remains constant, or (2) the reference density of postimpact sediments could be maintained at reasonable values but the thickness of the sedimentary infill reduced such that only 600 m of sedimentary rock overlie the peak ring. Here, we adopt the latter approach because it is more useful for predicting the gravity structure around the crater center. This approach will result in the model predicting a smaller magnitude anomaly in the terrace zone compared to the observational data due to the lack of low-density sedimentary rock at those radial distances.
The modeled gravity signal possesses a mass deficit in the range of 5.81-9.25 × 10 15 kg, considering the uncertainties presented in Table 2. This range includes the value of the mass deficit associated with the median gravity signal, 9.15 × 10 15 kg, and is fully within the upper and lower quartile gravity signals. The maximum modeled anomaly at the crater center is -9 to -28 mGal, consistent with the observed value of ∼ -15 mGal. Additionally, the maximum anomaly in the modeled profile occurs between 28-and 40-km radial distance and has a value of -22 to -34 mGal, consistent with the observed -26 mGal minimum at 35-km radial distance. Beyond the peak ring, within the annular trough and the terrace zone of the crater, the fit between the modeled and observed gravity profiles is poorer. This misfit can be explained primarily by the oversimplifications associated with modeling the postimpact sedimentary fill as a constant density material. The gravity anomaly is likely to be underestimated between 45-and 55-km radial distance, since deep sedimentary rocks in the annular trough adjacent to the peak ring are likely to have higher densities, as suggested by their increase in velocity with depth (Morgan et al., 2011). The overestimation of the observed profile from 55-to 75-km radial distance occurs, as previously stated, due to the lack of sedimentary infill at these radial distances. Additional causes of discrepancy between the model and observation may include: postimpact compaction, hydrothermal precipitation and/or dissolution, or a lack of initial porosity in the modeled preimpact sedimentary rocks.
The rocks of the Chicxulub peak ring, with porosities of 11.5 ± 4.7%, have low densities compared to typical crystalline basement rocks, which usually possess negligible porosity. Seismic velocities in the basement rocks of the region (Christeson et al., 2001) indicate that the preimpact porosity of the crystalline target rocks was small, ≪1%. Thus, compared to shocked target rocks in other craters (Pilkington & Grieve, 1992), the target rocks of the Chicxulub peak ring have among the highest impact-generated porosities of any impact structure (Table S4), that is, almost all of the measured porosity in the crystalline rocks of the Chicxulub peak ring was produced by impact processes or postimpact alteration. We attribute the remarkably high porosity to the relative size and age of the Chicxulub impact, sampling depth within the structure, and the complex deformation path that peak ring materials follow during dynamic collapse (Rae et al., 2019).
The main contributor to the porosity of target rocks in the Chicxulub impact structure is shock-induced microfracturing. Additional porosity is contributed from cataclasites, which are extremely porous but volumetrically small, and, possibly, large intergranular pores. The origin of large intergranular pores is almost certainly due to preparation-induced damage. Nevertheless, large intergranular pores are systematically prevalent near the top of the granitic rock section, coinciding with the largest deviation between the modeled porosity and observed porosity.
Most peak ring material never experiences tensile stress regimes during cratering (Rae et al., 2019). Consequently, the shear-induced dilatancy model used here should simulate the production of all styles of primary microporosity observed within the peak ring rocks. Secondary porosity may be generated by processes such as hydrothermal dissolution. The observational evidence shows that the porosity of the peak ring rocks is primarily hosted in intragranular fractures formed by shock metamorphism. In contrast, the dilatancy model predicts that very minimal amounts of porosity are produced during the passage of the shock wave and rarefaction. Instead, it predicts that, for the shallowest rocks in the peak ring, porosity is generated in two distinct phases of shear deformation, the first during transient cavity formation and rim collapse (0-100 s) and the second during central peak formation and collapse (175-350 s; Figure 8b). Deeper material within the peak ring structure experiences only one distinct phase of shear-induced dilatancy, which occurs more gradually through transient cavity formation and collapse. One explanation for the discrepancy between the 10.1029/2019JE005929 observation of shock-induced microfractures dominating porosity and the model prediction of later porosity generation by shear is that fractures form early, during high-pressure shock, without any porosity. Those fractures must then open during later shear-dominated phases of crater deformation. However, it is also possible that the shear-induced dilatancy model used here is an oversimplification of the process that generates porosity in para-autochthonous impactites. A process that may involve two distinct stages of porosity generation; first, a dynamic effect of fracturing during shock and the second due to true shear-induced dilatancy of a predamaged material. If this is the case, the GSI parameters used here for shear-induced dilatancy are too large.
The ability to accurately predict the porosity distribution of complex craters is important for understanding the gravity signal of large impact structures. Here, we have obtained a reasonable fit to the observed gravity signature of the Chicxulub structure, a large complex impact crater, using a numerical model that accounts for porosity generation/density reduction by dilatancy alone. Discrepancies between the modeled and observed gravity signals may be due to pore space compaction or displacement of preexisting crustal density variations not captured by the model. Additionally, changes to acoustic fluidization parameters may change the location of down-faulted sedimentary rocks in the collapsed transient cavity rim that underlies the peak ring and therefore affect the modeled gravity signature. Despite this, the overall consistency 10.1029/2019JE005929
between the gravity profiles and porosities between observation and model reinforces the significance of dilatancy on subcrater structure and deformation. Additionally, our results support the recommendation by Collins (2014) that low GSI dilatancy parameters should be used for impact cratering simulations.
In addition to the consistency between porosity and gravity measurements with the results of numerical impact simulations, it is also possible to predict the orientation of shock-induced microfractures in numerical impact simulations (Rae et al., 2019). Predicting the orientation of microfractures during cratering requires the assumption of two well-accepted concepts of shock physics and rock mechanics. First, the orientation of fractures during rock failure is predicted by the canonical model of Mohr-Coulomb rock failure (Anderson, 1905). In low-pressure regimes, rock failure is typically expected at 30 • to the orientation of 𝜎 1 , depending on the coefficient of internal friction. However, in a high-pressure regime, greater than the Hugoniot Elastic Limit, the coefficient of internal friction approaches 0, and consequently, the angle of failure approaches 45 • to 𝜎 1 . Second, the formation of PFs during shock occurs during shock loading, that is, once the Hugoniot Elastic Limit is exceeded but before peak pressure conditions are reached (French & Koeberl, 2010;Poelchau & Kenkmann, 2011). Support for this assertion derives from the observation that PFs commonly separate domains within grains that Planar Deformation Features, which are sensitive to the peak pressure of the shock wave, do not cross. Using these concepts, and by obtaining the orientations and magnitudes of deviatoric principal stresses, following the methodology described by Rae et al. (2019), it is possible to predict shock fracture orientations within target material during shock (Figure 9).
During crater excavation and modification, and to a much lesser extent, during shock, material is expected to rotate from its initial orientation. The orientation of material in numerical simulations can be tracked by calculation of infinitesimal strain tensors (see Rae et al., 2019). Here, the Dynamic Collapse Model has been used, and a single tracer selected to represent the recovered core material. This particle rotates by approximately 60 • clockwise during cratering (Figure 10c). Selecting tracers from deeper or at increasing radial distances within the peak ring material causes a systematic increase in the total clockwise rotation during cratering by up to an additional 30 • . Consequently, the final orientation of shock-related fractures, as predicted by the model, is at concentric strikes, dipping at 88 • inward-78 • outward, and the conjugate set at 12 • inward-2 • outward (Figure 10b). The observed orientations of microfractures within the shocked granitic rock samples can be subdivided into one set of steep apparent dips (69.4-89.9 • ) and a second set with shallow apparent dips (7.7-13.3 • ). Despite the observations only being capable of determining apparent dip, without any strike information, the modeled orientations are remarkably consistent with the observed micro-fracture orientations.
This result indicates several important points. First, that the orientation at which shock microfractures, that is, PFs, forms may be primarily controlled by the orientation of principal stresses during shock. The granitic rocks from Hole M0077A have minimal to no grain preferred orientation (Morgan et al., 2017;Riller et al., 2018); consequently, the preferred orientation of shock microfractures across a whole thin section of randomly oriented grains indicates that crystallographic orientation may not have primary control of the orientation of fractures. PFs are known to form on rational crystallographic planes (French & Koeberl, 2010); however, it is likely that the specific orientation of the plane that is activated during shock may be controlled by stress conditions.
Second, it is clear that a large rotation of the peak ring material is required between shock and the final emplacement of the peak ring. Regardless of the physical mechanism by which crater collapse occurs, the orientation of fractures due to shock will always be at moderate dip angles (approximately 45 • ) immediately after shock. Achieving conjugate sets of planar fractures in subvertical and subhorizontal respective orientations requires either a clockwise rotation of material by approximately 60-70 • or an anticlockwise rotation of 20-30 • . The Dynamic Collapse Model of peak ring formation precisely predicts the large clockwise rotation required.
Finally, stress-induced fractures significantly affect the permeability anisotropy of rocks. Canonically, the maximum permeability direction of a fault rock is along the intersection of minor faults, that is, parallel to the orientation of the 𝜎 2 direction that produced the fractures (Sibson, 2000). Combining this with the fact that shock microfracturing produces the dominant pore fabric in these peak ring rocks leads to possibility that numerical impact simulations may be able to predict the permeability distribution and anisotropy within impact craters, a critical parameter of hydrothermal models of impact structures.