A 1 m² area within our larger, long-term excavation project (totaling over 200 m² to date) was selected for detailed examination in this study (Figure 1A). The size and map-view location of each fossil in this square meter area were recorded by hand on a 1:1 scale paper map at the field site. After removing approximately 1 m of overburden, every fossil unearthed from 55 cm above and 5 cm below the Navesink-Hornerstown contact (NHC), regardless of size or completeness, was mapped (with a specimen’s stratigraphic position recorded based on its highest point for feasibility; see Supplementary Material for more details), collected, prepared, and identified to the least inclusive taxonomic unit possible. All fossils are reposited with unique identifiers in the Edelman Fossil Park Collection at Rowan University.

Plumb bob measurements (Figure 1C) were converted to positions above or below the NHC and graphed against specimen abundance using 1 cm intervals to track stratigraphic position. To examine fossil distribution patterns at several spatial scales, we explored trends in the data across the whole square meter as well as in spatial subsets of halves and quadrants (as portrayed in Figure 1B). Peaks in fossil abundance data for the entire square meter were delineated by counts >40; therefore, > 20 was used as a cutoff point when examining the data in spatial halves and ≥10 was used when examining the data in quadrants. Simpson Diversity Index (SDI) and taxonomic richness values were calculated using all fossil fragments because assignment of isolated and disarticulated fossils to the same individual is rarely possible at this locality (Boles, 2016) and many fossils cannot be identified to the species level due to incompleteness. For more information on how these metrics were calculated, see the Supplementary Material. For classification purposes, Xylophagella specimens were counted within the “plant remains” category (rather than the “bivalve” category) because they consist of trace fossils (boring casts) formed within driftwood without any portion of the pholadid bivalve remaining (Gabb, 1860).

The 3D position of each fossil was further analyzed in ArcGIS Pro in order to capitalize on the ability of GIS software to quantitatively evaluate and visualize three-dimensional (3D) spatial relationships. To enable these analyses, we first converted our paper field map to a digital data record through the process of georectification within a Geographic Information System. High-resolution digital photographs were rectified to a 1:1 scale cartesian grid using ArcMap 10.6. The visual centroid of each fossil was recorded and labeled with each specimen’s unique identifier, which were then joined by attribute to the data recorded in the field and laboratory. For all ensuing spatial analyses, taxonomic units were combined into broader taxonomic categories (e.g., Gastropoda, Chondrichthyes) to permit visualization across the volume. 2D visualizations of the spatial data were created to view all fossils, each taxonomic category separately, as well as by 1 cm stratigraphic bins.

Three-dimensional visualization was conducted in ArcGIS Pro 10.8 by creating a local scene for the same spatial reference grid and using the elevations of specimens in relation to the NHC. Concave hull polygons were created for the spatial data of each broad taxonomic category using the Minimum Bounding Volume tool, enabling visual comparison of the geometric trends of each group by transparently overlaying their enclosing geometries to examine internal vertices and those externally concealed by opaque overlay. A low stratigraphic peak in abundance (later termed the “oyster layer”) was visualized as a single concave hull of all oysters in the 2D-analysis defined range of 7–19 cm above the NHC. Additionally, we used the Near 3D Tool of ArcGIS Pro to calculate the 3D Euclidean distance from each Pycnodonte and identified those within a 1 cm search radius in any direction (see Supplementary Material for more details).

Since a stratigraphically-higher peak in abundance (determined to correspond to the MFL below) contains the highest concentrations of Cucullaea vulgaris clams, plant remains (primarily Xylophagella wood borings), and tetrapod bones (see below), these taxa were used to approximate the shape of this layer in 3D. To verify that these fossils comprise a single fossil horizon, as would be expected based on prior studies (e.g., Gallagher, 1993; Gallagher, 2003; Boles, 2016), we looked at the intersection of their individual concave hulls using the Intersect 3D Tool, exploring the apparent visual overlap of their spatial extents, and potential correlation of their relative slopes. Prior to this analysis, the highest Cucullaea fossil was identified as an outlier, and it was thus not considered. We intersected the tetrapod bones and plant hulls to create a new hull consisting only of the overlapping regions. By visualizing this new hull and the Cucullaea hull simultaneously, we were then able to compare their spatial extents with the 2D spatial data, individual concave hulls, taxonomic composition metrics, and specimen abundance data to delineate the number of distinguishable fossil layers and their stratigraphic bounds.

Minimum and maximum inter-quadrant slope models were calculated for each specimen abundance peak to describe their 3D shape. Slopes were calculated between the vertical endmember specimens of each quadrant (the topmost and bottommost points), with the minimum and maximum values in each quadrant thus constraining the upper and lower bounds of each fossil layer. For the stratigraphically-lower peak in fossil abundance, these slopes were calculated over several potential intervals suggested by the 2D statistical results in order to test the sensitivity of the stratigraphic cutoffs. These stratigraphic intervals included: 8–17 cm above the NHC, inclusive (including fossils found at either 8.0 or 17.0 cm); 8–17 cm above the NHC, non-inclusive (not including fossils found at either 8.0 or 17.0 cm, but all those within that range); 7–19 cm inclusive; and 7–19 cm non-inclusive. Slopes were also calculated for the interval suggested by the nearest neighbor calculations. For the stratigraphically-higher peak in fossil abundance, the slopes bounding the Cucullaea concave hull were compared with those of the tetrapod bones + plant intersection concave hull to determine whether these areas exhibit comparable slopes and thus could be considered to belong to the same fossil layer.

