Changing Movements in a Changing World: Modelling Early Pleistocene and Early Middle Pleistocene Climatic and Ecological Environments and Influences on Hominin Dispersal in Eurasia

2.1.1 Climatic reconstructions

Polvani et al. (2017) and Rocha (2023) presented the emulator Earth System Models (ESM) to simulate the Earth’s climate through biological, chemical and physical processes. Due to the heavy data load of such complex simulations, simpler climatic reconstructions have been produced by the Community Earth System Models (CESM) to facilitate applicability to analyses (Polvani et al. 2017; Rocha 2023). The latest version of CESM2 generates outputs from the atmosphere (i.e., precipitation, cloud forcing, land, ocean, wave models, river run-off, land and sea ice components), which can be combined depending on the research questions and software. The time step follows orbital sequences and can be down-scaled to 30-minute steps. The spatial resolution is approximately 1° × 1°, with varying resolutions across components. For pre-industrial periods, the reconstruction has been calibrated with two control simulations, respectively going back in time 500 and 1,200 years (Danabasoglu et al. 2020), and with PaleoCalibr, e.g., in the work of Zhu et al. (2022). The reconstructions have not been scientifically validated, but they include control simulation to ensure the model’s reproducibility (Yeager et al. 2006). The estimates show increasing uncertainties the further the model goes back in time (Zhu et al. 2022). However, as demonstrated by Lemmen (2007) and Zhu et al. (2017, 2022), they can be useful in studies of past environmental conditions (Lemmen 2007; Zhu et al. 2017, 2022).

Another model extending back to the Last Glacial Maximum (LGM) is the Hadley Centre coupled model (HadCM3), developed by Gordon et al. (2000). It emulates a 19-layer atmospheric component, a 20-layer ocean component with sea bottom topography, and sea ice estimations of ice concentrations and movements. It has monthly time steps and a spatial resolution of 3.75° × 2.5° for the atmospheric model, and 1.25° × 1.25° for the ocean model. The reconstruction has been validated by comparisons with observed ocean heat transport dynamics and temperature classes from the Atlantic and the Pacific. These generally concur, yet with a slightly overturning element in the HadCM3 emulator, causing overestimations of heat transported northwards by western boundary wind-driven components and southwards by interior components (Gordon et al. 2000). However, with the high spatial resolution of the ocean, the authors do not expect realistic boundary currents (Gordon et al. 2000; Valdes et al. 2017). The fact that the emulator was defined on the basis of post-industrial climatic data presents a challenge, and complicates its use in pre-industrial periods. However, Gordon et al. (2000) showed that the difference between reconstructions and real-world climate was 0.3 °C, which only represents minor uncertainty in long-term global estimates (Gordon et al. 2000).

Paleo-Clim by Brown et al. (2018) focuses on the Pleistocene climate, is complemented by atmospheric-ocean general circulation models (AOGCMs) from HadCM3, and simulates annual mean temperature, temperature of the coldest and warmest months, annual precipitation, the extremes of the wettest and driest quarters and the effects of topography. This model differs from other reconstructions by applying a snapshot approach, generating chronological windows at MIS19 (787 ka), during the mid-Pliocene warm period (3.264-3.025 Ma), and MIS M2 (3.3 Ma), with a pixel-to-pixel resolution of 0.04° × 0.04°. These reconstructions have been validated by comparing raw differences and anomaly ratios of precipitation data with corresponding values from HadCM3 reconstructions. Pearson correlations between the reconstructions, global visualisations and raw climate model values were generally high, although using raw differences may result in underestimated outputs. However, the authors argue that the snapshots offer detailed insights into certain moments in time, functioning as precise cases for the exploration and validation of climatic drivers of ecological change (Brown et al. 2018).

