Physiography, foraging mobility, and the first peopling of Sahul

The study investigates how anatomically modern humans dispersed across Sahul (the landmass of Australia, New Guinea, and Tasmania during the late Quaternary), focusing on routes and speeds of migration. The context includes the necessity of open water crossings from Southeast Asia and rapid post-arrival dispersal across varied environments amid MIS 3 aridification. Prior research proposes multiple dispersal directions, coastal pathways, and interior corridors, often using dated sites interpolation, least-cost path models, or genetic analyses. The authors aim to quantify how evolving physiography and ecological constraints shaped hunter-gatherer mobility by integrating time-evolving landscape models with Lévy-walk-based mechanistic movement to simulate foraging and migration dynamics, and to compare predicted routes, site visitation likelihoods, and peopling speeds against archaeological data.

The paper reviews approaches used to study Sahul peopling: interpolation between dated archaeological sites with first-order environmental factors; agent-based least-cost path models combining ecological niche and local decision-making; genomic analyses providing dispersal insights but limited by data availability and temporal resolution; and newer methods including ecomorphological agent-based models with machine learning, time-varying climate coupled to stochastic-ecological models, reaction–diffusion methods, and directionally supervised cellular automata and pedestrian network models. Prior hypotheses include simultaneous multi-direction dispersal, preferential coastal routes, and interior corridors influenced by biogeography, climate, ecology, ethnography, demography, and technology. Ethnographic studies on hunter-gatherer mobility (residential vs logistical moves) and evidence for Lévy-like foraging inform the movement modelling framework.

Physiography modelling: The authors use goSPL, a parallel landscape evolution model, to simulate Sahul's evolving topography from 75 to 35 ka at 1 km input resolution (resampled) with outputs every 1 kyr. Forcing includes a relative sea-level curve and paleo-precipitation from the HadCM3 coupled climate model at 2 kyr increments, interpolated to the model grid. The model includes river incision via stream power law and hillslope diffusion; tectonic uplift is set to zero over the period. The simulation tracks erosion, deposition, drainage networks, and hydrology across Sahul, showing stable major drainage patterns, spatially variable denudation, and hydroclimatic trends consistent with proxy records.

Physiographic resistance layers: From model outputs, four layers are derived—slope index, environmental constraints (barriers such as endorheic lakes and oceans; major rivers as partial barriers; low costs along river corridors; high costs in arid flat regions), physiographic diversity (multiscale topographic position index combined via Shannon equitability), and stream position/water flux. Net primary productivity (LOVECLIM) and reduced costs along coasts are added. These are amalgamated into a normalized resistance map (0–1) representing mobility costs.

Mechanistic movement simulations: Using SiMRiv, a spatially explicit multistate Markov model, movements are simulated as two-state Lévy-like walkers switching between random and correlated random walks (turning angle concentration 0.95; transition probabilities 0.01 random→correlated, 0.002 correlated→random), step length up to 1 km, perceptual range 10 km (walker minimizes resistance within that radius). Resistance maps are updated every 1 kyr to capture physiographic changes.

Entry scenarios and runs: Two hypothesized entry points are used—a northern route via the Bird’s Head (Vogelkop) Peninsula (~73 ka) and a southern route via the Timor Sea shelf into northwest Australia (~75 ka). For each route, 5000 realizations are run with 10 million steps per realization (upper limit ~250 km/yr if fully utilized). Simulated paths are aggregated to compute kernel density estimates (KDE) of occupation likelihood and probabilities of presence.

Linking to archaeology: Forty archaeological sites are used to assess visits (a site is reached if within 10 km). Walkers passing near sites are grouped by cumulative traveled distance using k-means to infer distinct visits (occupation numbers). To extract residential mobility from composite foraging/migration movement, displacements are aggregated within radii of 25 km and 50 km to define single residential moves, yielding traveled distance estimates to each site for migration speed calculations. Migration speeds are computed a posteriori by dividing predicted traveled distances by archaeological age ranges, assuming first arrival at 75 ka (southern) and 73 ka (northern). Additional analyses include random hypothetical site assessment (n=5000) to map potential residential likelihood and preservation context.