Finally, structural dip angles were calculated for upper and lower bounds based on the highest or lowest endmember points, respectively, of each stratigraphic peak in fossil abundance. Two methods were used to assign the oyster layer boundaries (fossil abundance and Nearest neighbor) and each of these intervals were used to calculate separate slope and dip values because the Nearest neighbor analysis provides an interval independent of any user biases for the region most densely filled with oysters.

Our excavation of a 1 m² by 60 cm volume yielded 1,306 specimens; 1,278 of these specimens are macrofossils, of which 15% are vertebrate remains. Fossil specimens range in size from ∼1 cm fish bones to a ∼11 cm Euclastes turtle costal, with most specimens ranging in thickness from 3 to 30 mm. These fossils were identified into 17 genus or species level taxonomic units (including four ichnospecies, with the oyster Pycnodonte dissimilaris being the most common), an additional seven broader taxonomic categories, as well as coprolites and unidentified burrow casts (Supplementary Data File S1). Among these fossils are many taxa known to be common in the MFL at Edelman Fossil Park (e.g., Pycnodonte, Cucullaea vulgaris, Odontaspis cuspidata; Gallagher, 1993; Gallagher, 2003) as well as less common taxa (e.g., Sphenodiscus lobatus, associated Euclastes wielandi remains, a partial hexanchid shark tooth), indicating that the fauna we recovered is reasonably representative of the fauna of the MFL as a whole.

When the distribution of fossil specimens is vertically collapsed into a single 2D plane (i.e., map view; Figure 1D), they appear to form a virtually uniform distribution across the entire square meter area. Only one instance of definitive fossil clustering was identified: a portion of an associated sea turtle was excavated entirely from the southwest quadrant, adjacent to where additional remains of this turtle were recovered during our ongoing excavation in the months preceding this project. Associated and partially articulated skeletons are known to be fairly common in the MFL (Gallagher, 1993; Boles, 2016; Ullmann and Carr, 2021), so this result is not unexpected. Additionally, the only regions exhibiting reduced fossil density or areas void of fossils occur primarily around the edges of our study area, almost certainly due to minor, inevitable erosion of these exposed unconsolidated sediments that occurred over the length of our excavation. Similarly, the boundary between the northeast and southeast quadrants also exhibits this pattern, as this region was unfortunately exposed as a wall over a winter and some minor spalling of the exposed greensand occurred despite our covering the area.

Figure 2 displays the stratigraphic distribution of fossil abundance in 1 cm bins above and below the NHC, which clearly displays a bimodal distribution within the lower Hornerstown Formation. Two peaks are present, in which fossil counts per centimeter exceed 40 for essentially continuous stratigraphic intervals (i.e., no more than a 1 cm gap wherein the specimen count drops below 40), suggesting the presence of two separate fossil layers. The stratigraphically-lower peak in abundance occurs roughly from 12 to 19 cm above the NHC and the upper peak occurs approximately from 23 to 32 cm above the NHC. However, the diversity of fossils is not consistent throughout these stratigraphic intervals: each peak in fossil abundance is composed of a distinct distribution of taxa (Figures 2, 3; Table 2).

Taxonomic groups such as plants (primarily Xylophagella wood borings) and tetrapod remains were found primarily near the top of the stratigraphically-higher peak in fossil abundance. Other, more common fossils, such as Cucullaea, tended to be found more evenly throughout this upper layer. Osteichthyan and chondrichthyan remains were found more or less evenly distributed throughout all 60 cm of our excavation. Although the oyster Pycnodonte was found throughout the entire volume of sediment excavated, its dominance in the lower abundance peak is readily noticeable during excavation and analysis (Figures 2, 4). These trends are in general agreement with anecdotal evidence from our larger excavation of over 200 m² across the same stratigraphic interval at this site.

Taxonomic richness values range from 0 to 10 in our study area across each 1 cm stratigraphic bin (Figure 3). Consistently-higher taxonomic richness values are present in the upper portion of the examined stratigraphic interval compared to in its lower portion. When stringent criteria are used to define peak intervals, then two primary richness peaks are recovered that correspond to the stratigraphic intervals of the lower and upper peaks in fossil abundance discussed above (see Supplementary Section S1.3 for further discussion of peak identification criteria). A subtle third peak in richness occurs near the top of the stratigraphically-higher richness peak, rebounding from a brief decline, wherein each 1 cm stratigraphic bin exhibits a richness ≥6. A three point running average of richness smooths some of the stochastic effects of varied fossil counts in each bin (Figure 3) and visually clarifies the extent of this subtle third peak in richness within the much more obvious, broad second peak.