The Paleo-PGEM climate emulator by Holden et al. (2019) generated the repository Paleo-PGEM-series (Barreto et al. 2023), presenting dynamic spatiotemporal climate estimates. The emulator generates detailed resolutions for detecting climate variability forced by orbital changes with PLASIM-GENIE, an AOGCM containing a 10-layered vertical climatology (Holden et al. 2019), BIOGEM, an ocean model (Ridgwell et al. 2007) and an efficient numerical terrestrial scheme (ENTS) for terrestrial carbon cycle modules (Holden et al. 2018). The simulation runs in orbit cycles, with varying concentrations of atmospheric CO₂ and ice sheet coverage, monthly temperature and precipitation estimates, and 17 bioclimatic variables. The model is validated by comparison with observation-based global temperature reconstructions for the LGM, the Mid-Holocene and the Last Interglacial, and proxy data from different periods and regions, showing a minor deviation in the spatial pattern of temperature anomalies around 127 ka (Barreto et al. 2023). The 1,000-year resolution yielded uncertainties in the feedback between atmospheric circulation and topography due to incomplete computational methods or a lack of knowledge on climatic variables. For this reason, the authors emphasise that the model should be considered as a realistic reconstruction for distinguishing spatial variabilities, and not as a complete high-resolution reconstruction (Holden et al. 2019).

Gamisch (2019) presents Oscillayers data to achieve high-resolution spatial and temporal climate estimates. Using 19 bioclimatic parameters, such as mean and extreme values of temperature and precipitation, Oscillayers provides proxy indicators of terrestrial climates and sea level changes. The spatial resolution of the model is 0.4° × 0.4°, with 10,000-year time steps, extending back to 5.4 Ma. After comparisons with general circulation models (GCMs), the reconstructions correlate with other Holocene, LGM and Pliocene climatic estimates. The author argues that reconstructions are useful for testing eco-evolutionary hypotheses when accompanied by ecological niche models. The main advantage of Oscillayers over dynamic climatic reconstructions is its detailed spatial resolution beyond 800 ka (Gamisch 2019). However, Barreto et al. (2023) express concerns about limits and strong assumptions on data, especially with regard to extrapolating changes linked to the extensions and contractions of the LGM ice sheet during the Pleistocene (Barreto et al. 2023).

PMIP is the palaeoclimatic segment of the Coupled Model Intercomparison Project (CMIP5), which generates high-resolution reconstructions of past climate. The current version of the emulator (PMIP4) encompasses atmospheric processes, ocean circulation, sea surface temperatures, sea and land ice, soil moisture and hydrology, carbon cycles and vegetation, and generates outputs of air temperature, precipitation, soil moisture, snow cover and depth, sea surface temperature, ocean heat and currents, the extent of sea ice and vegetation fraction. The spatial resolution of the reconstruction can be as low as 0.5° × 0.5 °, but for most reconstructions, resolutions are between 1° × 1° or 2.5° × 2.5°, depending on the resolution of the input data. As with Oscillayers, the PMIP4 emulator provides time-specific snapshots, which for now have been produced for the LGM (21 kyr), the mid-Holocene (6 kyr), and a transient simulation from 850–1850 CE (Kageyama et al. 2018; Taylor, Stouffer and Meehl 2012). The reconstructions have been validated against proxy reconstructions from ice cores, marine sediments, pollen and speleothems by direct comparison and forward modelling, and through benchmarking tools and standardised control simulations (Kageyama et al. 2018). However, the model struggles to generate more complex climate systems, such as the relationship between ice sheet volume, ocean heat transport, and cloud properties in the Arctic (Kageyama et al. 2021). This results in challenges for estimating accurate cooling in the northern and southern Arctic areas (Ohgaito et al. 2021).

All the reconstructions and emulators presented above (Table 1) combine global-scale temporal and spatial applicability, the ability to estimate past climates and potential for combination, calibration, and validation against climate proxies (Chevalier et al. 2020; Heiri et al. 2014; Krapp et al. 2021). Proxy records are routinely calibrated with instrumental records, thereby increasing their reliability and extending their range forward and backwards in time (Frank et al. 2010). However, extending back as far as the Early Pleistocene can be challenging for some proxies. To ensure reliable global reconstructions of the Early Pleistocene, isotopic data from ice core samples are often used for validation (e.g., Berends et al. 2021; Grant et al. 2025; Kageyama et al. 2018), while more regional proxies include pollen and faunal records (Bacon et al. 2023; Berto et al. 2024; Birks 2019; Herzschuh et al. 2023; van der Made, Morales and Montoya 2006). Having established the model’s reliability, the specific purpose of the reconstruction must be considered. Modifications to temporal limits, resolution, or the addition of other climatic variables highlight how the choice of emulators must be linked to the research question. Moreover, adding environmental elements to climatic reconstructions – such as net primary productivity variables, as shown by the work of Krapp et al. (2021), or vegetation feedback represented by leaf area index (LAI), as shown by O’ishi et al. (2021) – creates a broader impression of past environmental conditions. For this reason, the next section introduces emulators for environmental reconstructions that capture ecological variations. It is fitting to use such emulators alongside climatic reconstructions.

