Exploring the Role of Climate and Other Environmental Factors in the Distribution of Paleolithic Sites in France Between 34 and 21 cal ka BP: A GIS Analysis

1. Introduction

Since the pioneering work of Demars [1], the role of climate variability on the settlement of Western Europe during the late Paleolithic has been the subject of a great deal of research in recent decades (e.g., [2,3,4,5,6]). According to data on sub-actual hunter-gatherers, population density is closely linked to the mode of subsistence, which is itself strongly influenced by the environment and its productivity in terms of vegetation and animal biomass [7,8]. It follows that climate, on which the vegetal environment largely depends, would have influenced Paleolithic settlement, by modulating the quantity of available resources and the extent of the various biotopes exploited over time.

For the late Paleolithic, all studies based on the analysis of georeferenced databases of archaeological sites agree that there was a contraction of settlement into lower latitudes of Europe during the coldest periods of the Last Glacial Maximum (LGM), ca. 26–19 calibrated kiloanni before present (cal ka BP) [9]. A recolonization of northern latitudes and mountain ranges during Europe’s deglaciation followed this geographic contraction [1,10,11,12,13]. During cold periods, the most densely populated areas were Southwestern France, the southern European peninsulas and parts of Central Europe. At the same time, the frequency of sites (or dated sites) per period, considered by several authors to be an indicator of population density, decreased in Europe during cold phases, reaching a minimum during the Last Glacial Maximum [1,3,14,15,16,17]. This drop in density correlates with a genetic bottleneck in Western Europe and a genetic turnover in Southeastern Europe [18,19].

From the earliest studies [1,10], the decline in population has been interpreted as the result of a combination of lower temperatures and drier conditions. Without completely revising this view, subsequent studies have made extensive use of climate simulations and have made it possible to test the respective roles of various parameters. Modelling settlement by taking into account known relationships between climate and population density for hunter-gatherers, and incorporating evapotranspiration (ET), the hydrological budget (a proxy for ecosystem productivity), and the temperature of the coldest month (a factor influencing winter mortality), correctly reproduces the distribution of LGM sites in Europe [3]. Using new climate simulations, Burke et al. [4] indicate that the best explanatory variables for site distribution are altitude and climatic parameters, particularly spring rainfall variability, average autumn temperature, and temperature variability in autumn and winter. These factors, particularly temperature and precipitation variability, have been interpreted in terms of ecological risk that human groups faced [20], in which higher ecological risk reduced the predictability of animal and plant resources. Such factors would, thus, play an important role in population size fluctuations.

A number of studies point to the discontinuous nature of areas favourable to Paleolithic settlement in Europe during the Last Glacial. The ways in which settlement patterns were adjusted in response to the coldest periods remain a matter of debate (migration, local extinction of groups, persistence of cryptorefugia favouring rapid recolonization of northern regions during climatic improvements, changes in the cycles and routes of nomadism) [4,17,21]. Others, however, suggest that there was no fragmentation of settlement areas, allowing for an uninterrupted flow of cultural elements and genes between populations [3,22]. Eco-cultural niche modelling of the Solutrean by [23] also suggests that much of southern Europe was contained within its reconstructed niche, but these groups did not occupy it in its entirety due to the presence of another contemporary technocomplex situated in its eastern distribution (the Epigravettian).

The following assumption tends to form the basis of most, if not all, studies: the spatial distribution of known Paleolithic sites and their number provide a reliable indication of populated (or, conversely, unpopulated) areas and the density of occupation. Although known sites represent only a small fraction of the occupations that existed, it is often assumed that they provide a representative view of past settlement. Formal tests of this assumption are lacking. Some prehistorians who have empirical experience of various technocomplexes, however, argue that this premise needs to be rigorously evaluated and discussed. Intensive terrestrial surveys at a regional scale always lead to the discovery of new sites, some of which represent the first observations of a particular archaeological culture in those areas. To a certain extent, however, this assumption of representativity has been validated a posteriori at the European scale: the general pattern of populated areas across Europe has not varied significantly between the first published estimates [1,14] and those proposed subsequently [3,4,13,17,21,24,25], despite the increase in the number of known sites over the years. This situation suggests that other potentially relevant factors, such as bias linked to regional differences in knowledge or issues relating to site preservation or accessibility, do not appear to play an influential role. This hypothesis, however, warrants further research. In particular, the presence of thick loess cover in Northern France has been cited as a factor limiting the visibility and accessibility of sites, thereby explaining the low number of known occurrences in this region (e.g., [26,27]). A few sites have indeed been discovered in loess sequences several metres deep during exceptional works (e.g., [27,28]).

Although the prevailing hypothesis considers climate to be the primary factor responsible for the settlement pattern of Late Paleolithic Europe, the various studies suffer from limitations. These shortcomings are:

(1) The employed climate models have their own limitations. In particular, simulations often insufficiently reproduce temperatures derived from proxies such as pollen [29] or field evidence of permafrost [30]. The optimal conditions of temperature and precipitation for settlement scenarios proposed in some studies are, therefore, not reliable, and the boundaries of large biomes remain imprecise. However, these problems have been partially resolved in recent years with advances in modelling (e.g., [31]). Millennial scale climate variability (i.e., the Dansgaard-Oeschger events) is poorly accounted for in many simulations.

