Time-dependent taphonomic site loss leads to spatial averaging: implications for archaeological cultures

The study addresses how time-dependent, taphonomic loss of archaeological sites affects the spatial patterns archaeologists infer from material culture similarity. Cultural areas are often defined by perceived similarities across sites, yet the archaeological record is inherently incomplete and spatially sparse, especially in older periods. The authors hypothesize that spatially incomplete datasets (due to site loss and discovery bias) cause spatial averaging, whereby regions of similarity appear larger than they originally were. To avoid cultural-process assumptions, they test this using spatially autocorrelated, culturally independent soil data to examine how subsampling (mimicking site loss) alters detected distances to the first local minimum of similarity, which archaeologists often implicitly use to demarcate cultural zones. The work is important because very large prehistoric cultural regions (e.g., Aurignacian, Gravettian, Acheulian) may be partly artifacts of spatial averaging rather than solely sociocultural processes.

The paper situates its question within debates on cultural classification and the use of similarity to define archaeological cultures. Prior work critiques mapping material culture directly onto discrete ethnic or cultural groups and highlights the polythetic, cross-cutting nature of material culture traits. Quantitative approaches have examined distance-decay in cultural similarity, drawing on Tobler’s First Law of Geography, and similar distance-decay patterns are well known in ecology and biogeography. Recent archaeological studies emphasize how temporal processes (time averaging, taphonomic loss) distort behavioral signals and can inflate apparent spatial spread of cultural signatures. However, the spatial dimension of data loss (spatial averaging) has received comparatively less attention. This study builds on that gap, testing how reduced spatial sampling inflates inferred regions of similarity.

Data: The European Soil Database v2.0 (ESDAC) comprising polygons with soil classification attributes across Europe and Russia was converted to point data using polygon centroids in ArcGIS 10.4.1. Two spatial extents were analyzed: continental Europe (~25,000 points) and a national subset (Germany, ~2,500 points) to mitigate classification inconsistencies across countries. Thirteen categorical variables were selected (e.g., WRBFU, slope, depth to rock, accumulated temperature, surface/subsurface texture, packing density, structure, mineralogy, parent material hydrogeological type) to emulate trait-based comparisons analogous to archaeological practice. Similarity and distance: Geographic distances were computed in kilometers. Pairwise similarity between a randomly chosen reference point and other points was computed via simple matching across the 13 ordered categorical attributes as the proportion of matching categories (range 0–1). Decay curve and first minimum: For each reference point and sample, similarity was modeled as a smooth function of distance using a generalized additive model (GAM) in R 3.6.1. The first local minimum of the similarity–distance curve was identified as the operational boundary of a similarity region. Subsampling to mimic spatial incompleteness: Random subsets of the total points were drawn at 90%, 80%, …, 10% and then 9%, 8%, …, 1%. For each proportion, the procedure (random reference, GAM fit, first minimum extraction) was repeated 100 times to capture variability due to sampling and starting-point location. The absolute distance (km) to the first minimum was recorded for each run. Taphonomic “aging” scenarios: To approximate time-dependent site loss, an exponential decay was fit to a database of >14,000 reliable radiometric dates for European Palaeolithic sites (7000–900,000 BP). Resulting proportions of surviving sites for modeled ages 10 ka, 50 ka, and 100 ka BP were used to set sample sizes (491, 96, and 12 points, respectively) in the soil dataset, and distances to first minima were compared across ages using Mann–Whitney U-tests. Modeling relationships: For the Europe-wide data, nonlinear least squares was used to fit a negative power function d = a P^β to the relationship between distance to first minimum (d) and percent of points (P), and linear regressions were also fit on log-transformed distances in two ranges (100–10% and 10–1%). For the Germany subset, a log-linear model was found to best fit the data. All code and scripts were provided via GitHub.

• Across the European dataset, there is a significant negative power relationship between the percentage of sampled points and the absolute distance to the first minimum in similarity: as fewer points are sampled, the inferred distance to the first minimum increases, implying larger regions of similarity. Parameters for the power fit were significant at p < 0.001; β was estimated at −0.04148 with 95% CI [−0.052, −0.031]. Residual standard error was 513.1 km, reflecting variability due to different starting points.
• Linear regressions on logged distances confirm the pattern with negative slopes in both ranges: 100–10% slope ≈ −0.0005; 10–1% slope ≈ −0.001 (nearly twice as steep), indicating stronger effects at lower sampling proportions.
• Within-country analysis (Germany) shows a similar negative relationship; a log-linear model fits best. For Log(distance) as the dependent variable, the coefficient on percentage of points is −0.003 (SE 0.0004), p < 0.05; R² ≈ 0.036; n = 1,354.
• Modeled taphonomic loss (“aged” samples): Distances to the first minimum increase with modeled age (10 ka vs 50 ka vs 100 ka). Mann–Whitney tests show significant differences (p < 0.01) between 10 vs 50 ka and 10 vs 100 ka, but not between 50 vs 100 ka, suggesting diminishing returns at extremely low proportions (~0.004% and ~0.0005% of points).
• Starting-point effects: While the overall trend holds regardless of starting point, the size of the inferred similarity region can differ substantially with starting location; at lowest proportions, regions can be nearly twice as large depending on where the reference point is placed.
• On continental scales, each 1% reduction in sampled points can increase the distance to the first minimum by roughly 0.5–3%, producing large overestimates of similarity areas over thousands of kilometers.

Findings demonstrate that spatial averaging due to incomplete spatial sampling inflates the apparent extent of similarity in spatially autocorrelated datasets. For archaeology, this means that cultural areas inferred from sparse site distributions are prone to overestimation, especially at low sampling densities and in older periods where taphonomic loss is severe. Although explained variance (R²) is small because similarity varies widely at given distances and depends on reference-point location and mosaic-like spatial structure, the negative relationships are statistically robust. This indicates that the proportion of available points significantly affects the detected boundary (first minimum) of similarity regions. The effect intensifies at the smallest sampling proportions, making large-scale cultural attributions particularly vulnerable to bias. The study underscores that observed large cultural regions may partly reflect spatial averaging rather than purely sociocultural processes, and highlights the need to consider spatial and temporal information loss in interpreting material culture similarity and cultural taxonomy.

The study shows that spatial averaging increases as data become sparser, inflating inferred regions of similarity, and that the location of the initial reference point influences the estimated size of those regions. Consequently, older archaeological periods and regions with greater geomorphic disturbance likely suffer stronger spatial averaging, leading to overestimation of cultural areas—often compounded where lithics are the primary data. Because similarity-based cultural areas may not mirror original cultural geographies, archaeologists should be cautious in using their size to infer population and cultural dynamics. Future research should quantify the magnitude of these effects with archaeologically relevant datasets and examine the combined impacts of spatial and temporal averaging on cultural patterning.

• The analysis uses soil classifications as a culturally independent proxy; while appropriate for testing spatial autocorrelation and subsampling effects, magnitudes of effects may differ for archaeological datasets.
• Model fits exhibit low R² and high residual variance due to mosaic spatial structure and dependence on starting-point location; results are not intended to predict exact effect sizes in specific cases.
• International inconsistencies in soil classification warranted a national-scale test (Germany); broader generalizations across varying classification systems may introduce noise.
• The parameter estimates in the negative power model are sensitive to sampling variability and the chosen curve-fitting approach (GAM for minima detection); different similarity metrics or modeling choices could alter estimates, though the qualitative trend should persist.