2.2.1 Environmental emulators

LOVECLIM, by Goosse et al. (2010), is a coupled climate and environmental emulator, useful for both climatic and basic environmental estimates. It generates glacial temperatures, climatic anomalies (such as Heinrich and millennial-scale Dansgaard-Oeschger events), hydroclimatic variabilities, and includes a general ocean circulation model (CLIO model), a thermodynamic sea ice model (Goosse and Fichefet 1999), and a terrestrial vegetation model following annual mean values of climatic variables (VECODE by Brovkin et al. 1997). The emulator has a grid-size of 3° × 3°, but can be downscaled to a 1° × 1° resolution, with yearly timesteps. The reconstructions are validated against 63 proxy data reconstructions of the period 140–20 ka. The comparison showed a general correlation with other reconstructions, with minor underestimates of the reconstruction for the Last Interglacial Ocean warming. Otherwise, the estimated surface temperatures matched the ice-core data and the estimated anomalies in rainfall ratios. However, the reconstructions present biases linked to model formulation, and cannot be expected to simulate the same details as a GCM. Still, Goosse et al. (2010) argue for the emulator’s utility for the environmental reconstruction of past periods (Goosse et al. 2010).

A further development of this model is iLOVECLIM, which retains the same core elements (atmosphere, ocean, vegetation), in addition to the water oxygen isotope cycle, an oceanic carbon model, an iceberg trajectory module, and the interactive ice sheet model GRISLI. The spatial resolution of the emulator is the same as for LOVECLIM and has, for now, been used for reconstructions of the LGM. Estimates were evaluated against an observation-derived data set of mean annual surface temperatures and long-term mean climatologies of annual precipitation. In general, latitudinal cooling is reflected in the estimates, while estimations of mean temperatures for North America, Greenland and Western Europe are too warm (Quiquet et al. 2018; Roche et al. 2014).

The PRISM4 emulator by Dowsett and colleagues (2016) includes various variables to generate climate and vegetation reconstructions, including palaeogeography, sea levels, ocean temperatures, land and sea ice cover estimations, lakes and soil. The resulting output represents terrestrial vegetation using biome classes. It is a conceptual model of environmental conditions and landmass changes during the mid-Piacenzian (3.6–2.6 Ma), with a spatial resolution of 0.25° × 0.25°. The time-steps follow the chronology of marine isotopic stages by Lisiecki and Raymo (2005), with 1,000-year steps until 0.6 Ma, 2,000-year steps between 0.6 and 1.5 Ma, and 2,500-year steps between 1.5 and 3 Ma. Since the model is intended as a conceptual model for boundary conditions for other environmental reconstructions, it has not been scientifically validated. However, PRISM4 has proven useful when used to verify other models, with orbital forcing correlating with present-day forcing and isotope peaks. The authors argue that the multiproxy chronology is useful. However, it is still important to consider calibration uncertainties and the orbital-scale chronological resolution, which are more suited to regional than global-scale analyses (Dowsett et al. 2016).

The ORCHIDEE model, presented by Guimberteau et al. (2018) and Krinner et al. (2005), is a land-surface emulator which represents spatial variation in environments through different outputs: 13 plant functional type (PFT) classes, gross primary productivity (GPP), NPP and biomass values. The outputs are calculated from the internal dynamics of soil properties and feedback on water and atmospheric CO₂ concentrations. The reconstruction has a spatial resolution of 4° × 2.5° for global simulations, with 1-year time steps, or determined by the resolution of the climatic input data (Krinner et al. 2005). Simulation outputs were validated against real temperature gradients, basin-scale averages for the dominating waterbodies, other topsoil moisture models, and forest site data. These indicate that the model has a slight tendency to overestimate carbon density in some regions, including North-Western Europe. Still, the outputs of vegetational units are argued to be reliable representations of arctic vegetation, due to their separate calibrated parameters (Guimberteau et al. 2018).