(2) Non-linear effects between climate and population density may exist, particularly in relation to vegetation cover dynamics. Tallaavaara and Seppä [32] showed that, in Scandinavia, the gradual colonisation of the forest by spruce during the first half of the Holocene led to a significant decline in population density before the arrival of agriculture, as boreal forests have lower biological productivity than deciduous forests. It remains difficult to reconstruct vegetation in detail during the Last Glacial, both from proxies (pollen, macroremains) and from simulations. However, vegetation studies converge in considering that steppe-tundra with a few patches of boreal forest in valleys remained the dominant biome throughout the Late Paleolithic in mid-latitude Europe (e.g., [33,34]). The effects related to vegetation cover dynamics were likely to have been minor and, therefore, we consider them to have been negligible.

(3) Environmental factors not directly related to climate but determined by geological heritage and relief may have played an important role in settlement. This is particularly true of sand deserts, which developed extensively during the Last Glacial in areas without relief and where abundant sources of sand (particularly coastal and glacial) were available for deflation. Bertran et al. [35] demonstrated that the cover sands that developed in Southwest France provide a convincing explanation for the distribution of sites observed at a regional level. According to these authors, the sand deposits, which have yielded few traces of Late Paleolithic settlement, were low-productivity areas unattractive to animals and their predators (including humans) and likely acted as a physical obstacle modulating exchanges between settlement areas. This could, in part, explain the geographically differentiated cultural characteristics for some technocomplexes in Southwest France [35,36,37]. Other sand deserts existed in Europe, and they too may have played a role underestimated until now.

(4) Other factors favourable to settlement are mentioned frequently in regional studies, such as the proximity of lithic raw material (flint) sources [38,39] and relief, in particular large valleys (e.g., [35,40]). Although the transport of raw materials over long distances (over 150 km) has been documented in the late Paleolithic contexts (e.g., [41,42,43,44]), the proximity of sites to known source areas is often emphasised. In northern Hungary, Trajer [39] used statistical analysis of late Paleolithic occupations to show that the average distance between sites and raw material sources is less than 10 km. In Southwest France, preventive archaeology also has noted a wealth of late Paleolithic occupations near outcrops of high-quality flint such as the Bergeracois flint and conversely, the scarcity of occurrences on tertiary molasse substrate lacking flint. The influence of these factors on the distribution of sites at a larger scale seems likely and warrants further investigation.

5) While sociocultural dynamics and group organization were influenced to varying extents by climate variability [45,46,47], cultural choices and cultural evolution clearly played key roles in such dynamics. Consequently, spatial distributions and site frequency may have fluctuated due to seasonal activity cycles, mobility pattern, and group size (e.g., [48,49,50]). It is difficult to know, however, the extent to which such factors influenced occurrence distributions and whether such influences were consistent through time.

This study aims to analyse in detail the influence of various environmental factors sensu lato (geographical coordinates, topography, geomorphology, geology, current land use) on the distribution of Paleolithic sites, using an approach similar to that proposed by [51]. The period under consideration, 34–21 ka includes some of the coldest episodes of the Last Glacial, the GS-6 to GS-2.1a stadials, as well as more temperate interstadial intervals, GI-6 to GI-2 [52]. This chronological interval includes the Gravettian technocomplex, the development of the Solutrean, and the emergence of the Badegoulian.

The geographical scope of the study is restricted to present-day France in order to limit the time spent on a critical analysis occurrence contexts and lithic assemblage attributes. We compiled a database based on recent studies [36,50,53,54,55,56], as well as new data from preventive archaeology work. In practice, references were screened for precise information (such as descriptions and illustrations) that would allow for reliable cultural attributions: fossiles directeurs, technical markers and diagnostic chaînes opératoires, evaluations of assemblages’ archaeological or stratigraphic contexts. The various layers of information, which include data on site location and paleoenvironmental data, were cross-referenced using a Geographic Information System (GIS). The compiled data were analysed statistically to identify robust associations. Finally, we assessed and interpreted the relative influence of the various factors in order to propose hypotheses concerning observed archaeological distributions.

2. Methods

2.1. Database of Occurrences

In this study, Paleolithic technocomplexes are used as basic, temporally specific cultural entities for which variations in site distributions through time are examined. The term “technocomplex” describes a practical classification of industries defined by prehistorians and based on typological and technological similarities between assemblages [57]. Because archaeological technocomplexes are defined by the presence of discrete, temporally-diagnostic artefact types, they can be reliable means of attributing archaeological assemblages, both dated and undated, to specific past “cultural” entities. ¹⁴C-dated ornamented caves that have yielded either undiagnostic or no archaeological assemblages allowing for an attribution to a specific technocomplex to be made were excluded from the inventory.

Recent Bayesian statistical analyses of available radiometric ages, reliant on both strict archaeological associations between technocomplexes and archaeological contexts (i.e., stratigraphic constraints [37,58]), have made it possible to accurately place these technocomplexes chronologically and in turn correlate them to recorded climate variability. Based on currently available data, we consider that across our region of study, all occupations belonging to a given technocomplex are approximately contemporary and fall within the intervals established by [37]. In other words, we consider negligible any delay linked to the gradual expansion or contraction of a technocomplex across the territory due to the chronological resolution that characterizes our period of study.