Though the SDI becomes highly biased at low sample sizes, such biases are only seen in the stratigraphically highest and lowest ∼10–15 cm, outside of the intervals exhibiting peaks in fossil abundance (Figure 3). Of note, the stratigraphic interval exhibiting the lowest diversity values (<0.3) occurs from 12 to 17 cm above the NHC. This interval wholly overlaps with the stratigraphically-lower peak in fossil abundance found above. Conversely, when stratigraphic bins biased by “low” sample size (i.e., <15) are excluded, the interval exhibiting the highest diversity values (>0.5) occurs from 26 to 33 cm above the NHC. This interval wholly overlaps with the stratigraphically-higher peak in fossil abundance found above. Although these stratigraphic correlations are intriguingly consistent, it must be noted that stratigraphic bins exhibiting high fossil counts and high richness do not always also exhibit high diversity values (e.g., the bin from 18 to 19 cm above the NHC).

Percent composition of Pycnodonte, which was calculated as an additional means of tracking faunal diversity by 1 cm stratigraphic bins (see Methods), gradually increases down section (Supplementary Data File S1). Within the stratigraphically-higher peak in fossil abundance found above, the percentage of Pycnodonte varies from 28 to 62% of the fossils recovered. The stratigraphic interval exhibiting the consistently-highest percentage of Pycnodonte (>80%) occurs from 12 to 17 cm above the NHC; this interval overlaps with the stratigraphically-lower peak in fossil abundance found above.

As expected, the exact stratigraphic ranges of peaks in fossil abundance vary slightly among the examined spatial subsets compared to the whole square meter (Supplementary Data File S1). Two peaks were recovered in the south half, and they coincide in position with the two peaks identified within the full dataset (Figure 2). Though two peaks were also recovered in the north half, they occur slightly stratigraphically higher, and data from the east half exhibit three peaks. Data from the west half form a trend suggestive of the presence of two peaks, but the apparent lull between them does not satisfy the criterion to split the trend into two separate peaks (see Supplementary Section S1.1 for further discussion of peak stratigraphic intervals). To more finely characterize these spatial variations in the stratigraphic bounds of the two apparent fossil layers, we further parsed the specimen occurrence data into quadrants (Supplementary Data File S1). The stratigraphic intervals of peak fossil abundance recovered in each quadrant are detailed in section 1.1 of the Supplementary Material. As when the abundance dataset was partitioned into halves, quadrants also show variation in the number and stratigraphic position of peaks, with one to four peaks being identified among quadrants. Cumulatively, these results from spatial subsets reveal that the northern portion of our excavation exhibits both broader stratigraphic dispersion of fossils and slightly stratigraphically-higher peaks in fossil abundance than found across the entire square meter as a whole.

Here, we first consider concave hulls created from the positional data of specimens belonging to our broad taxonomic categories. Each concave hull calculated for a taxonomic category exhibits a unique 3D geometry (Supplementary Figure S2). Pycnodonte specimens form a relatively amorphous shape with two higher peaks in the northwest and southeast quadrants (Supplementary Figure S2A). Compared to Pycnodonte, chondrichthyan specimens form a visually “opposite” amorphous shape with lower protrusions in the northwest and southeast quadrants (Supplementary Figure S2B). The geometry of the tetrapod bone concave hull, in contrast, is relatively flat at the top of its stratigraphic bounds, with the highest specimens in the northwest and lowest specimens in the southwest quadrant (Supplementary Figure S2C). Similarly, plant fossils form a planar hull that is tilted stratigraphically upward to the northwest (Supplementary Figure S2D). The concave hull calculated for Cucullaea, both with and without the single high outlier removed, are also similar to those of tetrapod bones and plant fossil hulls in that they are overall flat on the bottom and curve along their tops with a northwest-trending positive slope (Supplementary Figure S2E). Burrows (Supplementary Figure S2F) and osteichthyan bones (Supplementary Figure S2G) form concave hulls that are irregularly dome-shaped, each with a gently-sloping base. Additionally, the osteichthyan bones hull exhibits a corrugated shape that gradually slopes upward toward the north. The gastropod (Supplementary Figure S2H) and other fossil specimens (Supplementary Figure S2I) hulls exhibit similar characteristics in that both are generally convex polygons with the lowest specimens present in the northeast.

Across the entire dataset, the average number of oysters whose nearest-neighboring oysters is within 1 cm is two, with a standard deviation of 2.3. The most dense horizon of oysters in 3D space, as identified from the stratigraphic interval exhibiting the highest frequency of oyster nearest neighbors, occurs at 16.0 cm above the NHC, with 18 oysters being within 1 cm of the nearest neighboring oyster (Figure 5); in other words, this identified the greatest abundance of closely-packed oysters occurs at that stratigraphic height. The next-highest frequency (12 oysters as nearest neighbors) occurs immediately below this maximum. From 10.6 to 19.1 cm above the NHC, frequency values consistently exceed seven, greater than two standard deviations above the mean; this is the only examined stratigraphic interval to consistently exhibit such high frequencies of neighbors. Stratigraphically beneath this region, Pycnodonte neighbor counts almost never exceed the mean, and immediately above this region occurs a 4 cm thick interval of oysters whose number of neighbors rarely exceeds one standard deviation above the mean. Above this 4 cm interval, the frequency of nearest neighbors occasionally spikes just over two standard deviations above the mean but never remains at these levels for successive heights. Thus, the stratigraphic interval from 10.6 to 19.1 cm above the NHC is statistically significant regarding its oyster density, even among other intervals within the study area which include abundant Pycnodonte fossils.