The environmental emulator BIOME4 by Harrison and Prentice (2003) uses atmospheric CO₂, temperatures, mean monthly precipitation, soil physical properties, and solar radiation at different latitudes to infer changes in PFTs grouped into 28 biome categories. The model has a spatial resolution of 0.5° × 0.5°, and 1,000-year time steps going back to the LGM. For validation purposes, the reconstruction is accompanied by a site-based vegetation map, derived from pollen and plant macrofossil data and divided into PFT groups. Overall, the reconstruction matches the observed forest cover with minor overestimates in higher latitudes and underestimates in the tropics (Harrison and Prentice 2003). Even though the time range of the emulator does not correspond to the chronology of this paper, it is included because of its use in Early Pleistocene studies, e.g., in Romanowska et al (2017). BIOME4 was the only available emulator for dynamic environmental reconstructions in 2017, despite the fact that the emulator was constructed on the basis of LGM data (Romanowska et al. 2017). However, other climatic and environmental emulators have been published since then, such as the Paleo-PGEM Series, enabling the use of suitable chronological reconstructions.

The Carbon Assimilation in the Biosphere (CARAIB) model is a process-based dynamic vegetation emulator used to study global vegetation and its role in the past and present carbon cycle. The emulator contains modules for hydrological budget, heterotrophic respiration, soil carbon dynamics, photosynthesis and stomatal regulation, carbon allocation, plant growth, plant competition, biogeography and fire disturbance. Combined, these modules produce outputs of soil hydrology, biomass, GPP, NPP, LAI and biome distribution by bioclimatic affinity groups (BAGs). The model can produce output down to a spatial resolution of 0.5° × 0.5°. However, this is determined by the resolution of climatic data, which also determines temporal resolution (Dury et al. 2011). Still, given the spatiotemporal scope of this paper, the emulator is problematic, namely in terms of the resolution requirements of input data. The emulator requires climatic input data (air temperature, precipitation, relative humidity, relative sunshine hours, wind speed) at a daily resolution (François et al. 1998), which can be challenging to obtain for the Pleistocene. However, reconstructions have been generated back in time, for example by François et al. (1998, 2006) for the LGM and the late Miocene. These reconstructions were generated with climatic data from present GCMs, and their outputs have been validated against other simulated reconstructions and real CO₂ measurements (François et al. 1998, 2006).

Environmental emulators and reconstructions, with different resolutions, variables, and outputs, clearly provide various levels of insights into past conditions (Table 2). Reconstructions vary in terms of temporal resolution, time range and output data, which are important aspects to consider when choosing a suitable reconstruction.

Although Timmermann and Friedrich (2016) argue that climate is the dominant driver of hominin dispersal, the influence of vegetation must not be overlooked. A comparison of the study by Timmermann and Friedrich (2016) with that of Zeller et al. (2023), which adds the biome factor, shows how additional information influences the dispersal options of hominin populations. The latest models of environmental impacts on hominin dispersals out of Africa have demonstrated that high-resolution environmental data could enhance our understanding of the influence of vegetation on dispersal behaviour (Hughes et al. 2007; Romanowska et al. 2017).

Yet, linking climatic and environmental reconstructions does not satisfactorily reflect their influence on hominin behaviour, and further theoretical and methodological tools are necessary to connect them. The next section presents different theories and simulation techniques for analysing the hominin-environment relationship in detail.

2.3.1. Static models

In addition to palaeoenvironmental reconstructions, another important element is missing when investigating early hominin dispersals into Eurasia during the Early Pleistocene: the hominins themselves. Different hominin species are known to have dispersed outside of Africa. However, various obstacles arise when we attempt to address hominin diversity. Firstly, there are limited hominin fossils in the Early and early Middle Pleistocene archaeological record, which complicates precise species divisions (Arzarello et al. 2007; Gabunia et al. 2000; Gibert et al. 1998; Huguet et al. 2025; Larick et al. 2001; Lorenzo et al. 2015; Palmqvist 1997; Swisher et al. 1994). Secondly, there are no clear definitions of physiological and cultural differences between species and thus no insights into how such differences influenced behaviour (Antón and Josh Snodgrass 2012; Duda and Zrzavý 2013). Therefore, the approach taken in this paper assumes that the hominin population is homogeneous.