The database used in this study is a relational database that combines different tables associated with geographical data and cultural attributions, radiometric ages, and raw material and technological data. From this relational database, we extracted the data presented in Table 1. In this table, each row corresponds to an archaeological occurrence. Each occurrence corresponds to one or multiple, superimposed archaeological levels. In the absence of refitting tests between layers at many sites excavated in the past, it is impossible to determine whether the levels distinguished in the field are behaviourally relevant and accurately reflect distinct occupations, a single occupation whose objects have been dispersed by taphonomic processes, or a palimpsest of multiple occupations (e.g., [59,60]). The occurrences selected are for which we could make a reliable attribution to one of the technocomplexes concerned by a study period based on lithic typology or technology. A total of 440 occurrences are retained (Table 1).

Despite critical scrutiny of the available data, our approach has potential biases and limitations that should be kept in mind:

• It is impossible to directly verify each occurrence and its associated assemblages. We, therefore, relied largely on published sources and grey literature reports, which can be variable in the details provided.
• Excavation, post-excavation and curation methodological biases, especially across the 19th and 20th centuries, may have obliterated certain diagnostic features necessary for making a reliable attribution to a technocomplex (e.g., [56]).
• The definition of technocomplexes has evolved over time (particularly with regard to different analytical approaches) and research efforts dedicated to defining and examining technocomplexes are unequal.
• Finally, the varying “detectability” of technocomplexes, either bibliographic in nature or due to difficulty in identification, may have influenced our dataset, as some technocomplexes present marked and easily recognizable diagnostic artefact types, whereas others possess characteristic artefacts that are present in other technocomplexes as well, making them harder to recognize (e.g., Late Gravettian, Badegoulian without “raclettes”).

For the purposes of analysis, the technocomplexes were grouped together to divide occupations into periods of roughly comparable duration and thus obtain statistically meaningful sample sizes. These groupings are as follows: Early and Late Badegoulian (= Badegoulian), Middle and Upper Solutrean (= Late Solutrean), Protosolutrean and Lower Solutrean (= Early Solutrean), and Recent and Upper Gravettian (= Late Gravettian). The duration of these periods, which correspond to the 95% confidence intervals calculated by [37], is 2,300 years (Badegoulian), 1,604 years (Late Solutrean), 2,271 years (Early Solutrean), 1,304 years (Final Gravettian), 2,495 years (Late Gravettian), and 2,900 years (Middle Gravettian). The Lower Gravettian has not been the subject of chronological examinations comparable to those carried out on other technocomplexes, therefore, its duration remains imprecise.

2.2. Preventive Archaeology Database (Caviar)

To evaluate whether uneven survey coverage produced a biased distribution of known sites in France, we consulted the Institut National de Recherches Archéologiques Préventives’ (Inrap) Caviar database. It lists areas surveyed and provides information on the extent of the surfaces investigated (trenches) and the number of stratigraphic units (SU) containing Paleolithic remains. This database is not homogeneous across the entire territory, as some regions have been investigated more intensively and data records began being maintained later in some regions than others. Nevertheless, it provides a picture of discoveries made through “blind” survey and, given that the methods used for data collection and analysis are consistent throughout the country, does not show any major bias in favour of any particular region. Due to the nature of preventive archaeology, areas close to large urban centres and transport routes are better surveyed than more remote areas. However, the most thoroughly studied areas are spread across all regions of the study area.

Using this database (last accessed in November 2024), the ratio between the number of stratigraphic units yielding Paleolithic artefacts and exposed surface area was calculated. Calculations were carried out according to the major administrative regions (ARA: Auvergne-Rhône-Alpes, BFC: Bourgogne-Franche Comté, CIF: Centre-Ile de France, GE: Grand Est, GO: Grand Ouest, HdF: Hauts de France, NAOM: Nouvelle Aquitaine, MM: Midi-Méditerranée). As late Paleolithic discoveries remain few in number, the entire Paleolithic period was included in the analysis.

2.3. Environmental Information Layers