The stratigraphically-lower peak in fossil abundance, as defined by our 2D analyses, displays a relatively flat geometry. Minimum and maximum inter-quadrant slope models of this layer (Table 3) show the most consistent geometry from the 2D model occurs from 7 to 19 cm above the NHC. Slope values for the top and base of this layer are all very close to zero (Table 3). Dip angles calculated based on the concave hull and the nearest-neighbor-defined lower fossil layer (see Methods) were also found to be zero or near zero (Table 4), in general agreement with these low slope values.

The upper fossil layer was found to be more difficult to define in 3D space due to its greater diversity; therefore, we calculated the slope of this layer using multiple definitions (Table 3). First, the stratigraphic interval encompassed by this layer was estimated by analyzing the distribution of all tetrapod bones, plant remains, and Cucullaea specimens in the entire square meter. The lower boundary of this interval was placed one standard deviation below the mean height of all specimens of these taxa (31.3 cm above the NHC). Its upper boundary was placed at two standard deviations above the mean rather than one to account for the apparent upward skew of fossil occurrences of taxa commonly or exclusively found in this higher peak in abundance. Using these bounds, our 3D analyses found the stratigraphically-higher abundance peak to occur from 23.9 to 41.4 cm above the NHC. As two alternative definitions, we also considered: 1) the geometry of the Cucullaea hull alone (Supplementary Figure S2E) as a general representation of this layer, and; 2) the sum of the volumes encompassed by the Cucullaea hull and the tetrapod bones + plant intersection hull. On visual inspection, all of these hulls exhibit similar slopes that extend downward from the northwest to southeast along the top with an overall flat bottom, as evidenced by the higher maximum dip angles on the top compared to its base (Table 4; see Supplementary Section S1.4 for more details).

Together, these findings suggest that the shape of the upper fossil layer (assigned to be the MFL, see Discussion) is not a purely tabular layer, but rather it exhibits subtle variations in “topography” which include a slightly upward slope toward the northwest within our study area. In contrast, the lower fossil layer lacks such slope variations within our study area.

Data from our microstratigraphic analyses strongly suggest the presence of two distinct fossil assemblages exhibiting stratigraphic and taxonomic heterogeneity and, presumably, differing formational histories within sediments comprising the basal Hornerstown Formation. Unlike prior studies, our explicit intent was to characterize the centimeter-scale 3D distribution of fossils within this stratigraphic interval, and this unique focus allowed us to distinguish these two distinct fossil assemblages. However, we are not the first to propose differentiation of fossils within the basal Hornerstown Formation. Gallagher (1993) proposed that the MFL is composed of several layers including a basal oyster bed with overlying vertebrate remains. Our findings quantitatively corroborate that qualitative assessment as we have identified additional taxonomic and spatial distinctions between what have been occasionally regarded as “lower” and “upper” portions of the MFL.

Evidence from our 2D analyses for the presence of two fossil horizons in the basal Hornerstown stems from all three aspects we considered: fossil abundances, taxonomic richness, and taxonomic diversity. The strong bimodal distribution of fossil abundances across the entire square meter, with a lower peak from 12 to 19 cm and an upper peak from 23 to 32 cm above the NHC, is indicative of the presence of two fossil layers (Figures 2, 3). Bimodal distributions were also found in the north half, south half, and southwest quadrant, and the general character of fossil abundances in many other data partitions (e.g., the east and west halves) is also bimodal (Supplementary Data File S1). Although the west half and northwest quadrant each exhibit just a single peak in abundance, those broad peaks appear to meld the upper and lower peaks (i.e., the assemblages) observed in other data partitions.Where a bimodal distribution is apparent, the base of the lower peak in fossil abundance consistently occurs within a 3 cm stratigraphic range (10–13 cm above the NHC) whereas its top may vary in stratigraphic position by up to 6 cm (18–24 cm above the NHC). This consistency indicates that the lower assemblage is generally structured as a tabular bed of nearly uniform thickness in stratigraphic conformity with the structural plane of the NHC. This conclusion is further corroborated by our 3D analyses of fossil distributions. Slope and dip values each remain essentially constant near zero across the top and base of this assemblage, regardless of the quantitative data examined (Table 3). Our nearest neighbor analysis identified the stratigraphic interval from 10.6 to 19.1 cm above the NHC to also form a tabular region encompassing the highest density of oysters within our excavation (Figure 5). That interval overlaps stratigraphically with the lower peak in fossil abundance (Figures 2, 3), further supporting characterization of this interval as a distinct fossil layer. Cumulatively, these findings demonstrate that this lower layer is generally tabular in form and that it dips (relative to the NHC) by less than 1°.