Niche modelling, or ecological niche modelling (ENM), is a computational tool for analysing species’ distribution potential and actual geographical distributions. These algorithms derive from the discussion centring on Grinnellian niches and the principles of behavioural ecology, and use climatic or environmental data and presence/absence data to predict the geographical distribution of a species or population (Araújo and Guisan 2006; Austin and Smith 1989; Fitzhugh et al. 2019; Pearson and Dawson 2003; Soberón 2007).

Examples of these concepts are species distribution models (SDMs), which predict the presence or abundance of a species depending on climatic and/or geophysical characteristics, substrates, and nutrients. Simulations have been used to identify ecological variables that predict species distribution probability and spatial preferences, trained on environmental and species data in either absence, presence/absence or abundance data sets (Elith and Leathwick 2009; Guisan and Thuiller 2005; Kearny and Porter 2009; Simoes et al. 2020; Soberón 2007). The MaxEnt simulation tool is used to generate SDMs from presence-only data and environmental parameters relevant to habitat suitability. Such studies show that presence-only data provide a foundation for distribution analysis, which can be compared to other case studies (Hijmans and Graham 2006; Valavi et al. 2022). Furthermore, SDMs can be improved owing to the growing number of occurrence records and high-resolution spatial environmental data sets (Phillips, Anderson and Schapire 2006). Another tool for simulating SDMs is tidysdm, an R package containing correlative SDMs which benefit from the existing functions of tidymodels in an easy-to-adapt modelling framework. With the tidysdm package, it is possible to use climatic data from several periods, either in conjunction with presence data or using the pastclim tool, and it is thus ideal for addressing archaeological questions (Leonardi et al. 2024). Franklin et al. (2015) emphasise how SDMs can accompany niche models to discover and validate environmental parameters when location data are missing (Franklin et al. 2015). Despite the complex interplay between scale, modelling method and data quality, SDMs have been used for large-scale species distribution mapping (Elith et al. 2006, 2011), and in archaeological research on hominin dispersal and habitat preferences (Franklin et al. 2015; Padilla-Iglesias et al. 2022; Yaworsky, Nielsen and Nielsen 2024). Close relatives to SDMs are habitat suitability models (HSMs), which identify influential environmental parameters limiting a species’ spatial range using environmental and presence data (Guisan and Thuiller 2005; Thuiller 2024). However, neither HSMs nor SDMs allow for the observation of the temporal evolution of the system. Without a temporal perspective, it is impossible to observe how population dynamics (i.e., population structure), dispersal, abundance and local extinction evolve depending on their adaptation or preference for a set of climatic or environmental variables. To overcome this limitation, authors such as Thuiller and Münkemüller (2010) suggest linking HSMs with individual-based modelling methods, such as ABMs (Thuiller and Münkemüller 2010).

2.3.2. Dynamic models

Since population dispersal is a dynamic process evolving over time, static insights from niche models and SDMs are not suitable for studying the causal mechanisms between climatic and environmental changes and the process of hominin dispersal. Instead, this question demands dynamic models (Romanowska, Wren and Crabtree 2021). One of these simulation methods is cellular automata (CA), which uses a grid-based world where a numerical state value is attributed to each cell, often ranging within a scale or a Boolean value. Spatial scale is determined by the number of cells, cell-size, and interaction-reachability of the agent within the cell-neighbourhood. In dispersal models, the numerical value of the cell represents the probability of occupation, and the Boolean value indicates the presence or absence of a population. This type of simulation is appropriate for large-scale studies, where the population is viewed as a single coherent entity, and less suitable for investigating inter-population dynamics (Bithell and Macmillan 2007; Romanowska 2015), as shown through the use of CA models in studies of global-scale dispersals of hominin populations.

The first Out of Africa CA by Mithen and Reed (2002), named Stepping out, featured a landscape arranged in hexagons, a climate following glacial-interglacial cycles, two sea-level heights, and a single hominin population with various dietary and environmental preferences. At each time step, the population could stay and adapt to the current cell and/or colonise a neighbouring one, starting their dispersal in East Africa. Simulations showed how adaptive specialisation and the colonisation rate influenced dispersals and thereby how environmental factors could evaluate dispersal routes, conditions, and chronology (Mithen and Reed 2002).