• Simulations reproduce dispersal radiating across Sahul, with movement preferentially along riverine corridors and some coastal segments. KDE of occupation varies by entry: northern entry emphasizes routes on both sides of Lake Carpentaria (stronger east), while southern entry favors western interiors, south of Lake Carpentaria toward the east, and southern Western Australia.
• Archaeological site reachability: Best single realizations reach up to 39/40 sites for each entry. Across 5000 runs per entry, all sites are reached >30% of the time except Ivane Valley (Papua New Guinea highlands) <5%. Several sites exceed 60–70% visitation (e.g., Madjedbebe, Nauwalabila, Nawarla Gabarnmang; Ngarrabullgan Cave; Puritjarra; Riwi; Carpenter’s Gap 1). Average number of walkers near sites typically 100–500; inferred occupation numbers generally 10–30, with >50 at Nepean River sites.
• Traveled distances (residential-move radii): For northern route, minimum distances to sites are ~10^2–6×10^3 km (25 km radius) and 6×10^3–4×10^4 km (50 km radius). For southern route, minima are 4×10^2–3×10^3 km (25 km radius) and 3×10^2–2×10^4 km (50 km radius). Overall mean traveled distances are marginally lower for the southern route (2.32×10^5 and 1.91×10^5 km for 25 and 50 km radii) than for the northern route (2.47×10^5 and 1.94×10^5 km).
• Migration speeds (mean across sites): Southern entry 8.28 km/yr (25 km radius) and 5.95 km/yr (50 km radius); northern entry 9.88 km/yr (25 km) and 6.70 km/yr (50 km). Using minimum traveled distances yields much lower mean velocities: southern 0.36 and 0.25 km/yr; northern 1.15 and 0.83 km/yr.
• Probability-of-presence maps align with many proposed least-cost ‘superhighways’ but diverge in some areas (e.g., limited support for high-elevation New Guinea corridors and the Great Australian Bight coastal superhighway), reflecting physiographic costs, sparse paleo-rivers, and low NPP.
• Predictive archaeology: From 5000 random hypothetical sites, regions with high residential likelihood are identified across wetlands, river/lake banks in Pilbara and Western Desert, waterholes in arid zones, and incised valleys near Atherton Tablelands. Ten high-potential areas (sites A–J) are highlighted as targets for future archaeological surveys.

The results address the peopling question by demonstrating that evolving physiography and Lévy-like foraging can generate realistic dispersal patterns and site visitation consistent with archaeological records. Rather than a single optimal ‘superhighway’ or purely coastal route, the simulations suggest diffuse spreading guided by river systems and interior corridors, differing by entry route. Calculated migration speeds, derived a posteriori from traveled distances and site ages, fall within ranges reported for other regions and models, supporting plausibility of the inferred pathways. Comparison with least-cost and supervised cellular automata models shows broad concordance in main corridors but key differences where physiographic resistance and low productivity limit movement (e.g., high-elevation PNG, Great Australian Bight). Two methodological advances—explicit time-evolving landscape constraints and Lévy-walk foraging dynamics—improve realism over static-cost or Gaussian diffusion models, potentially refining expectations for traveled distances and timing of dispersal. The ability to estimate residential likelihood and propose new survey targets underscores the archaeological value of the approach.

By coupling a climate-driven, time-evolving physiography model with a mechanistic Lévy-walk movement framework, the study provides a unified, process-based view of early human dispersal across Sahul. It predicts interior-oriented migration along major drainage networks with route differences by entry point, reproduces visits to known archaeological sites, and yields migration speeds comparable to independent estimates. The framework also identifies areas of high residential likelihood to guide future fieldwork and offers explanations for the prevalence of small, shallow early sites consistent with repeated short-term occupations by highly mobile groups. The approach complements least-cost and cellular automata models and is transferable to other regions. Future work could integrate demographic dynamics, social-cultural decision rules, refined temporal treatment of movement, site taphonomy, and higher-resolution, locally calibrated environmental proxies to further constrain routes, speeds, and site preservation potential.

• Demographic processes and social, cultural, and economic decision-making are not explicitly modeled.
• Movement time is not simulated directly; velocities are inferred a posteriori from distances and archaeological dates.
• Entry scenarios are limited to two hypothesized routes and a 75–73 ka arrival framework; multiple waves and alternative entries are possible.
• Archaeological dates represent minimum peopling ages and may postdate first presence in a region.
• Taphonomic processes and preservation biases are not modeled; preservation potential is only indirectly inferred from landscape evolution outputs.
• Spatial and temporal resolutions (e.g., 2 km resistance grids, 1–2 kyr environmental updates) may smooth local variability; climate and sea-level reconstructions and model parameters entail uncertainties.
• High-elevation corridors and coastal shelves have uncertainties due to physiographic diversity and paleo-river reconstructions.