Since this stratigraphically-lower fossil assemblage is predominately comprised of oysters and exhibits a correspondingly-low Simpson Diversity Index, we elect to describe it as an “oyster layer” near the base of the Hornerstown Formation. Allowing for single-centimeter intervals of stochastic variation in fossil composition, this oyster layer exhibits a Simpson Diversity Index <0.3 due to >80% of its fossil specimens pertaining to the oyster Pycnodonte. Gallagher (1991) reported a similar high percentage of 78% for representation of Pycnodonte oysters within the basal Hornerstown, with this slightly lower percentage likely owing to pooling of data from the oyster layer (as defined above) and overlying MFL. The striking consistency of percent Pycnodonte and SDI throughout the oyster layer make them the best, most-reliable parameters for determining the upper and lower boundaries of the oyster layer, which in turn constrain its overall thickness within our study area. Based on these metrics, the oyster layer is 9 cm thick and extends from 8 to 17 cm above the NHC. In terms of faunal composition, the most consistent interval within the oyster layer (i.e., the portion most likely to be recognized while in the field) occurs from 12 to 17 cm above the NHC, thus pertaining to the upper portion of the assemblage. This span continuously meets the above SDI and percent Pycnodonte criteria (i.e., without single-centimeter “violations”) and contains over 40 fossils per 1 cm stratigraphic bin within the excavation area. Fossil abundances, taxonomic richness, and proportional abundance of Pycnodonte all decrease gradually down-section toward the base of the oyster layer (i.e., from 12 cm down to 8 cm above the NHC; Figures 2, 3, Supplementary Data File S1).

Multiple taphonomic scenarios could account for the gradational base of the oyster layer. Because the genus Pycnodonte is an epifaunal recliner adapted to living on soft substrates (Dhondt, 1984; Macellari, 1988; Callapez et al., 2015), the oyster layer appears to represent an autochthonous assemblage. Therefore, the gradational base of this layer could reflect: 1) a rapid mass death event followed by minor disturbance of fossil distributions via intense bioturbation of these sediments by thalassinoid crustaceans (Wiest et al., 2016); 2) a temporally-building mortality event over the order of 36 kya (cf. Esmeray-Senlet et al., 2017, based on fitting a best-fit sediment mixing model to the iridium profile), ∼75 kya (cf. Gallagher, 1993, based on dinoflagellate biostratigraphy), or ∼90 kya (cf. Obasi et al., 2011, based on glauconite maturity levels across this interval), and/or; 3) that the oyster layer assemblage formed by more random attritional accumulation. The significant taphonomic implications of each of these non-mutually-exclusive scenarios should be considered in a future study. For the moment, however, we note that the rarity of taxa other than Pycnodonte within the oyster layer indicates potentially selective die-off of oysters in an ecological event, and this character clearly contrasts what is seen in the upper fossil assemblage (the MFL, see below).

In contrast to the oyster layer, the positions of the boundaries (and therefore thickness) of the upper peak in fossil abundance are more variable. Certain spatial subpartitions, namely the south half and southwest quadrant, exhibit a position for the upper peak of fossil abundance with its base at 23 cm above the NHC and top at 28 or 32 cm above the NHC. Other subpartitions exhibit slightly greater variation, with the base and top of this upper peak in fossil abundance occurring as high as 30 and 37 cm above the NHC (in the north half), respectively. Thus, depending on how the data are partitioned, the base of this upper peak in abundance may vary in stratigraphic position by up to 7 cm and its top may vary in position by up to 9 cm.

Presumably this upper fossil assemblage should correspond to the MFL. However, prior characterizations of the MFL relied heavily on attributes beyond fossil abundance (e.g., presence of diverse vertebrate remains; Gallagher, 1993; Gallagher, 2003), so it will be helpful to briefly review the origin, evolution, and varied uses of the name Main Fossiliferous Layer over the last 50 years before assigning stratigraphic bounds to the MFL herein. The phrase “main fossiliferous layer” was originally used by Parris (1974) as an informal term for the most prolific fossil-bearing assemblage in the lower Hornerstown Formation at Edelman Fossil Park, then an Inversand Company mining quarry. Gallagher (1984) informally synonymized Parris (1974) “main fossiliferous layer” with other informal names historically applied to this layer, such as the lower Hornerstown “bonebed” (White, 1972; Richards and Gallagher, 1974; Gaffney, 1975), while also capitalizing the name MFL for the first time; this began continued use of the term Main Fossiliferous Layer as a formal name for this prolific assemblage (e.g., Gallagher et al., 1986; Parris et al., 1986). However, to our knowledge, no study has ever used precise, quantitative methods to define the stratigraphic extent or faunal characteristics of the MFL. Therefore, below we propose unique, spatially-consistent trends within our data as defining attributes of the MFL in order to help standardize use of this name in future studies, which is direly needed given the frequent inconsistencies in application of this name in recent years (summarized in Table 1 and Supplementary Figure S1).