Nikitas and Nikita (2005) expanded Stepping out with three landscape types and three population dispersal modes, the first with a stable total population size, the second with population growth based on movement, extinction and colonisation rates, and the third as a combination of both scenarios. In this way, they tested the impact of natural obstacles and dispersal dynamics on early dispersal timing and routes into Eurasia. The results demonstrated that dispersal dynamics had a greater impact on dispersal timing than environmental barriers, and that a stable population size was more probable than a high population growth rate (Nikitas and Nikita 2005). Hughes et al. (2007) added environmental elements from BIOME4 to Stepping Out and modified the simulation into a null-hypothesis model by testing the influence of the environment on the dispersal pattern of a simplified non-developing hominin population. They observed an increase in the chronological predictability of occupation, thus arguing for a greater focus on ecological niches in hominin dispersal (Hughes et al. 2007).

This perspective was included in the latest modelling attempt by Romanowska et al. (2017), who investigated whether dispersal and population dynamics could explain the distribution of archaeological assemblages in relation to the Movius line. Cell biome values were generated using an equation-based modelling technique, with an environmental setup including temperature, sea-level curve, and eight biome categories connected to these variables. The environment was built with the structures of a CA, with grids containing biome-specific values for carrying capacity, population dispersal probability and population growth rates. The simulations emphasised that local environmental conditions and levels of hominin adaptations were more relevant to hominin population density than an increased distance from the population centre (Romanowska et al. 2017).

Agent-based modelling (ABM) offers inter-population insights by adding a third component; agents, along with a spatially structured environment and rules of behaviour defining the interactions between the two. Agents and the environment possess different attributes, making it possible to observe changing properties as the result of their interplay, thereby providing a bottom-up insight into large-scale patterns (Romanowska, Wren and Crabtree 2021; Wilensky and Rand 2015). Used in a number of disciplines, the impact of this method in archaeology has steadily increased over the last couple of decades (Cegielski and Rogers 2016; Chliaoutakis and Chalkiadakis 2016; Davies et al. 2019; Graham 2020; Lake 2015; Reeves et al. 2023; Romanowska et al. 2019). ABM’s bottom-up approach inspired the study of hominin dispersals, since it allows for the observation and testing of different scenarios, a rare opportunity in archaeological research (Hölzchen et al. 2016; Hughes et al. 2007; R. Vahdati et al. 2019; Romanowska et al. 2017). Agent-based models have been used to test different drivers of dispersal, such as environmental conditions (d’Alpoim Guedes et al. 2016; Hughes et al. 2007; Perry et al. 2016), demographic dynamics (Coto-Sarmiento et al. 2023; Cucart-Mora, Lozano and Fernández-López de Pablo 2018; R. Vahdati et al. 2022; Romanowska et al. 2017), and cultural and cognitive abilities (Barton et al. 2011; Reeves et al. 2023; Wren et al. 2014).

Bioco et al. (2022) presented SDSim, a simulation approach connecting the principles of SDM and ABM. Based on the use of important variables for species ecology, the SDSim generates a raster habitat suitability map and combines it with agents’ birth, reproductive, death, and distribution rates (Bioco et al. 2022).

The growing number of simulation studies reveals increasing confidence in the methodology and its potential for investigating climatic and environmental effects on dispersal dynamics. However, in the same way as for palaeoenvironmental reconstructions, simulated dispersal patterns require validation. Given the limited number of sites, it is vital to take stock of the risk of overfitting the dispersal simulation to the available data. This is often a challenge when working with archaeological data, which is why various approaches have been developed to evaluate the proximity between simulation results and archaeological data, e.g., matching the model to stylised facts (i.e., generalised patterns, e.g., Romanowska 2014; Romanowska et al. 2017), scenario testing (e.g., Griffin and Stanish 2007; Wren et al. 2014), and parameter sweeping (e.g., Hughes et al. 2007; Nikitas and Nikita 2005). These examples show how to heighten awareness of biases in the archaeological record, thereby strengthening the use of data. In addition to comparisons to archaeological data, simulation results can be compared to the output of previous simulations, such as those presented above (Hughes et al. 2007; Romanowska et al. 2017). Comparisons with simulations covering the same chronology and scale would reveal whether finer resolution data increase or decrease the precision of the simulation’s predictive accuracy. Another intercomparison could be made with equation-based models, comparing results across different modelling approaches (Callegari et al. 2013).

2.3.3 The optimal method for modelling hominin dispersals out of Africa