Various layers of environmental information (topography, geomorphology, geology, current land use) were included and processed using QGIS software version 3.40.6. For the relief data, we used the European Digital Elevation Model (EU-DEM) with a horizontal resolution of 30 m (https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1). Rivers come from the shapefile provided by the Institut Géographique National (https://www.ign.fr/carto). We considered only rivers with a channel width greater than 15 m, as they are associated with large valleys.

For aeolian deposits, we employed two classes of data. The distribution of sand and loess proposed by [61], derived from a map of the topsoil texture produced as part of the European LUCAS project [62], takes into account all aeolian deposits, even those of low thickness. The distribution of loess proposed by [63], extracted from the 1:50,000 geological map of France, is complementary and only depicts areas where deposits exceed approximately 2 m in thickness and constitute extensive cartographic units.

We obtained the distribution of raw material sources (flint and related knappable stone) from the mapping efforts provided by a consortium of researchers, which is available as a shapefile at https://www.cartosilex.fr/ (last accessed November 2024). Since only a portion of France is covered, we produced a new map for the entire territory by extracting geological formations that contain flint in a primary position (mainly limestone) from the 1:1,000,000 geological map of France (https://infoterre.brgm.fr/) and verified the consistency between the two maps in areas where they overlapped. Due to contrasting resolution and the grouping of geological strata, the new map is less accurate and overestimates the extent of source areas. We also took into account formations containing flint in a secondary position, particularly weathering products of limestone. Only alluvial formations associated with rivers that cross primary or secondary flint sources over a significant portion of their course were included. We removed areas with thick loess cover [63] from the broad areas of weathering clay on chalk depicted on the 1:1,000,000 map in the northern half of France.

For climatic factors, we chose various indirect indicators: (1) altitude and latitude measured using the EU-DEM, (2) exposure (considered significant when the slope exceeds 5°), also measured using the EU-DEM, and (3) the extent of permafrost during the LGM. Due to the punctual nature of the permafrost indices recorded in the field [30] and the difficulty in precisely delimiting its extent, we used the output of a simulation carried out using the Max Planck Institute’s Global Climate Model ([31]; updated in [64]). The results show good agreement with field evidence and we, therefore, considered them robust. Permafrost indicates areas with null or negative average ground temperatures during the LGM. This distribution may have changed throughout the period under consideration (34–21 cal ka BP), particularly during long interstadials, but it is not possible to identify specific variations due to permafrost traces’ general lack of precise chronological data.

We derived the current land use type (in 2012) in France from the European Environment Agency’s CORINE Land Cover (CLC) map(https://www.statistiques.developpement-durable.gouv.fr/corine-land-cover-0), which classifies land into different categories (urban areas, agriculture, forest, etc.) based on satellite imagery. We used these data to identify the contexts in which recorded sites are currently situated.

2.4. Statistical Analyses

For the statistical analysis, we assumed that each occurrence of a technocomplex at a site was independent of previous occurrences at the same site. Thus, a site can be counted several times when if it contains multiple levels representing different technocomplexes. In such situations, we consider that the different occupations were sufficiently distant in time such that the choice of play by human groups of a given technocomplex did not play a role in the subsequent occupation of the site by human groups of a different technocomplex.

Our statistical assessments evaluated whether observed distributions deviated from random. To do so, we first defined the “credible” settlement area to be the portion of the landscape with an elevation below 600 m (see below). Next, we defined n to be the number of sites in the database. We repeated (B = 200) a procedure of randomly placing n locations (points) within the credible settlement area. We, thus, generated B distributions for the variable of interest (e.g., distance to the closest river, latitude, elevation…). We then compared these B random distributions to the observed distribution of archaeological occurrences (i.e., sites). We plotted the empirical cumulative distribution function of each of these distributions, along with the 95% confidence region of the observed distribution [65]. This allowed us to evaluate whether observed settlement patterns significantly deviated from random.

We applied a similar protocol for examinations of aeolian deposits and their relationship to archaeological occurrences when conducting a correspondence analysis that included distributions of the different types of aeolian deposits (observed distribution and B random distributions) as inputs. Observing the position of the actual distribution relative to the 95% confidence ellipse of the B artificial distributions allowed us to determine whether the association between occurrences and aeolian deposits significantly deviated from random expectations.

We employed R 4.6.0 [66] for all statistical analyses. The R codes produced for this study, along with the list of employed R packages, additional results and detailed elements of interpretation, are available on Software Heritage (https://archive.softwareheritage.org/swh:1:dir:bdb7416341238d5e0b372bc7287861a090d500f9).

3. Results

3.1. Distribution of Occurrences Listed in the Database

Overall, occurrences (all technocomplexes considered) are distributed unevenly across the territory (Figure 1) and show a concentration in Southwest France that extends into adjacent areas of the Iberian Peninsula (Cantabria, not documented here), as has already been noted in previous studies [1,4,13]. Areas of lower site density also appear around the Mediterranean Sea and, to a lesser extent, in the Paris Basin, while scattered occurrences are present in Northeast France and neighbouring regions (Belgium, Germany). For the broad technocomplex groups (see above), two types of distributions are apparent in the maps and the graph showing mean latitude as a function of longitude for each technocomplex (Figure 2). The first pertains to the Lower and Late Gravettian, during which a large part of the territory and neighbouring regions in the northeast were occupied. The other corresponds to the other technocomplexes, whose settlement contracted southwards, with only one observed site in the HdF and GE regions, i.e., of the northern and northeastern portions of the study area. The Late Gravettian has the highest average latitude of all the considered technocomplexes.

The frequency of occurrences by our technocomplex groupings exhibits significant variability, with a minimum during the Early Solutrean and Final Gravettian (Figure 3A). When weighting the number of occurrences by the duration of each technocomplex (taken from [37]), this significant variability through time remains. The pattern that emerges is that the number of occurrences per 1,000 years for technocomplexes that include an interstadial is significantly higher than those that were present only during a stadial (Figure 3B). The Late Gravettian, which spans a period including the GI-3 interstadial and possibly GI-4 (two brief periods of climate improvement following the extremely cold GS-5 stadial), stands in an intermediate position. This indicates that warmer periods are characterized by a higher frequency of archaeological sites. This is consistent with field observations in the loess deposits of Northern France where the relatively rare archaeological occurrences are associated with paleosols, which formed during interstadial events (e.g., the site of Renancourt [67]).

3.2. Caviar Database

The analysis of the Caviar database highlights the following points:

(1) Open surfaces (trench areas) vary greatly in extent by region (Figure 4A). They exceed at least 1.2 km² (HdF) and reach up to 17.0 km² (GE). This variability is due to the variable inclusion of such data, as many older operations are not included in the database.

(2) The ratio between the number of stratigraphic units containing Paleolithic remains (all periods combined) and investigated surface area (in km²) varies considerably by region (Figure 4B) and ranges from 0.833 (HdF) to 60.982 (MM). In general, the ratio decreases with latitude. The regions that have yielded the most remains are MM, NAOM, and ARA, i.e., the southern half of France. Conversely, the region that yielded the least is HdF, which is also the northernmost. The very high value for the MM region is likely an artefact of the criteria used to define and record stratigraphic units, especially when compared to NAOM. In the MM study area, some operations defined a stratigraphic unit on the basis of at least one observed Paleolithic artefact. Conversely, in the NAOM study area, stratigraphic units were created only when a large number of artefacts were observed, with smaller numbers or individual finds considered to represent “background noise”.

The picture that emerges from this analysis is that the distribution of Paleolithic sites found during large-scale archaeological surveys is not homogeneous, but rather concentrated in the southern regions of the study area. This distribution of archaeological occurrences from the Inrap database is similar to that visible in our database.

3.3. Distribution of Occurrences According to Land Use Type

A comparison of our database with the current land cover (CLC) shows that occurrences are mainly found in cultivated areas (45.6%), closely followed by forested areas (43.8%) and, to a minor extent in artificial areas (9.8%), particularly peri-urban areas (discontinuous urban fabric) (Figure 5). When these types of zones are weighted according to their relative surface area in France (< 600 m), discoveries are overrepresented in peri-urban areas, vineyards, and deciduous and mixed forests, while underrepresented in cultivated areas and, even more so, in pastures. Mineral extraction sites (quarries) are also overrepresented, but their overall importance remains minimal.

The distribution of site types (caves and rockshelters versus open air) indicates that 55.1% of the former are located in forested areas, while the latter are primarily observed in cultivated areas (including vineyards) and, more rarely, in forests (24.7%).

3.4. Distribution of Occurrences and Climate-Related Variables

The analysis of the distribution of occurrences according to altitude shows an overrepresentation of sites between 60 and 180 m, with a peak between 80 and 140 m (Figure 6A). A single site (Le Rond-du-Barry, Badegoulian, 768 m) has an altitude above 600 m, which empirically represents a “limit of habitability” with respect to sustained occupation. The comparison between the observed distribution and that of 200 random distributions of the same number of points in areas below 600 m indicates that the observed distribution deviates significantly from random (Figure 6B). None of the 200 random distributions falls within the confidence interval of the observed site distribution.

The analysis of the distribution of occurrences according to latitude also reveals a concentration between 43 and 46°N, and especially between 44 and 45°N (Figure 7A). The pattern, like with altitude, differs significantly from random expectations (Figure 7B).

Since local climate correlates strongly with altitude and latitude, it is probable that the factors behind these patterns are climatic in nature. Figure 8 shows the relationship between latitude and altitude of occurrences. As expected if climate played an important role, the maximum altitude of occurrences decreases as latitude increases. The distribution of occurrences based on site type (open air, rockshelter, cave) according to latitude and altitude does not reveal any clear trend, except for the fact that many cave sites are located between latitudes 43°N and 45°N, i.e., in the limestone-dominated landscapes of Southwestern France. Rockshelters and cave sites are more rare above 47°N, due to the reduced extent of karstic limestone formations in which they naturally occur. For all regions, rockshelters and cave sites are not the highest with respect to altitude.

We also tested the potential influence of site exposure using data from the EU-DEM. Orientations were grouped into classes of 20°. When only sites with a slope greater than 5° are considered, south-southwest orientations between 160°N and 260°N are predominant, with a peak at 180°N (southerly exposure) (Figure 9). In the northern hemisphere, such an orientation ensures a maximum exposure to solar energy at midday.

Figure 10 illustrates the distribution of occurrences, classified by technocomplexes, with the limits of permafrost during the LGM. These limits are based on simulations and field data and taken from [64]. Overall, permafrost does not appear to have limited land occupation, consistent with previous observations (e.g., [13]). However, if we consider the technocomplexes individually, it appears that the geographic extensions of the Middle Gravettian, Final Gravettian, Solutrean, and Badegoulian are limited to the northeast by the boundary of the continuous permafrost zone (i.e., by an air temperature of -3 to -5 °C). We do not observe this boundary for the Lower and Late Gravettian, for which sites exist in Northern France and neighbouring countries in the continuous permafrost zone, particularly in Belgium and Germany.

3.5. Comparison with Geological and Geomorphological Data

Overall, considering all periods together, very few occurrences are located in areas of aeolian deposits as mapped by [61] (Figure 11). Since the deposits cover only a small portion of France, the observed distribution must be compared to what it would be if the sites were randomly distributed across the territory, at least in the area identified as favourable for Paleolithic settlement, i.e., at elevations below 600 m (Figure 12A, B). This distribution depends on the proportion of surface area occupied by the different types of aeolian deposits. A comparison of the two bar plots clearly shows that occurrences in sand and loess are underrepresented. To evaluate whether this pattern deviates from random, we produced 200 random samples of the same number of points in France and evaluated the locations of the points in each sample with respect to aeolian deposits. We performed a correspondence factor analysis based on the results of these evaluations (Figure 13). It highlights the opposition between the distribution of occurrences and that of aeolian deposits. In other words, the former deviates from random expectations and reflects an inverse relationship between archaeological occurrence and aeolian deposits. The loess map proposed by BRGM [63] only covers areas with thick loess (> 2 m) and is therefore more limited in extent than in the previous map. Thick loess is mainly found in the northern third of France. The analysis produces a result consistent with that of our correspondence factor analysis (Figure 12C, D), that is an inverse relationship between loess and archaeological occurrences.

Regional studies often highlight the proximity of occupations to large rivers. To test this potential association, we calculated the distance to the nearest watercourse (wider than 15 m) for each occurrence. The result (Figure 14) shows that (1) the majority of occurrences (57%) are located less than 2 km from a river; (2) their number decreases rapidly with distance from the river; (3) the distribution of 200 random samples of the same number of points differs significantly from this pattern. For the randomization experiment, on average, only 19% are less than 2 km from a river, and the proportion of remaining sites decreases only slightly with distance. Proximity to large waterways is, therefore, a parameter that influences site distribution.

The analysis of occurrence location in relation to raw material source areas indicates that approximately 66% of them are located on flint (or related rocks) formations, while random points fall on average on such formations only approximately 41% of the time (Figure 15). When sites are surrounded with a 10 km buffer, it is apparent that almost all occurrences (97%) are located on or adjacent to flint formations. This proportion rises to over 99% when a 30 km radius buffer zone is employed. A radius of 30 km corresponds to the distance that covered in one day when walking over relatively flat terrain. Examination of the map (Figure 16) shows that the areas where occurrences are least numerous, or absent, correspond to areas where the Hercynian plutonic and metamorphic basement outcrops (Brittany, Massif Central, Vosges Mountains), but also to a portion of the Pyrenean foothills consisting of Cenozoic sandy and clayey detrital formations (molasse). This map also shows that the association between occurrences and flint-containing formations, mainly carbonate rocks, is not directly linked to the presence of karstic cavities. While areas in Southwestern France where formations that contain flint outcrops are associated with high frequencies of occurrences situated in caves and rockshelters, this is not the case in the Paris Basin and northern France, where open-air occurrences predominate. Based on the fact that this relationship holds true regardless of site type (cave, rockshelter, open-air), indicates that a clear relationship between raw material sources and sites exists.

4. Discussion

4.1. Interpretation of the Various Factors Explaining the Distribution of Occurrences

The relationship between current land cover and occurrences clearly indicates that certain environments are more favourable than others for the discovery of Paleolithic sites, particularly peri-urban areas, vineyards, forests, and, to a lesser extent, quarries. This finding largely concurs with that already made by [51] in Europe. Peri-urban areas and quarries are regularly subject to archaeological surveys, particularly in the context of preventive archaeology, and are characterized by abundant discoveries. The overrepresentation of finds in vineyards reflects the fact that the ground between the rows remains bare throughout the year, which is favourable for archaeological survey (at least until farming practices changed in recent decades to favour the development of grass cover). Conversely, occurrences are underrepresented in cultivated areas, where the length of the period without vegetation is reduced, and even more so in pastures. The abundance of discoveries in forested areas is more counterintuitive. The reason is that caves and rockshelters are mostly located in areas that are not suitable for agriculture (rock walls, scree) and, therefore, dominated by natural forest vegetation. While this variability in discoveries based on land use does represent a bias in the archaeological record [51], its overall impact on the representativeness of the database of occurrences is likely to be insignificant, mainly because the areas of under-representation or, conversely, over-representation of sites are distributed across the study area. Preventive archaeology, which has developed extensively in France since the 1990s, has sampled across the entire study area and has, thus, mitigated biases generated previously.

The analysis also shows that other variables contribute significantly to explaining the distribution of occurrences during the period 34–21 cal ka BP, including indirect climate proxies, as well as environmental parameters not directly related to climate. All variables that strongly influence local climate, such as altitude, latitude, and exposure, reveal a consistent pattern characterized by a marked decrease in occupation with latitude and altitude, with a potential continuous habitability limit of around 600 m, and orientations favourable to maximum sunlight. This finding is consistent with observations made previously [1,3,4,13,51].

Overall, low temperatures limit occupation. The ratio of occurrences per thousand years shows a sharp decrease during the GS-3 stadial, one of the coldest phases of the Last Glacial. This phase coincides with the transition between the Gravettian and the Early Solutrean. It is worth noting that it coincides with the genetic bottleneck identified in Western Europe by [18,19], thus suggesting a sharp decline in population. In turn, this suggests that, among all the factors likely to influence the number of occurrences per period (population density, seasonal activity cycles, mobility patterns, group size), population density is likely the most influential.

The southern limit of continuous permafrost, indicating a mean annual air temperature of ca. -3 to -5 °C, effectively delimits the northeastward extension of the Middle Gravettian, Final Gravettian and subsequent technocomplexes, whereas the Early and Late Gravettian are not affected by this limit and have a wider distribution. We propose two possibilities for this contrasting pattern. The first is that Early and Late Gravettian occupations in the northeast portion of our study area are largely contemporary with interstadials (relatively milder conditions). Stratigraphic distributions of Paleolithic levels in Northern France, Belgium and Central Europe support this (e.g., [68,69,70]). Furthermore, numerous researchers propose a hiatus in the occupation of Northern and Central Europe during portions of the LGM [11,12,71,72]. The second possibility is that the Solutrean and Badegoulian correspond to Western European technocomplexes present in a relatively less rigourous oceanic climate, while another technocomplex, the Epigravettian, occupied Central Europe under cold, continental conditions. This probably reflects a different culture-environment relationship between the North-Central European Epigravettian and the Solutrean-Badegoulian.

Among the environmental parameters unrelated to climate, proximity to sources of lithic raw materials appears to be a determining factor. This finding is surprising at first glance, given that the transport of raw material blocks, cores, or tool blanks over distances ranging from one to several hundred kilometres has been documented at numerous sites (e.g., [40,73,74,75,76]). The displacement of lithic raw materials, as well as ornaments [77], connected distant regions and involved movement across areas that have yielded no evidence of occupation. Aubry et al. [76] make a similar observation in Portugal and Western France. First, site visibility depends on the quantity and representativeness of the material culture remains present, which is influenced by both mobility and site function. In a given region, sites resulting from brief occupations associated with specific activities coexisted with sites where occupations were longer and activities more diverse. The latter are probably overrepresented in the database because they are more easily identifiable during archaeological surveys. The lack of such sites across large areas, however, raises questions. To explain this pattern, one hypothesis is that settlements located far from sources of raw materials are less visible in the archaeological record. Raw materials were used at sites close to sources without concern for economy. Large quantities of waste material were discarded during occupation, making them easy to detect during surveys. Conversely, raw materials transported over long distances were often curated, i.e., resharpened and reused rather than being discarded [78,79]. Sites dominated more by curation behaviours are less visible and may contain few diagnostic elements, making them difficult to attribute to a particular archaeological culture. A number of occurrences attributed to undetermined late Paleolithic technocomplex exist in the study area and have not been included in the database. According to this assumption, the distance to raw material sources would not constitute a limitation to settlement; rather, the absence of occupations more likely reflects a bias in the archaeological record related to their visibility.

The proximity of a large river is also a predictor of site distribution. There may be multiple reasons for this, including (1) access to drinking water, although a small river or springs would suffice. Therefore, this explanation is unlikely; (2) greater plant biomass in valleys compared to plateaus. Available few pollen analyses and vegetation simulations suggest the persistence of patches of shrub vegetation in valleys during the LGM (e.g., [33,80]). This would have increased available resources; (3) at least seasonally abundant animal biomass, with large valleys serving as migration corridors for large herbivores; and finally, (4) the presence of karstic cavities (rockshelters, caves) conducive to the establishment of sheltered camps along valleys in limestone terrains, particularly in Southern France. These different factors are not mutually exclusive, but together likely created conditions favourable to the occupation of large valleys.

The analysis also highlights a negative association between aeolian deposits and sites. The reasons for this may be different for areas dominated by sand and loess. For sand deposits, which are mainly located in Southwest France (“Sable des Landes”) where they cover an area of approximately 13,000 km², the hypothesis of a desert environment poor in plant and animal resources and, therefore, little used by hunter-gatherers [35] remains valid. Outcrops of cretaceous flint-bearing limestone exist in a small valley (the Ciron valley) that crosses this sand region (Villagrains flint: [81,82]). These flint sources were exploited during the Middle Paleolithic and the Late Magdalenian periods [83,84], but no LGM occupations are known in the vicinity. The aeolian sands of Sologne situated south of Orléans (Paris Basin), albeit a much smaller area (approximately 2,600 km²), provide a similar picture, with no occurrences recorded in this area.

On the other hand, with respect to loess deposits, which are present mainly in Northern France (Northern European Loess Belt, NELB) and deposited in a steppe-tundra environment capable of supporting large herbivore populations, another explanation seems plausible. Preventive archaeology data show that few Paleolithic sites have been discovered in this broad area, mainly from the Grand-Ouest (GO) and Hauts-de-France (HdF) administrative regions. Similarly, there are few sites in the Grand Est (GE) region, where loess cover is virtually non-existent. This suggests that there is no causal relationship between the presence of loess and the absence of sites. The hypothesis that a thick layer of loess limits the detection of sites during archaeological surveys alone does not explain the observed pattern. The scarcity of sites in the NELB is likely due to other factors, particularly climate, as indicated above. In contrast, it is worth noting that there are numerous occurrences in the loess band situated around the “Sable des Landes” in Southwest France (i.e., at lower latitudes in relatively close proximity to the Atlantic Ocean and thus milder climatic conditions), and this for the entire late Paleolithic [35].

4.2. Implications

The results of this study have several implications. We can reject, with a reasonable degree of confidence, the hypothesis of a significant bias in the distribution of sites linked to regional differences in archaeological survey and discovery. In particular, discoveries made over the last 30 years in the context of preventive archaeology have only marginally altered the distributions of sites described in early studies [1,13]. Southern France consistently appears to be the region with the highest concentration of Paleolithic sites.

We can also reject the hypothesis of an important bias due to issues of site visibility or discovery, particularly in areas with a thick loess cover. Although the analysis reveals a negative association between loess and occurrences, we can establish no causal link between the two, as one can observe areas of low occupation density in Northern France both within and outside the NELB. Conversely, the loess areas of the Aquitaine region in Southwestern France have high numbers of archaeological occurrences like other regions of loess in Central Europe.

Nevertheless, one cannot state that the distribution of known sites is free of any bias and thus provides a complete picture of past settlement. The highly significant relationship between the distribution of areas yielding raw materials suitable for knapping and the distribution of occurrences strongly suggests that the latter’s visibility during surveys plays an influential role. Thus, it is likely that large “void” areas on the map do not reflect a total absence of settlement, but are likely a testament of our ability to recognize or characterize occurrences with few lithic artifacts. This is because the frequency of flint artefacts in an assemblage and distance from lithic source areas are inversely related [85]. This hypothesis contradicts the idea of fragmentation of settlement areas during the coldest periods of the Last Glacial.

Analysing Paleolithic settlement in France using climate as the sole explanatory factor seems reductive, even though climate had a significant impact. Other locally and regionally expressed environmental factors, such as sandy deserts and large valleys, had significant influences that impacted broader geographic patterns to varying degrees. Climate simulations potentially offer an incomplete means of evaluating human settlement during the Pleistocene and, therefore, should be complemented with other data so that a broader range of relevant factors are considered when evaluating past settlement dynamics.

Since present-day France represents only a small portion of the European continent, it is pertinent to ask whether the importance of the parameters we identified remain relevant at this broader geographic scale. Examinations of available maps provide some means of assessing this question. Figure 17 illustrates the distribution of European Gravettian-age sites (those from this study along with those compiled by [17]) superimposed on a map of aeolian deposits [61]. The map highlights the role as a barrier to settlement played by the northern European Sand Belt (ESB). With the exception of one site in Poland, no occupations are recorded in the ESB, whereas sites are known to exist in the United Kingdom at higher latitudes. Boemke et al. [51] made a similar observation. As in Southwest France [35], they also showed that Late Magdalenian and Epipaleolithic occurrences, which are contemporary with the Lateglacial recolonization of these sands by vegetation, are, conversely, well represented.

Since limestone is the main flint-bearing formation, the geological map of Europe at a scale of 1:5,000,000 (IGME 5000) [86] provides a general, albeit approximate, overview of the distribution of raw material areas. Comparison with the distribution of LGM occupations highlights their preferential location on or near limestone outcrops (Figure 18). This pattern is particularly clear in the Iberian and Italian Peninsulas, as well as in Eastern Europe, as already noted by [39]. Boemke et al. [51] also identified an overrepresentation of sites on Jurassic and Cretaceous limestone and favoured the hypothesis that the choice of site location was tied to the presence of karstic cavities. With respect to the French Late Paleolithic record, our results allow us to reject this hypothesis (see above).

5. Conclusion

The analysis of occurrences belonging to the Gravettian, Solutrean, and Badegoulian technocomplexes in France, via GIS data and statistical methods, provide a means of exploring the possible role of numerous environmental factors in the observed distribution of Paleolithic sites between 34 and 21 cal ka BP. Our results indicate that occurrence distributions deviate significantly from random expectations. Climate, inferred here from indirect proxies (latitude, altitude, exposure), influenced observed site distributions and, particularly, explains the north-south gradient in site frequency, with southern France consistently being the region with the highest concentration of Paleolithic sites. The lowest frequency of sites corresponds to GS-3 (ca. 27–24 cal ka BP), one of the coldest events of the Last Glacial. The analysis also highlights positive or negative associations between raw material sources, major rivers, aeolian deposits and the site distribution. These correlations can be explained either by a bias in the archaeological record, due to the difficulty of identifying and characterizing occurrences poor in lithic material located far from raw material sources, or by these factors being favourable or, conversely, unfavourable to the settlement of particular areas by Paleolithic groups. Assessing the relative role of individual factors is not a trivial undertaking, and future work will necessitate the development of appropriate statistical approaches. Ultimately, this study suggests that approaches based solely on paleoclimatic simulations do not provide a complete understanding of Late Paleolithic settlement patterns. We propose that future studies should take into account additional data, particularly geological and geomorphological in nature, when examining settlement dynamics.