At least a portion of what most authors have previously called the Main Fossiliferous Layer clearly corresponds to the upper fossil assemblage identified herein. Based on the distribution of fossils in our excavation area and their taxonomic identities, we distinguish the MFL as solely pertaining to this upper assemblage. Specifically, throughout the entire interval from 23 to 33 cm above the NHC, at least one of the following criteria are met: fossil count >40 per 1 cm stratigraphic bin; SDI ≥0.5, and/or; taxonomic richness ≥7 (Figure 3). At least two of these metrics are met from 24 to 32 cm above the NHC, signifying that this interval represents the most consistent portion of the MFL and the most likely interval to be recognized while in the field. Fossil abundances and taxonomic richness each decrease rapidly below and above this 8 cm stratigraphic interval, indicating that the MFL exhibits relatively sharp basal and upper bounds. SDI values also decrease away from this interval, but at a lesser rate. Owing primarily to fossil occurrences within the north half of the excavation area, richness and SDI (but not fossil abundance) values for the entire square meter each increase again to levels characteristic of the MFL over a thin stratigraphic interval just above the MFL as defined above. This is particularly evident from 35 to 37 cm above the NHC, which corresponds to 3–5 cm above the MFL, wherein both taxonomic richness and SDI values meet the metrics used to define the MFL. We therefore term this thin interval a “shoulder” of the MFL (Figure 3).

To further explore this question, we used spatial analysis tools in ArcGIS Pro to examine patterns in the 3D occurrences of fossils with greater (0.1 cm) resolution. Our study is one of the first to use 3D GIS methods for the purpose of not only analyzing fossil assemblages in high-scale resolution, but of also defining their planimetric and stratigraphic extents (see Supplementary Section S1.5 for more on our digital methodological advancements). Several lines of evidence indicate the upward-sloping “shoulder” of the MFL is a genuine feature of this layer’s geometry within our excavation area rather than an artifact of our data processing methods. First, the concave hulls of several key groups which may be considered representative of the fauna of the MFL, such as tetrapod bones, plant remains (wood and Xylophagella boring casts), and Cucullaea clams, each display a dramatically-positive, northerly tilt, as does their combined intersection hull, created to represent the general form of the MFL in 3D (Figure 6). Concave hulls of other, more stratigraphically-distributed taxonomic categories such as fish bones also display northern tilts, but of gentler slopes (Supplementary Figure S2). Visual examination of the entire dataset in 3D (Figure 4) also reveals an upward tilt of diverse fossils toward the north half of our excavation area, which cumulatively create the “shoulder” of the MFL seen in our 2D analyses. Additionally, our dip calculations of the MFL from our 3D analyses are consistent with the understanding that the top of the MFL is more steep than its base, which is consistent with our 2D findings but inconsistent with the near-horizontal dip of the NHC and oyster layer (Table 4).

Such variations in stratigraphic position and thickness of the MFL suggest that it does not form a simple, tabular layer but instead exhibits small-scale undulations or a tilt that do not match the structural dip of the Hornerstown Formation in southern New Jersey (0.08–0.12°, Steckler et al., 1999; Esmeray-Senlet et al., 2015). In contrast, the oyster layer displays considerably less variation in its stratigraphic position and thickness. Thus, the two layers do not form planes parallel to one another within our excavation area, and the MFL in particular appears to exhibit a slight overall tilt relative to the NHC and oyster layer that cannot be accounted for by structural dip. This implies that the sedimentary and/or taphonomic processes which created unevenness within the MFL are temporally constrained to have occurred after deposition of the oyster layer. Although sediments at this site have been heavily bioturbated by decapod crustaceans (presumably cf. Hoploparia sp.; Gallagher, 2003), it is highly unlikely that bioturbation can account for the stratigraphic separation and contrasting diversity metrics of these fossil assemblages because: 1) crustaceans are documented to have been burrowing throughout this entire timespan (not just during deposition of the MFL; Wiest et al., 2016), and; 2) the vast majority of macrofossils collected in this study (e.g., Cucullaea steinkerns, Pycnodonte valves) are much larger than the average diameter of the Thalassinoides burrows within this interval (∼1.1 cm; Wiest et al., 2016). Therefore, other explanations must be examined. Indeed, few options remain that could explain these small (<10 cm of stratigraphic variation over 1 m of distance) but noticeable (during excavation) variations in stratigraphic geometry of the MFL.