Given the various options for climatic and environmental reconstructions, emulators, and simulation techniques, the next question is which combination is the most suitable for the research question examined here: how and which climatic and environmental factors influenced spatiotemporal hominin dispersal patterns during the Early and early Middle Pleistocene? Here, a conceptual workflow proposal presents the general steps of connecting environmental data with a dispersal simulation (Figure 1). Detailed tutorials of building a simulation can be found in Romanowska (2015), and instructions for specific software for building reconstructions are often provided on accompanying webpages or documents, e.g., Paleo-PGEM Series (PALEO-PGEM-Series: a spatial time series of the global climate over the last 5 million years (Plio-Pleistocene)) or BIOME4 (GitHub – jedokaplan/BIOME4: The BIOME4 equilibrium global vegetation model).

Figure 1

As demonstrated in Figure 1, the process starts with the best suited climatic or environmental reconstruction to the research question, generated from the emulator and data set. In Tables 3 and 4, the options are listed along with considerations on the suitability and combinations of data and emulators for answering this question.

Factors limiting the use of the presented data sets or emulators are low spatial or temporal resolution or the model’s temporal limits. In terms of our research question, a temporal range back to the Early Pleistocene is indispensable. This criterion rules out HadCM3 for climate estimates and BIOME4 as an ecological emulator. Studying environmental influences on global-scale hominin spatiotemporal dispersal and comparability to archaeological arrival points requires a high resolution. Here, the resolution of climate estimates from the CESM2, iLOVECLIM, and Paleo-PGEM, and reconstruction from ORCHIDEE are the most suitable. A combination of any of the three climatic reconstructions and the land-surface emulator ORCHIDEE would be the optimal solution for a palaeoenvironmental reconstruction applicable to the current research question. It should be noted that combining different data sets and emulators entails formatting data, which can be challenging depending on the systems in use. However, most systems offer support desks or have online communities, which can provide help and solutions.

When working with the Early and early Middle Pleistocene, caution should be paid to glacial and interglacial periods, which comprise substantially different climatic and environmental conditions (Clark et al. 2006; Herbert 2023; Lisiecki and Raymo 2005). This is important when establishing a habitat suitability map through niche modelling. Niche modelling has been used in a number of simulation studies to indicate population presence probability based on environmental conditions (i.e., Hughes et al. 2007; Romanowska et al. 2017; Yaworsky et al. 2024), and software tools have been developed to perform these calculations, further increasing the scope of application of the method (Elith et al. 2011; Merow, Smith and Silander Jr 2013; Muscarella et al. 2014). Since niche models provide static representations based on climatic and environmental data from a given point in time, it is essential to generate two independent habitat suitability maps; one projecting hominin presence in glacial periods and another for interglacial periods (see Figure 1). In this way, differences between glacial and interglacial periods, which may influence hominin dispersal, are highlighted.

The final step is defining how hominin populations react to the habitat suitability value of the cell. In previous modelling attempts of hominin dispersals out of Africa (Hughes et al. 2007; Mithen and Reed 2002; Nikitas and Nikita 2005; Romanowska et al. 2017), this reaction was tailored to the characteristics of the data and the research question. However, hominin reaction patterns overlapped with similar demographic aspects, such as the colonising rate, extinction rate, ecological preferences, and population growth, all of which depend on the environmental condition of the cell. Population growth is vital for the population’s large dispersal movements, when accompanied by a movement strategy, especially for a continental dispersal model. Therefore, the main aim of this process is to link climatic or environmental variables to population growth, so that each cell contains a value representing the suitability of the variable-defined environment for population prosperity. This approach has been applied in other studies to estimate population dynamics (Schmidt et al. 2020, 2025; Wren and Burke 2019), and is useful for studying past populations when empirical indications of population dynamics are limited (French 2016, 2021; Riede 2019). The output of this modelling approach thus provides a raster layer with each cell containing a value of population growth, based on the climatic/environmental variable chosen at the beginning of the process.

The above section presents the conceptual steps necessary for data and emulator selection, data formation, and linking environmental variables to population dynamics. The workflow presented here is merely one possible approach to linking environmental conditions with hominin dispersal behaviour. Depending on the data, scale, methodological preferences, and disciplinary backgrounds, other workflows can be chosen (Hölzchen et al. 2016). However, a workflow as presented here integrate approaches from multiple disciplines, providing a robust research foundation for addressing various concerns and aspects.