One potential explanation for the complex geometry of the MFL is that we have identified natural variations in paleobathymetry over which the fossil organisms lived and died. Centimeter to multiple decimeter-scale variations in seafloor topography within a single square meter area are common on the shallow shelf owing to currents and the actions of bioturbators (e.g., Wang and Tang, 2009; Schönke et al., 2017), both of which would have also exerted influences on the seafloor topography at Edelman Fossil Park at the time of deposition of the MFL. Therefore, we hypothesize that the subtly-undulatory geometry of the MFL may have arisen from these processes, with bioturbation likely playing a major role (e.g., Hoploparia and/or other small crustaceans likely made spatially-heterogeneous mounds of sediments on the seafloor as they excavated Thalassinoides burrows; cf. Felder, 2001; Rodríguez-Tovar et al., 2010). Although this hypothesis requires confirmation from additional 3D mapping of fossil occurrences over a larger area, it could plausibly explain the comparatively even/uniform character of the southern half of the MFL and its gradual slope upward toward the north. An alternative hypothesis of a local east-west trending fault is not supported by our data because the MFL fossil data points form a gradual slope rather than a stepped offset (Figures 4, 6). Additionally, any such fault would have to have been an extremely small-scale feature (i.e., a “microfault”) as fossils within the underlying oyster layer do not exhibit a similar slope/step, nor does the NHC. Still, other alternative hypotheses remain plausible, such as attritional accumulation of organic remains over a short period of time (<250 kyr; Gallagher, 1993), resulting in them being preserved at multiple stratigraphic horizons within the MFL. The continuous, gentle nature of the slope in MFL fossil positions (Figures 4, 6), however, seems difficult to account for by attritional accumulation alone, which (by definition) would usually be expected to yield more erratically-positioned fossils. Repetition of this study over a larger geographic area would be required to fully evaluate any of these taphonomic hypotheses.

To examine the faunal profile of each fossil assemblage, we combined fossil occurrence data from the stratigraphic bins encompassed by each layer and then calculated an average taxonomic richness, Simpson Diversity Index, and percent Pycnodonte for each pooled dataset. This reanalysis found the 9 cm interval encompassing the oyster layer (8–17 cm above the NHC) to be composed of 87% Pycnodonte with an SDI of 0.25 and taxonomic richness of 12. In contrast, the 10 cm interval encompassing the MFL (23–33 cm above the NHC) contains a lower percentage of Pycnodonte (71%) and substantially higher SDI (0.49) and taxonomic richness (22) values. Thus, contrasts in taxonomic richness, diversity, and proportional abundance of the oyster Pycnodonte further enhance distinction between these two assemblages initially identified by peaks in our fossil abundance data.

Finally, although our thorough analyses have elucidated many contrasting features about these two fossil layers that will aid ongoing efforts to understand the formative history of the MFL and faunal changes occurring across the K/Pg boundary in the northern Atlantic Coastal Plain, several of these distinguishing characteristics of each layer would be difficult to identify in the field. We therefore propose the following “field definitions” of the oyster layer and MFL to facilitate their recognition in the field. The oyster layer begins roughly 10 cm above the NHC, is approximately 5–10 cm thick, and is marked by a great abundance of fossils, primarily Pycnodonte oysters with occasional specimens of other taxa. The MFL, by contrast, begins a little over 20 cm above the NHC contact, is roughly 10 cm thick on average, contains more fossils than the oyster layer (N = 420 versus 366, respectively; Supplementary Data File S1), and fossils within it pertain to a considerably more diverse suite of taxa, including plants (primarily Xylophagella wood boring casts) and large vertebrates (e.g., turtles, crocodilians, and mosasaurs). Specimens of Cucullaea vulgaris are likely to be encountered frequently within the MFL (Gallagher, 1991), but this taxon is not especially numerous nor uniformly distributed spatially; therefore, wood and tetrapod remains (i.e., bones and teeth) are also useful fossils for identifying the MFL while in the field. Further, Cucullaea vulgaris specimens can also be found both above and below the MFL, so recovery of a single individual of this bivalve is not enough to identify the position of the MFL. It is true, however, that within strata near the NHC, steinkerns occur most commonly in the MFL [Figure 2, and Gallagher (1993)]. The gap between the base of the MFL and the top of the oyster layer can vary from approximately 5–10 cm, which can be recognized during top-down excavation by a relative drop in the number and diversity of fossils recovered. These field definitions, based on the 1 m² area studied herein, are additionally supported by anecdotal observations from our larger systematic excavation of over 200 m² through the same stratigraphic interval at Edelman Fossil Park.

Comparing the stratigraphic extent of the MFL as defined herein to its position and thickness as described in previous reports (Table 1, Supplementary Figure S1) reveals that much of the variation among studies can likely be attributed to the use of varied stratigraphic thickness and position for this layer and the diverse methodologies that researchers have used to study it. For instance, studies which have reported the MFL to be 20 cm or more thick (Gallagher, 1993; Staron et al., 2001; Obasi et al., 2011) have clearly included the oyster layer within their assignment of the MFL. This is not surprising given how the two assemblages occur in very close stratigraphic proximity and together form an obvious, broad peak in fossil abundance compared to the strata overlying and underlying them which contain considerably fewer fossils (as also evidenced by our data; Figures 2, 3). However, our findings strongly support stratigraphic, paleoecologic, and taphonomic distinction between these two fossil assemblages. Consideration of the oyster layer as a basal portion of the MFL would also explain most differences in the reported placement of the base of the MFL by prior studies, as this would result in a lower base for the “assemblage” as a whole (as interpreted by Gallagher, 1993; Staron et al., 2001; Obasi et al., 2011). For example, the 30 cm thickness reported by Staron et al. (2001) suggests that they included the oyster layer within the MFL, and their report of the base of the MFL as being 10 cm above the NHC further supports this conclusion as that is roughly the location of the base of the oyster layer as defined herein. Therefore, we suggest that, stratigraphically, the majority of prior reports actually agree more than the reported values suggest, and the differences reported in Table 1 simply stem from varied assignment of the MFL.

Based on our data, we agree with Wiest et al. (2016) and Horner et al. (2016) that the MFL is quite thin and that small magnitude variations in thickness of the MFL exist (e.g., 10–20 cm; Horner et al., 2016). We add that this variation is likely due to the weakly-undulatory, non-tabular geometry of this layer identified herein. Additionally, this new, formal definition restricts the name Main Fossiliferous Layer to pertain specifically to the highly-diverse fossil assemblage in the lower Hornerstown Formation which most researchers have been interested in studying due to its potential connections with the K/Pg mass extinction. The only prior report whose claims we cannot reasonably reconcile with these conclusions is that by Esmeray-Senlet et al. (2017), who reported uniquely higher numbers for the base and thickness of the MFL (Table 1). We speculate that the uniqueness of their results might have arisen from uncertainties involved in assigning the placement of the MFL within the limited sample volume of a drill core.

Defining the stratigraphic boundaries and thickness of the MFL is important for characterizing both its formational history and chronostratigraphic relation to the K/Pg boundary. For example, assuming a relatively constant sedimentation rate [based on glauconite maturity (Obasi et al., 2011) and dinoflagellate biostratigraphy (Koch and Olsson, 1977)], the 10 cm thickness definition proposed herein suggests that the MFL was deposited over a shorter timeframe than inferred from many prior, thicker estimates (e.g., Gallagher, 1993; Staron et al., 2001; Obasi et al., 2011; Esmeray-Senlet et al., 2017), which in turn suggests that [in agreement with Boles (2016) and Wiest et al. (2016)] fossil remains in the MFL are less time averaged than some have thought. Additionally, when the definition of the MFL proposed herein is applied to prior reports (listed in Table 1) of the location of the K/Pg boundary, the chronostratigraphic location of this boundary changes from the base of the MFL in Gallagher (1993) and Obasi et al. (2011) to actually occur at the base of the oyster layer (Obasi et al., 2011) and in the middle of the oyster layer (Gallagher, 1993). Therefore, they no longer agree with the placement of the K/Pg boundary by Wiest et al. (2016) and Horner et al. (2016) which remain at the base of the MFL. Furthermore, using this definition of the MFL, Esmeray-Senlet et al. (2017) actually placed the K/Pg boundary well above the top of the MFL, in general agreement with Staron et al. (2001). Prior reports that placed the K/Pg boundary at the NHC (e.g., Minard et al., 1969; Self-Trail and Bybell, 1995) are not affected by any definition of the MFL. It is not possible to reconcile any other prior claims about the position of the K/Pg boundary (e.g.; by Olsson, 1963; Petters, 1976) into the precise centimeter-scale stratigraphic framework examined herein as they did not report a precise stratigraphic position nor thickness of the MFL. Nevertheless, the discrepancies among the studies reconciled above emphasize the need for all future researchers to employ the same definition of the MFL. Universal use of the same definition will allow any precise geochronologic indicators (e.g., shocked quartz, spherules) or other signals of environmental disturbance found in this interval in the future to be placed in their exact chronostratigraphic positions relative to this fossil assemblage. Only then can the origin of the MFL and any connection it may have with the Chicxulub impact (Obasi et al., 2011) truly be discerned with certainty.

Since herein we did not search for nor collect any precise geochronologic indicators, we can only attempt to synthesize prior claims about the position of the K/Pg boundary. Prior researchers have used a wide variety of methods to discern the position of this boundary (i.e., geochemistry, ichnology, biostratigraphy), and few studies have employed the same techniques. For example, among the studies reconciled in the previous paragraph, only Esmeray-Senlet et al. (2017) made their assignment based directly on a precise indicator of the Chicxulub impact event (an iridium spike). Wiest et al. (2016) alternatively placed the K/Pg boundary based on an abrupt reduction in the average diameter of Thalassinoides burrows, which Horner et al. (2016) found to be accompanied by a shift in the elemental composition of the sediments infilling burrows. One possible hypothesis to reconcile these data and the associated contrasts in placement of the K/Pg boundary among studies could be the following: given that 1) iridium-rich dust was likely settling through the water column to the seafloor for a period of several years [Vellekoop et al. (2016), and references therein] after the moment of impact [which defines the geochronologic K/Pg boundary; Molina et al. (2006)] and 2) iridium can be remobilized in unconsolidated sediments [e.g., Hull et al. (2011), Racki et al. (2011), Wallace et al. (1990)] whereas burrows cannot, it is plausible that iridium-rich dust predominately fell after decapod crustaceans were already constrained to smaller body sizes by post-impact environmental disturbances. This scenario would be consistent with a hypothesis that the K/Pg boundary occurs at the base of the MFL, in agreement with the conclusions of Wiest et al. (2016), Horner et al. (2016).