ARTICLE

Received 2 Jun 2016 | Accepted 20 Apr 2017 | Published 30 May 2017 DOI: 10.1038/ncomms15696 OPEN

Modelling the role of groundwater hydro-refugia

in East African hominin evolution and dispersal
M.O. Cuthbert1,2,3,4,5, T. Gleeson6, S.C. Reynolds7, M.R. Bennett7, A.C. Newton7, C.J. McCormack4
& G.M. Ashley8

Water is a fundamental resource, yet its spatiotemporal availability in East Africa is

poorly understood. This is the area where most hominin first occurrences are located,
and consequently the potential role of water in hominin evolution and dispersal remains
unresolved. Here, we show that hundreds of springs currently distributed across East Africa
could function as persistent groundwater hydro-refugia through orbital-scale climate cycles.
Groundwater buffers climate variability according to spatially variable groundwater response
times determined by geology and topography. Using an agent-based model, grounded on the
present day landscape, we show that groundwater availability would have been critical to
supporting isolated networks of hydro-refugia during dry periods when potable surface water
was scarce. This may have facilitated unexpected variations in isolation and dispersal
of hominin populations in the past. Our results therefore provide a new environmental
framework in which to understand how patterns of taxonomic diversity in hominins may have
developed.

1 School of Earth and Ocean Sciences, Cardiff University, Main Building, Park Pl, Cardiff CF10 3AT, UK. 2 Water Research Institute, Cardiff University, The Sir

Martin Evans Building, Museum Ave, Cardiff CF10 3AX, UK. 3 Department of Geography, University College London, Gower Street, London WC1E 6BT, UK.
4 School of Geography, Earth and Environmental Sciences, University of Birmingham, Birmingham B15 2TT, UK. 5 Connected Waters Initiative Research

Centre, University of New South Wales, Sydney, New South Wales 2052, Australia. 6 Department of Civil Engineering and School of Earth and Ocean
Sciences, University of Victoria, PO Box 1700, ECS 316, 3800 Finnerty Road, Victoria, British Colombia, Canada V8W 2Y2. 7 Institute for Studies of
Landscape and Human Evolution, Faculty of Science and Technology, Bournemouth University, Bournemouth BH12 5BB, UK. 8 Department of Earth &
Planetary Sciences, Rutgers University, Piscataway, New Jersey 08854-8066, USA. Correspondence and requests for materials should be addressed to
M.O.C. (email: cuthbertm2@cardiff.ac.uk).

NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696

E
stablishing the link between environmental change and sustain a small wetland—see Methods). In the present day,
resource availability in the East African Rift System (EARS) despite a likely mapping bias towards underestimation especially
is a matter of intense debate, the resolution of which is in the more humid and upland areas, over 450 such springs occur
fundamental to understanding hominin evolution as well as in the region (Fig. 1). Of these over 85% are fresh, and provide the
patterns of dispersal1. Most authorities currently suggest that only naturally available potable water year-round within semi-
climate change was a key factor for hominin evolution2–4 based arid and arid areas. In contrast, the majority of modern lakes and
on the coincidence of climate shifts with major evolutionary streams in the eastern branch of the EARS are either alkaline/
events such as the appearance and extinctions of hominin saline or ephemeral, with the exception of catchments that drain
species5–7. Some focus on the idea that evolutionary adaptation the extensive Ethiopian Highlands, an area of above average
was caused by multiple climate cycles and enhanced pulses of rainfall compared to the wider region. Of the lakes that are fresh
climate variability4,6. The apparent coincidence of enhanced lake (8 out of 34), all are small, with the exception of Lake Turkana
levels with evolutionary pulses is also considered significant by (Kenya) (Fig. 1). While brackish and relatively alkaline in the
others8,9. In all cases, large-scale climatic events are held to be present day, this lake nevertheless supports a fresh water fish
responsible for modification of local habitats and resource population and is therefore considered here as a ‘fresh’ lake,
distributions causing evolutionary consequences10, although the despite not being able to sustain modern humans on a sustained
detailed evolutionary mechanisms remain unclear. The variability basis. Notably, fresh water springs often occur within the
in the hydrological landscape is seen as having an important catchments of, and sometimes directly adjacent to less potable
role4,11. However, the potential for widespread groundwater saline/alkaline lakes, such as Lake Natron (Tanzania, Fig. 1).
hydro-refugia, such as springs and groundwater-fed perennial
streams, has long been neglected; yet it may challenge prevailing
views regarding the environmental context for hominin evolution Controls on the presence and persistence of springs. Within
and dispersal. We use this perspective to stimulate fresh thinking the drier parts of the study area where spring mapping is most
around the climate-forcing hypotheses, by focusing specifically on reliable, in comparison to what would be expected if springs were
how hydrological aspects of the landscape interact with climate randomly distributed across the landscape, there is a slight bias
change to control water availability, a key resource for survival. for fewer springs to be located in the most arid areas (Fig. 2).
In the EARS where most hominin first occurrences are However, there is also a counter-bias such that those springs
located12, potable water features in the form of surface water modelled as being most persistent are more likely to occur in the
which are persistent on greater than seasonal timescales, are driest areas. This runs counter to the intuition that catchments
scarce. For example, lakes are often alkaline, saline and thought receiving more groundwater recharge will have springs which
to have been increasingly ephemeral during the dry parts of persist for longer during dry periods and yields the surprising
precessionally forced climate cycles in the Plio-Pleistocene, within result that modern climate is not the primary control on spring
critical periods for hominin survival5,13,14. Present day conditions persistence in this context. This is explained by spring persistence
in much of the EARS are analogous to relatively dry periods being strongly correlated with groundwater response time,
[B70% is arid to semi-arid with groundwater recharge of which is a function of subsurface hydraulic properties (storage
o50 mm y  1 (ref. 15)] and therefore provide a way of exploring and transmissivity) and geometry (length scale and topographic
the likely hydrological conditions experienced by early hominins. gradient) (see Methods, Supplementary Table 1). In general
Here, we show how hydrogeological modelling of the present terms, the groundwater response time is a measure of how long
landscape coupled with agent-based modelling of hominin an aquifer takes to respond to a change in boundary conditions,
movement yields new insight into potential correlates of hominin such as rates of groundwater recharge varying due to climate
survival and dispersal. Because groundwater acts to buffer climate change. Therefore, given the same variations in groundwater
variability it could have provided, via springs and baseflow to recharge through time, spring discharge from an aquifer with a
perennial streams, hydro-refugia which persisted through long large groundwater response time will be relatively temporally
dry periods. In past transitions to wetter periods, trans-rift stable in comparison to spring flows issuing from an aquifer with
dispersal routes may have become active before those along the a much shorter groundwater response time. In wetter areas, while
rift axis, and under the wettest scenarios modelled hominin more groundwater recharge may be available to eventually dis-
dispersal (and therefore gene flow) may have been widely possible charge at springs, the resulting higher water table leads to more
across the region. The hydro-refugia model shows that early intersection between the water table and the topography. This
hominins, and later Homo, survival and dispersal is likely to have reduces the distances between points of groundwater discharge,
been facilitated under drier conditions than previously thought greatly reduces the value of groundwater response time, and thus
possible. increases the responsiveness of springs to variations in climate.
The persistence of a spring during a dry period is thus a complex
function of the timescale of climate variability, the topography
Results which determines the catchment of the spring, and the hydraulic
Present day fresh water distribution in East Africa. Our properties of the aquifer, all of which determine the temporal
hydrological mapping and modelling focused on groundwater relationship between the recharge input and the spring discharge.
manifested as springs but also assumed that major regional rivers, The models we have used integrate all these factors to determine
for which there is geological evidence of persistence through dry the primary controls.
periods16, were groundwater fed. Here, we therefore define both Whereas previous research has focussed on climate variability
these types of hydrological features as potential groundwater being the dominant control on the availability of water, the data
hydro-refugia. We quantified the distribution and persistence of and models presented here show that geology and topography act
springs by first mapping existing ‘permanent’ springs through the to greatly buffer the impact of climate variability. For example,
eastern branch of the EARS, then modelling the temporal our results indicate that a majority of springs (B40–60%) would
persistence of each spring above a flow threshold required for a remain productive for periods of hundreds of years (Fig. 3)
spring to act as a significant water source, within an otherwise and around 30% (that is, 4100 springs) would still remain as
dry landscape (that is, 1,000 m3 y  1—enough running water to hydro-refugia in even drier parts of precessional cycles, assuming
provide drinking requirements for hundreds of animals and to transitions from modern recharge conditions to conditions of just

2 NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696 ARTICLE

100
km 90

Cumulative % of active springs in areas

80

receiving <60 mm y–1 recharge

70

60
All springs
50 P = 300 y
P = 3000 y
40
P = 23000 y

is
30

ax
Regional recharge rate

y
le 20
al

Semiarid
ETHIOPIA
v

Humid
ft
Ri

10

Arid
0
0 10 20 30 40 50 60
Recharge (mm y–1)

Figure 2 | Cumulative frequency distribution of springs as a function of

Lake the average catchment recharge compared with the regional recharge
Abaya distribution in areas receiving o60 mm y  1 recharge. Springs were
considered as ‘active’ if they maintain flow 41,000 m3 y  1 for gradual
recharge variations sinusoidally fluctuating from zero to the average
Lake catchment recharge with period (P). The ‘regional recharge rate’ curve is
Turkana indicative of the cumulative frequency distribution that would be expected if
springs were randomly distributed across the landscape. The deviations
between the actual spring distributions indicate that spring persistence
(as opposed to spring presence) is not controlled primarily by modern
KENYA climate, consistent with statistical results which show that spring
persistence is not significantly correlated to any individual catchment
characteristic including groundwater recharge (Supplementary Table 1).

1 mm y  1 recharge (that is, arid-hyperarid). For a sudden

hydrological transition to arid conditions, we see an initially
quicker decline in the number of springs still remaining
productive, as would be expected (Fig. 3 and Supplementary
Fig. 1). Sensitivity analysis and Monte Carlo experiments indicate
N that while the hydraulic parameter uncertainty for our estimates
of spring persistence for an average spring is ±37%, the
Lake
maximum error in the combined modelled percentage of springs
Natron
persistent on any timescale is only ±5% (Fig. 4, Supplementary
Hominin sites Figs 2 and 3).
Main rivers
Olduvai Perennial streams The abundance of persistent groundwater hydro-refugia. There
Gorge Ephemeral streams are only a small number of hydro-refugia which are known from
Persistent springs the geological record to have survived during the driest periods of
TANZANIA
All other springs precessional climate cycles and are also associated with hominin
fossils and stone tools. Lake Turkana c. 2–1.85 My16 (Kenya,
Figure 1 | The pattern and mode of hydrologically available water in Fig. 1) is an important example, supported by the paleo-River
present day eastern Africa. The distribution of water is controlled by Omo, which has a catchment in the Ethiopian Highlands.
geology, topography and climate. Hominin sites are closely associated with Although there is an increasing recognition of paleo-springs in
the rift valley axis. Fresh water lakes (dark blue), alkaline/saline lakes (red), the geological record, they are still relatively rare. This is likely
wetlands (pink), background is groundwater recharge from Döll & Fiedler15 due to their low preservation potential; groundwater discharge
coloured yellow (250 mm y  1) to white (0 mm y  1). ‘Persistent Springs’ only leaves a direct geological record under specific geochemical
are those modelled as productive (41,000 m3 y  1) at precessional (23 ky) conditions that lead to mineral precipitation (such as tufa17).
minima under gradual climate change. The number of such springs is However, a spring located at paleo- Lake Olduvai c. 1.84 My17
considered conservative, since at least some persistent springs are likely to (Tanzania, Fig. 1) is thought to have continued to flow during a
be present during dry periods in areas currently mapped as having perennial prolonged period of aridity18 at the precessional minimum. Our
streams. Streams, lakes and marshes digitized from map series as described results suggest that there would have been orders of magnitude
in the Methods; National borders and main rivers 45 km3 y  1 from more groundwater hydro-refugia during arid phases than the
‘Major rivers of the World, classified by mean annual discharge’, GRDC geological record currently indicates. Furthermore, groundwater
http://grdc.bafg.de. Projection: WGS 1984. hydro-refugia would have been much more abundant than fresh
water potable lakes during periods of prolonged aridity. Since
some recharge still occurs even in arid areas in the present

NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696

100 0

an
n l /IOD
sp
90

ife
Hu NSO
Fig 5c 0.2

ma
E
80

Gr (p

Cumulative frequency
ad eri
70

ua od
0.4

l tr ic)
% of springs still active

an
Su

sit
60
dd (ste

i on
en p)
tra
0.6
50
ns
itio
Fig 5a
40 n
0.8
30

cyc l &

al
les
DO lenia

bit
Or
20 l 1
Mi
100 102 104
10 Groundwater response time (y)

Figure 4 | Cumulative frequency distribution of groundwater response

0
times for 10,000 Monte Carlo realizations indicates that the combined
0 3 30 300 3,000 23,000
parameter uncertainty in the CDF of modelled spring longevity across the
Period of oscillation or time needed for spring
discharge to fall below 1,000 m3 y–1 (y) EARS is ±5%. The result using an average set of parameters (plotted in
red) is superimposed on all other realizations (plotted in black). The
Figure 3 | Many groundwater hydro-refugia persist through even maximum spread occurs at the steepest part of the curve for groundwater
precessional orbital climate cycles. Modelled frequency of springs that are response times of 100–200 years. Since groundwater response times vary
still active (that is, flow rate still above 1,000 m3 y  1) at the driest part of by five orders of magnitude, changes in groundwater response times can
climate changes over a range of timescales under gradual or sudden have more impact on spring persistence than changes in recharge.
boundary condition changes (ENSO—El Nino—Southern Oscillation; However, only a small amount of recharge is needed to maintain flow to a
IOD—Indian Ocean Dipole; DO—Dansgaard–Oeschger). Recharge is spring; for example, a recharge of just 1 mm y  1 over an area of 1 km2
assumed to vary between the modern day value to a minimum value of could maintain flows of 1,000 m3 y  1. Therefore, while a small amount of
1 mm y  1. Sensitivity of persistence to the recharge occurring at the driest recharge is a necessary, but not always sufficient, condition for a persistent
point in a climate cycle is explored in Supplementary Fig. 1. spring, the geological and topographical characteristics of the EARS are also
fundamental controls.
day15,19–21, the 1 mm y  1 recharge applied during the driest
periods used in our models leads to a conservatively low number
of predicted active springs during these periods. The fact that at system’s constituent units (the agents) and their interactions with
least some persistent springs are likely to have been present the terrain. The landscape is heterogeneous in terms of land
during past dry periods in areas mapped as having perennial cover, slope and roughness, and this affects the ability of agents to
streams in the present day also suggests our predictions are traverse an area. This is explicitly incorporated by including
conservative (Supplementary Fig. 4). Nevertheless, while the variation in transit time in response to terrain properties within
presence of spring-based refugia has been demonstrated for an the model. There are many potential routes between two points in
isolated site22,23 our results suggest their presence and a landscape, and our model allows agents to explore a wide range
importance over a much wider geographical area. of different routes, according to the decisions that would be made
by an individual walking across a landscape. Only those
successful transits between water bodies are considered in our
Modelling the dispersal of hominins. We examine the evolu- analysis as evidence of linkage between ‘nodes’ (or water bodies).
tionary implications of this improved hydrological understanding We have modelled a range of climate scenarios along the
using an agent-based model (ABM). The specific advantages continuum of a precessional orbital climate cycle. In the
of using an ABM in this type of context are summarized by absence of specific data on relative spatial changes in effective
Bonabeau24 who identifies three particular benefits compared to precipitation through a precessional cycle (for example,
other modelling approaches, namely their ability to: (i) capture regionally, or say between highlands and lowlands), it is assumed
emergent phenomena; (ii) provide a natural description of a that as the climate becomes more arid, decreases in effective
system; and (iii) be relatively flexible. We apply the principal rainfall and groundwater recharge occur proportionally every-
of uniformitarianism and consequently use the contemporary where leading to progressive shrinking and fragmentation of the
landscape as a method of exploring potential constraints on hydrological network. Conversely as rainfall/recharge increases,
hominin movement imposed by the availability of water (that is, we assume that water tables increasingly intersect stream and lake
springs, streams, lakes or wetlands). Our aim is to understand the beds and there is enhanced potential for springs to occur at
role that changes in the generic hydrological resource network, higher elevations, all of which leads to an expanded hydrological
through a simulated climate cycle, may have played as one network. We recognize four hydrological components in our
control on hominin movement. Here, we do not use the ABM to hominin mobility model: (1) mapped springs (seasonal, perennial
examine explicit hypothesizes related to past events or scenarios, and geothermal); (2) mapped streams (seasonal or perennial),
rather we examine generic principles that might apply in this and major rivers using flow thresholds set at 40 or 45 km3 y  1
past landscapes. The connectivity between two water sources in (GRDC data, see Methods); (3) mapped wetlands (seasonal,
our study does not depend on fixed properties based on their perennial); and (4) mapped lakes (fresh, saline and seasonal).
location, but arises as an emergent property of the simulated Four hydrological scenarios were run representing a dry to wet

4 NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696 ARTICLE

Table 1 | Climatic scenarios modelled using the agent-based model.

Run Id Climate Hydrological components

state
Run-1 Future wet Modern perennial/seasonal springs, geothermal springs, fresh water and saline lakes, perennial/seasonal rivers, major rivers
(flow 40 km3 y  1) and perennial/seasonal wetland/marsh.
Run-2 Present wet Modern perennial springs, geothermal springs, fresh water and saline lakes, perennial rivers, major rivers (flow 40 km3 y  1) and
perennial wetland/marsh.
Run-3 Present dry Modern perennial springs, geothermal springs, fresh water and saline lakes, major rivers (flow 40 km3 y  1) and perennial
wetland/marsh.
Run-4 Future dry Persistent springs (that is, those modelled as being persistent through precessional cycles), geothermal springs, fresh water lakes
and major rivers (flow45 km3 y  1) with saline lakes excluded as either desiccated or too hyper-saline to be potable.

continuum (Table 1). Run-4 corresponds to the driest phases prevailing climate would result in decreased streamflow and a less
of a precessional climate cycle (for example, 23 ky), where expansive stream network than observed in the present day even
predominantly arid conditions existed across eastern Africa, with during periods of relative (for example, seasonal) wetness.
only the most persistent springs, deepest lakes and highest order
and/or groundwater-fed sections of the stream network present.
Discussion
One of the fundamental problems in explaining human evolution
Controls on hominin mobility and gene flow. The modelling by invoking climate in East Africa is that precessional climate
results demonstrate: (1) That some hominin movement may have forcing occurs on a timescale that appears too short for allopatric
been possible between spring networks and along major rivers speciation to occur. To counter this, one has to invoke the
(groundwater hydro-refugia) that would have allowed hominin variability itself as the key factor4 or look to longer eccentricity
populations to survive in specific regions, even during the most cycles and their modulation of the precessional cycle amplitude5.
extreme arid climate phases (Figs 5a and 6a). The presence Our work suggests three alternative possibilities, however:
of hydro-refugia during the driest of scenarios is robust taking (1) a different geographical distribution of landscape elements
into account uncertainty in the model input parameters at various times in the past may have favoured longer term
(Supplementary Fig. 6). (2) That under the ‘present dry’ condi- periods of isolation; (2) that the population density was such that
tions modelled hominin mobility occurs transverse to the rift axis while connectivity was possible it was not exploited; or, (3) more
rather than along it (Figs 5b and 6b). In specific cases modelled, likely in our view, that climate may not play such a primary
modern springs high on the rift margins play a role in connecting role in human evolution, as is commonly asserted. In fact, the
the rift floor hydrological system with those of the rift flanks. potential for frequent, widespread dispersal as illustrated here,
Note the presence of cross-rift connections in the Ethiopian and the resulting potential episodes of genetic admixtures, might
Highlands is sensitive to the input parameters used in the model explain the lack of phylogenetic diversity in the hominin lineage
(Supplementary Fig. 7) but still holds for other areas of the rift noted by some32.
(Fig. 5b and Supplementary Fig. 8). The importance of dispersal In conclusion, the hydro-refugia model therefore points to the
routes transverse to the rift axis is at odds with the common need to evaluate a range of parameters and variables beyond
assumption of along-rift dispersal23,25 but shows agreement the current paradigm of climate-driven environmental change
with westward dispersals observed in some hominins26 and to explain hominin evolution. It provides a basis for palaeoan-
genetic studies of several other species27,28. (3) Under the wetter thropology to explore the possible mechanisms by which
scenarios modelled (Figs 5c and 6c and Supplementary Fig. 8), in taxonomic diversity in hominins arose32, and may help to
which the fluvial network becomes dominant, the potential for explain levels of genetic exchange identified in ancient African
widespread hominin dispersal as well as associated gene flow is populations33. Our hydrological results have wide global
evident. These conclusions are summarized in Fig. 7 and have applicability (drylands cover around 45% of the Earth’s
been found to be statistically robust through model repetition (see landmass34) and the importance of groundwater for the
Methods; Supplementary Table 2 and Supplementary Fig. 15). survival of our hominin ancestors when faced with dramatic
In addition to the long timescale changes in climate expected climate changes in the past could also inspire and inform
through a precessional cycle, shorter term variations (for strategies for human resilience to future climate change35,36.
example, seasonal dry periods or multi-year droughts) would
have altered the availability of fresh water. For the ‘present’ Methods
scenarios, the effect of seasonality on the potential connectivity Data compilation and hydrological mapping. The study covers an area of
of hydro-refugia is incorporated in the analysis by comparing 2,093,280 km2 stretching from northern Tanzania (  5.390°S, 33.996°E) to
the ‘present wet’ and ‘present dry’ scenarios (Fig. 5b and Ethiopia (13.797°N, 42.799°E) and focused along the eastern African Rift.
Supplementary Fig. 10). Such patterns of variation in the location Taxonomic occurrences of early (pre Homo sapiens) Hominidae in Ethiopia,
Kenya and Tanzania, representing 178 sites in total, were collated from all available
of available spring water are in accordance with the modern records in the Paleobiology Database (Fossilworks: http://fossilworks.org),
experience of East African communities29,30. For the ‘future wet’ and their co-ordinates digitized.
scenario, mobility is already so easy that seasonality would make Stream locations were digitized from 1:500,000 British Army maps dating from
little difference (Fig. 5c). For the ‘future dry’ scenario, the impact the 1940s with gaps in-filled from: (1) modern 1:250,000 topographic maps,
Kenyan Government, 1981; (2) East Africa, Series 1501, Joint Operations Graphic
of seasonality is harder to constrain but in the driest parts of the (Air) 1:250,000 US Air Force; or (3) East Africa 1:250,000, Series Y503, British
precessional cycle envisaged, seasonal expansion of the drainage Overseas Mapping circa 1963 (http://www.lib.utexas.edu/maps/ams/east_africa).
network is likely to have been much less than that during Second-order rivers were digitized for both seasonal and non-seasonal cases.
the present day. This is because both runoff and recharge are Modelled flow data was used to define major perennial rivers as published by
GRDC37, based on WaterGap 2.1 global hydrological model output. Lake outlines
strongly controlled by antecedent moisture and water table were obtained from Digital Chart of the World (http://www.diva-gis.org/gdata)
conditions21,31. Hence while seasonal mobility may be enhanced with salinity for named lakes determined via literature. If salinity information was
to some extent even during short wet periods, a sustained drier not available then the lake was classed as ‘fresh’ by default.

NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696

a b c

N
Relief Path density
5,884
High

0 0 140 280 km
Low
Persistent springs Geothermal springs Perennial springs Geothermal springs

Figure 5 | Modelled dispersal of hominins between hydro-refugia through a wet–dry climate cycle. As the climate changes from the driest (a) to wettest
(c) state envisaged during a precessional cycle, the modelled dispersal of hominins between hydro-refugia increases until dispersal is almost ubiquitous.
Agent-based modelling results based on three hydrological scenarios, using a maximum three-day travel distance of 150 km and surface roughness as the
cost layer scaled according to Supplementary Fig. 11. The black lines shown represent the tracks of agents in the model. It is important to note that just
because the model shows a potential pathway between two water sources it does not necessarily mean that it was actually followed. (a) Driest scenario
using persistent springs (Fig. 1), geothermal springs, fresh water lakes, major rivers with a flow 45 km3 y  1 (Run-4). Note the networks of springs acting
as hydro-refugia, which persist irrespective of the topographic cost layer and scaling used in the model (Supplementary Figs 9–12). (b) Dry scenario using
modern springs (perennial þ geothermal), lakes (fresh þ saline), perennial wetlands and major rivers with a flow 40 km3 y  1 (Run-3). Note the potential
for dispersal transverse to the rift axis (white dashed line) and the absence of along-axis dispersal routes. Springs act to connect rivers on the rift flank with
those of the rift floor. The cross-rift movement in the vicinity of Lake Abaya (Ethiopia), where the rift cuts the Ethiopian Highlands, is sensitive to the
topographic cost layer and scaling used in the model (Supplementary Fig. 9); however, the principle holds in southern Kenya. The next stage in the
increased water availability continuum is shown in Supplementary Fig. 10, which includes perennial streams and reveals a progressive increase in cross-rift
movement and the start of dispersal along the rift (Run-2). (c) Wettest scenario uses modern springs (seasonal, perennial þ geothermal), wetlands
(perennial þ seasonal), lakes (fresh þ saline), major rivers with a flow 40 km3 y  1 (Run-1). Note the potential for widespread dispersal of hominins and
genes. Background relief map is based on a 30  30 m SRTM model: http://earthexplorer.usgs.gov/. Projection: WGS 1984.

The location of the main East African Rift Valley axes were digitized from We are also aware of possible bias of maps underestimating the number of springs
Hayes et al.38 Locations of springs were digitized from a variety of available map in the wettest areas of the study area, where perennial streams are frequent and
sources to enable 100% coverage for the study area on as consistent a basis as springs may thus be overlooked. This is consistent with known springs used for
possible as follows: Ethiopia (East Africa Y401—GSGS 4355, 1:500,000 scale), water supply in the Ethiopian Highlands (for example, Calow et al.29) which were
Kenya (East Africa JOG 1501 AIR, Y503, 1:250,000 scale; Tanzania (East Africa not represented on the large-scale maps. Given the likely underestimation of the
Y401—GSGS 4355, 1:500,000 scale and Geological Map 1st Edition Quarter Degree total number of springs, to ensure the sample of springs we have mapped is
Sheets, 1:125,000 scale). Mapping at smaller scales inevitably leads to a larger nevertheless representative for making inferences regarding spring persistence, we
number of mapped springs. A comparison of map sheets from Tanzania where we have randomly resampled the total set of springs for increasing subsample sizes as a
have two scales of maps to compare indicates that the 1:500,000 maps record, on proportion of the total spring set. For each subsample, we calculated the RMSE for
average, 40% of the number of springs which are present on the 1:125,000 maps. the cumulative frequency distribution of groundwater response times of the
The 1:125,000 maps have been ‘ground-truthed’ by the authors (unpublished) subsample against the full spring set. The plot of RMSE versus proportion of the
across several map sheets in northern Tanzania with the input of local Masai total sample size (Supplementary Fig. 4) demonstrates that the set of springs we
guides, suggesting they are remarkably accurate in representing the main sources of have mapped is representative as the errors reduce to zero for a subsample size of
water used by local people. However, the scale of mapping across the modelled area around 75% of the full set.
generally decreased in resolution from 1:125,000 in Tanzania, to 1:250,000 in ArcGIS was used to delineate for each spring: the distance to the nearest
Kenya, and 1:500,000 in Ethiopia. Thus, we expect that the undersampling due to watershed boundary (B) average slope (a) using the Hydrosheds database39 and a
the changing map scales leads to a loss of accuracy in the absolute spring count 90 m spatial resolution Digital Elevation Model40. The Hydrosheds database was
from close to 100% accuracy in the south of the area to around 40% in the north. also used to derive upstream contributing catchment areas (A). Hydraulic

6 NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696 ARTICLE

conductivity (k0) and porosity (ne) values for each spring location were sourced uncertainties involved in modelling long timescale climatic changes in this context.
from Gleeson et al.41, which is based on the highest resolution mapping digitally Thus our choice of modelling approach needed to be generic enough to cover the
available42. While, we recognize that there may be local variations in groundwater most important features of groundwater flow systems across a range of contexts in
recharge43–45, we used the distribution of potential groundwater recharge from which EARS springs are found, while being mathematically simple enough to
Döll and Fiedler15 consistent with other regional African groundwater studies42. apply existing analytical solutions to the governing flow equations to enable the
The range of input parameter values are plotted in Supplementary Fig. 2. uncertainties to be easily explored. Although there is debate about the distribution
ArcGIS was used to analyse relationships between spring persistence, climate and precise role of transverse faulting in controlling groundwater flow in the
conditions and catchment properties. EARS46, an emerging pattern is one of nested flow patterns whereby groundwater
age increases from rift flank to graben46–48. Due to the heterogeneity and
anisotropy of many EARS lithologies which in some locations determines spring
Development of models of spring persistence. East Africa is a very diverse locations, it would be virtually impossible to predict the locations of all springs on
but under-researched hydrogeological environment. There are also inherent the basis of regional scale geological mapping and analytical flow models alone.
However, by mapping the locations of the present day springs directly, it is then
reasonable to apply a simple flow model and bulk hydraulic parameters to
15 determine the likely variability of the springs to climate variations.
a Each spring was therefore modelled using a one-dimensional (1D) linearized
Boussinnesq equation as follows:

@Z k0 Z0 cos a @ 2 Z k0 sin a @Z
¼ þ ; ð1Þ
@t ne @x2 ne @x

10 where t is time [T], x is distance along the aquifer base [L], Z is hydraulic head [L],
k0 is aquifer hydraulic conductivity, ne is aquifer specific yield [  ], a is the slope of
the aquifer [  ], Z0 is the average water table height equal to pD, where D is the
maximum saturated thickness of the aquifer [L] and p is a linearization constant
normally assumed to equal 0.3.
This approach assumes that a homogeneous, isotropic, sloping aquifer
extending to a watershed boundary at one end (no flow boundary condition at
5 x ¼ Bx ¼ B/cos(a)) receives groundwater recharge uniformly over its surface and
transmits groundwater discharge to a spring at its lower end (Dirichlet boundary
condition at x ¼ 0). In reality, preferential flow through high-permeability fracture
features may dominate hydraulic head distributions at a local scale in some spring
systems. However, it has been shown that even in such cases, similarly idealized
analytical spring flow models to those used here can correctly simulate the
observed spring discharge dynamics of the bulk groundwater flow system49.
0 Our models neglect heterogeneity and anisotropy and assume isotropic hydraulic
–2 0 2 4 6 8 10 12 14 16 18 parameters—hence, they could be improved in the future if higher resolution maps
of the required hydrogeological parameters become available for the region.
30 Two scenarios were modelled for ‘sudden’ and ‘gradual’ climate change end
b members as follows.
25 First, models for sudden climate change were implemented assuming recharge
(R) ceases entirely after a period of climatic steady state. Following Brutsaert50,
spring discharge (q) is given by:
20
     2 Hi2  
X1
zn2 1  2 cosðzn Þexp Hi exp  zn þ 4 t þ
qðt Þ ¼  2Bx Rcos a  2
2   2 ;
15 n¼1;2;3 ... zn2 þ Hi4 þ Hi2 zn2 þ Hi4
ð2Þ
10
with zn being the nth root of tan(z) ¼  2z/Hi and:

5 Bx tan a
Hi ¼ ; ð3Þ
Z0
0
 
k0 Z0 cos a
tþ ¼ t: ð4Þ
–5 ne B2x
–5 0 5 10 15 20 25 30 35 40
Second, models for gradual climate change were implemented using
35 sinusoidally varying recharge as a top boundary condition as follows:
c
R ¼ Rav ð1  cosðot ÞÞ; ð5Þ
30
with recharge thus varying between zero (that is, hyperarid conditions) and a
maximum value (that is, 2Rav) for periods (P ¼ 2p/o) ranging from 1 y to 23 ky.
25

Figure 6 | Network shapes derived from successful journey matrices

20
output from the PATH model across the modelled climatic cycle based on
ten repeated runs. (a) Future wet scenario (Run-1). (b) Present dry
15
scenario (Run-3). (c) Dry (23 ka) scenario (Run-4). Note the progressive
increase in the number of sub-networks. The axes are dimensionless and
10
are expressed in non-geographic units. The networks simply portray the
relationship between nodes linked by common edges (that is, one or more
5 successful journey). Each sub-graph represents a distinct network that is
unlinked to any other by an edge; the more the sub-graphs, the poorer the
0 connectivity between the sum of the nodes present. The networks are
derived from matrix of successful/unsuccessful journeys and were created
–5 by exporting the successful journey matrix from the PATH model and
–5 0 5 10 15 20 25 30 35 40
plotting nodes and edges within Matlab.

NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696

a Drier b Dry c Wet

2
2

1 Connecting spring
2 Hydro-refugia network
2
3 Rift-flank river
Dry spring
Active spring
River

Cross Rift - No Cross Rift - Yes Cross Rift - Yes

Along Rift - No Along Rift - No Along Rift - Yes
[Isolation] [Dispersal + Isolation] [Dispersal]

Figure 7 | Conceptual model showing the role of springs across various climate scenarios. Under the drier scenario (a) hominin survival is focused on
single springs (or spring clusters) and movement between springs (or spring clusters) is limited. As climate improves (b) the availability of water sources
increases particularly as the water table intersects rift flank rivers. Springs high on the rift sides may act to link rift flank rivers with water sources in the rift
facilitating transverse rift movement. As water becomes widely available (c) hominin movement occurs in all directions including along the rift axis.

The solution for this case was derived by Cuthbert & Ashley22 as: and for the sudden end member as the time taken for spring flow to fall below this
X1   threshold following the step reduction of recharge from steady-state conditions.
A A The spring discharge threshold was defined, conservatively, as 1,000 m3 y  1. For
qðt Þ ¼ Rav 2 2
ð  o sinðot Þ þ Ccosðot ÞÞ  ; ð6Þ
n¼1;2;3 ...
C þo C example, this could be envisaged as enough running water to provide drinking
requirements for 100s of animals of B2,000 l d  1 (B750 m3 y  1) with the rest of
"    #
2k0 Z0 cos a zn2 1  2 cosðzn Þexp Hi the discharge sustaining a small wetland of B100 m2 transpiring at an average
A¼   2  2 ; ð7Þ annual rate of B2,500 mm y  1 typical of (semi)arid environments
ne B x zn2 þ Hi4 þ Hi
2 ( ¼ 250 m3 y  1).
 Springs likely appeared or disappeared due to local geological or topographic
k0 Z0 cos a 2 Hi2 changes, such as faulting and volcanism, or blockage of groundwater discharge by
C¼  z n þ : ð8Þ
ne B2x 4 sedimentation or accumulation of spring precipitates. However, we are examining
We implemented these equations in MATLAB setting the necessary numerical the combined regional distribution of springs through time and space in the EARS,
tolerances to yield an error of o1% in the resulting spring persistence metrics we a region which is still undergoing extensional tectonics that began B30 Ma ago52.
output. The models assume a 1D geometry with the output, q, being a discharge per It has therefore, throughout the time period in which hominins evolved, provided a
unit width of aquifer [L2T  1]). We have therefore rescaled the results to suitable hydrological setting by creating the necessary relief and landscapes for
approximate the actual spring flow (Q, [L3T  1]) for each entire spring watershed developing active groundwater flow systems; topographic highs that trap moisture
whereby Q ¼ A*q/B. In reality, the catchment geometry will often include non- to provide groundwater recharge and drive groundwater flow by gravity to
uniform (predominantly convergent) flow fields, which will affect the groundwater discharge in topographic lows53.
response times as outlined below.
Available recharge values from Döll and Fiedler15 were derived from a global Testing models of spring persistence. It was not possible to calibrate the models
hydrological model which takes no account of the underlying geology. It was directly across the region due to the sparsity of available spring flow time-series
therefore important to ascertain whether this recharge could be accommodated by data in East Africa and the long timescales considered in this paper. However,
the groundwater flow systems feeding each spring by ensuring that the modelled using what data were available, a range of model tests were carried out over
water table stays below the ground surface. To do this a steady-state solution to timescales from years to millennia as follows. One published multi-decadal spring
equation 1 was used for the previously stated boundary conditions as follows51: flow time series was found within the study area (Mzima Springs, Kenya54) and
 i 
tanðaÞx tanðiÞ h
tanðaÞx
 1 digitized along with the next nearest available long record from South Africa
Z ¼ RpD cosðaÞ 1  e pD  1  e pD Bx þ x k sin2 ðaÞ : ð9Þ (Uitenhage Spring55) with a comparable fractured-rock hydrogeology. For these
pD
sites, the rainfall time series (published with the spring flow records) were
Further, the position of the maximum water table height is given as follows:

 characterized using a superposition of two sinusoids whose relative amplitudes
pD ln 1 þ tanðaÞL were used as input data to the periodic model equations (6–8) to generate spring
pD
xmax ¼ ð10Þ flow output. For Mzima, hydraulic parameters were used directly from the GIS
tanðaÞ mapping described above, and for Uitenhage, the parameters were set using the
These equations were used to calculate a refined estimate of the recharge values published range of literature values for that site. The results for both springs show
needed to ensure physically realistic water table conditions. These refined values that the model performs well with respect to simulating the observed degree of lag
were less than or equal to those derived by the global hydrological model. and attenuation between the rainfall and spring flow (Supplementary Fig. 5) with a
The gradual and sudden end member climate scenarios were both modelled slight tendency for over-responsiveness in the model. The total spring flow for
using a range of maximum recharge values including the modern day distribution Mzima is underestimated by nearly 30%, but since it is based on uncalibrated
of calculated actual recharge, as well as with spatially constant values from the parameters defined on the basis of regional GIS mapping, this is a very
current regional potential recharge average (49 mm y  1). Recharge minima were reasonable result.
defined as a 1 mm y  1 rate considered to be a conservative value for arid The only way of directly testing the model validity over much longer timescales
conditions 15,19–21, and scenarios with minima of 0 mm y  1 recharge were also run is through geological evidence, and output from an identical analytical model has
as the most extreme possible end member. We assigned hydraulic properties to the previously been shown to be consistent with the geological record at Olduvai
model based on the results from mapping as described above, assuming D ¼ 100 m Gorge, Tanzania within the study area22. However, another approach is to compare
and that specific yield was equal to drainable porosity, and tested the parameter the analytical model used here with a more complex numerical model simulation
uncertainty as described below. operating over longer timescales, as a proxy for actual data. Such a model
Spring persistence was defined for the gradual end member as the maximum developed for North Africa56 was successfully calibrated to present day data and
period of variation for which spring discharge did not fall below a threshold value, then hindcasted to simulate a 20 ka period of groundwater discharge recession with

8 NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696 ARTICLE

no groundwater recharge since the last ‘pluvial’ period. Our analytical cessation cell size is 1 by 1 km. The terrain is rarely flat and each cell traversed has a cost
equations (2–4) have therefore been applied using the range of parameters from related to the topography which consumes the agent’s energy proportionately
the published ‘base case’ model and the outputs compared (Supplementary Fig. 5). thereby reducing the total travel distance before the agent’s death. Supplementary
The recession from the published study has an initially steeper recession than the Fig. 9b shows the variation in network connectivity with different travel distances
analytical model, but yields a similar overall recession on the 20 ka timescale and also the intra-variability at any chosen distance (see also Supplementary
modelled. Fig. 10c). The inter-run variability is small compared to that between different
Despite the lack of data for calibration, these tests thus demonstrate the travel distances. Other parameters included in the model are ‘maximum turn
reasonableness of modelling approach applied for the aims of quantifying, to angle’. Each agent places its back to the resource-patch from which it is hatched
within an order of magnitude accuracy, the response times of springs in East Africa and moves forward. If the turn angle is low then there is little deviation in agent
over a range of timescales from years to millennia (Supplementary Fig. 5). motion to the left or right. As the turn angle increases there is greater capacity for
lateral movement and in theory for an agent to return to its source.
Experimentation revealed that an angle of 20° gave the best compromise in terms
Groundwater response times and uniform flow field assumption. In general of exploration to the right and left versus directed forward motion. The results are
terms, the groundwater response time (GRT) is a measure of how long an aquifer not sensitive to turn angle set unless it exceeds about 60° when a significant
takes to respond to a change in boundary conditions, such as rates of groundwater number of agents begin to lose their way and return home.
recharge varying due to climate change. It is given by the equation GRT ¼ B2/D Other model inputs subject to parameter choices include the cost-map and
with hydraulic diffusivity for a sloping aquifer, D ¼ T/(necos(a)) (ref. 50), and cost-scaling. The impact of different cost-maps, variations of slope and roughness
where T is the aquifer transmissivity. Our models indicate a good linear is shown in Supplementary Figs 10–13. The simulation in Supplementary Fig. 10
relationship (R2 ¼ 0.64, Supplementary Table 1) between GRT and spring flow illustrates the differences between slope and roughness as cost layers. Roughness
recession timescales as expected. However, we note that GRTs for non-uniform provides a more uniformly variable cost layer that is impacted less by the macro-
flow fields increase for convergent and decrease for divergent flow geometries57. scale geomorphology of the rift and more by local ground conditions that would be
Thus our results, which assume uniform flow, are conservative with respect to the experienced by an agent moving over the terrain. The conclusions reported here
modelled timescales of spring persistence. use surface roughness as the primary cost layer. Slope restricts mobility to a greater
extent, providing a more conservative set of results with respect to agent mobility.
It is important to note that by using roughness our results favour mobility over
Monte Carlo experiments and groundwater model sensitivity. Sensitivity of isolation. The overall geographical pattern of results reported here, however, do not
modelled spring persistence to parameter uncertainty was tested by varying each change in broad terms (Supplementary Figs 6,7,11) if slope is used, although the
parameter by ±25% in turn and interrogating the model output (Fig. 4). Monte number of successful journeys does decline and the cross-rift linkages in particular
Carlo experiments (MCEs) were prohibitively computationally expensive to run on around the Ethiopian Highlands cease under the present day dry scenario
the full models. However, we felt it important to explore the full range of parameter (Supplementary Fig. 7). Similarly by varying the scaling of cost used impacts on
uncertainty. Thus, we used the result that spring longevity was well correlated with agent mobility (Supplementary Fig. 12); the scaling that most approximates
GRT, and ran the MCE on GRT for each spring for 10,000 combinations of Naismith’s Rule was used (Supplementary Fig. 9). The contribution of each
parameters sampled randomly from the range defined as follows. The range of hydrological component was also modelled separately and these results are shown
parameter uncertainty was assumed to be well constrained for mapped parameters in Supplementary Figs 13 and 14. At each tick or time-step the model hatches for
(catchment length, area and slope) and allowed to vary through a standard 4,000 agents who then attempt a journey. The model was routinely stopped after
deviation equal to 10% of the mean. Hydraulic properties k0 and ne were free to 100 ticks or after 400,000 attempted journeys. The results are not sensitive to the
vary through one standard deviation as defined by Gleeson et al.41 The saturated run-time beyond 100 ticks±20, but connectivity can be reduced if the model is
flow thickness was deemed to be more uncertain and free to vary through a normal stopped prematurely, although most connections are established by 50 ticks. In
distribution with a standard deviation of 25% of the mean (Supplementary Fig. 3). summary, the model parameters include: (1) maximum travel distance in 3 days
without water; (2) cost-scaling of the landscape reducing travel distance as the
terrain become rougher; (3) the maximum turn angle set for an agent which
Agent-based modelling. ABM was implemented using the Pathway Analysis
determines the degree of lateral exploration versus forward motion; and (4) the
Through Habitat (PATH) algorithm run in Netlogo58,59. A copy of the original
model run-time which can impact on the probability of a connection being found
model code is available from: http://extras.springer.com/2012/978-1-4614-1256-4
by an agent between two resource points. The conclusions presented here are
and a copy of the modified code is included in Supplementary Information. The
robust across a range of parameter values.
algorithm involves launching agents from defined habitat patches, in this case
Outputs from the PATH model were analysed in three ways: (1) visually maps
water, to simulate the journey made by individuals through a landscape until they
of successful agent journey are exported from Netlogo as raster files and uploaded
either arrive at another suitable habitat patch or die. There are four model
for visualization and analysis in ArcMap; (2) within ArcMap the PATH output can
requirements: (1) a value for total travel distance in a unit of time; (2) a map of
be converted into a binary land cover map (crossed or not crossed) for analysis in
water-patches; (3) a map showing the energetic cost associated with the terrain
Fragstats68 which produces a range of connectivity variables for habitat patches;
travelled; and (4) a map showing the potential for death (lethality) associated with
and (3) in addition the PATH model provides a matrix of successful journeys
the terrain travelled. While the original PATH model has a capacity to include a
between named habitat (water) patches which can be exported and plotted within
spatially variable measure of lethality, in the interests of parsimony, we have chosen
Matlab (Fig. 6). Supplementary Tables 2–6 provides output metrics for the main
to set this to zero. This gives the agents the maximum chance of reaching their
model runs and associated sensitivity analysis. Supplementary Fig. 15 reports the
destination and achieving dispersal and thereby tests rigorously the potential for
results of a Principle Components Analysis on the all the metrics derived from the
population isolation. Furthermore, sensitivity tests using a range of possible
different model runs across the four climate scenarios. The 95% confidence ellipses
lethality do not materially alter the potential for dispersal reported here. Modelling
clearly show the statistically significant differences between the four climate
was to a resolution of 1 km, with relief roughness60 used as the cost layer, scaled to
scenarios modelled. These differences are reinforced in Supplementary Table 6,
reflect the variation in walking speed with slope/roughness as defined by
which reports the mean model runs and the 95% error margins.
Naismith’s Rule61 (Supplementary Fig. 9a). A maximum travel time of 3 days
without water was used62. Daily walking distances were based on modern human
walking speeds. Derived from on a large sample (N ¼ 3,500) Tarawneh63 quotes
Data availability. Locations of springs, Matlab model code for estimating spring
pedestrian walking speeds for different age categories of the order of o20 y,
persistence, and the NetLogo ABM are available at: https://doi.org/10.6084/
1.29 m s  1; 21–30 y, 1.49 m s  1; 31–41 y, 1.47 m s  1; and 465 y, 1.1 m s  1.
m9.figshare.c.3721141 (ref. 69)- GIS shapefiles for the digitized hydrological
While modern humans are capable of walking at speeds upwards of 2.5 m s  1,
features used in the analysis are available from the authors on request.
especially for short distances, they typically choose to walk at their ‘preferred
walking speed’ which is variously defined as 1.3 and 1.4 m s  1 (refs 64–67).
We acknowledge that there may be other aspects of terrain, which might impede
movement such as vegetation type and surface albedo or ‘going characteristics’,
References
but on the basis of the range of data, one can assume that daily walking distance for 1. Behrensmeyer, A. K. Climate change and human evolution. Science 311,
a modern human is likely to vary between 40 and 50 km over a 10 h period. Taking 476–478 (2006).
the least conservative value, the model scenarios were run with a 3 days distance of 2. Cerling, T. E. et al. Woody cover and hominin environments in the past 6
150 km representing a maximum possible travel distance. Values are likely to have million years. Nature 476, 51–56 (2011).
been less for pre-Homo agents and a hominin troop may have moved at the pace of 3. deMenocal, P. B. Climate and human evolution. Science 331, 540 (2011).
the slowest member. All conclusions were tested specifically against travel distances 4. Potts, R. Variability selection in hominid evolution. Evol. Anthropol.: Issues,
of 120, 150 and 180 km (Supplementary Table 5, Supplementary Fig. 6). While News, Rev. 7, 81–96 (1998).
lower travel distances of say 80 km or less in 3 days are conceivable, they simply 5. Maslin, M. A. et al. East African climate pulses and early human evolution.
reduce the connectivity between points (that is, springs) in the driest scenario and Quat. Sci. Rev. 101, 1–17 (2014).
delay the switching on of cross-rift routes as climate ameliorates. 6. Maslin, M. A. & Trauth, M. H. in The First Humans–Origin and Early
Implementation of travel distance is based on a simple ‘energy quota’. Each Evolution of the Genus Homo 151–158 (Springer, 2009).
starting agent has an energy value of 150. If the cost is zero (that is, the terrain flat) 7. Vrba, E. S. Environment and evolution: alternative causes of the temporal
then they can travel 150 units (or cells) which is equivalent to 150 km given that distribution of evolutionary events. S. Afr. J. Sci. 81, 229–236 (1985).

NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696

8. Shultz, S. & Maslin, M. Early human speciation, brain expansion and dispersal 43. Taylor, R. G. et al. Evidence of the dependence of groundwater resources on
influenced by African climate pulses. PLoS ONE 8, e76750 (2013). extreme rainfall in East Africa. Nat. Clim. Change 3, 374–378 (2013).
9. Trauth, M. H., Maslin, M. A., Deino, A. & Strecker, M. R. Late Cenozoic 44. Cuthbert, M. O. & Tindimugaya, C. The importance of preferential flow in
moisture history of East Africa. Science 309, 2051–2053 (2005). controlling groundwater recharge in tropical Africa and implications for
10. Marean, C. W. et al. A new research strategy for integrating studies of modelling the impact of climate change on groundwater resources. J. Water
paleoclimate, paleoenvironment, and paleoanthropology. Evol. Anthropol.: Clim. Change 1, 234–245 (2010).
Issues, News, Rev. 24, 62–72 (2015). 45. Bonsor, H. C., MacDonald, A. M. & Davies, J. Evidence for extreme variations
11. Potts, R. Environmental and behavioral evidence pertaining to the evolution of in the permeability of laterite from a detailed analysis of well behaviour in
Early Homo. Curr. Anthropol. 53, S299–S317 (2012). Nigeria. Hydrol. Process. 28, 3563–3573 (2014).
12. Wood, B. Human evolution: fifty years after homo habilis. Nature 508, 31–33 46. Bretzler, A. et al. Groundwater origin and flow dynamics in active rift
(2014). systems–A multi-isotope approach in the Main Ethiopian Rift. J. Hydrol. 402,
13. Ashley, G. M. Orbital rhythms, monsoons, and playa lake response, Olduvai 274–289 (2011).
Basin, equatorial East Africa (ca. 1.85–1.74 Ma). Geology 35, 1091–1094 (2007). 47. Kebede, S. et al. Groundwater origin and flow along selected transects in
14. Deocampo, D. M., Cuadros, J., Wing-Dudek, T., Olives, J. & Amouric, M. Ethiopian rift volcanic aquifers. Hydrogeol. J. 16, 55–73 (2008).
Saline lake diagenesis as revealed by coupled mineralogy and geochemistry of 48. Rango, T. et al. The dynamics of central Main Ethiopian Rift waters:
multiple ultrafine clay phases: Pliocene Olduvai Gorge, Tanzania. Am. J. Sci. evidence from dD, d 18 O and 87 Sr/86Sr ratios. Appl. Geochem. 25, 1860–1871
309, 834–868 (2009). (2010).
15. Döll, P. & Fiedler, K. Global-scale modeling of groundwater recharge. Hydrol. 49. Swanson, S. K. & Bahr, J. M. Analytical and numerical models to explain steady
Earth Syst. Sci. 12, 863–885 (2008). rates of spring flow. Ground Water 42, 747–759 (2004).
16. Joordens, J. C. et al. An astronomically-tuned climate framework for hominins 50. Brutsaert, W. Hydrology: an Introduction (Cambridge University Press, 2005).
in the Turkana Basin. Earth Planet. Sci. Lett. 307, 1–8 (2011). 51. Verhoest, N. E. & Troch, P. A. Some analytical solutions of the linearized
17. Magill, C. R., Ashley, G. M., Domı́nguez-Rodrigo, M. & Freeman, K. H. Dietary Boussinesq equation with recharge for a sloping aquifer. Water Resour. Res. 36,
options and behavior suggested by plant biomarker evidence in an early human 793–800 (2000).
habitat. Proc. Natl Acad. Sci. 113, 2874–2879. 201507055 (2016). 52. Paul, J. D., Roberts, G. G. & White, N. The African landscape through space
18. Magill, C. R., Ashley, G. M. & Freeman, K. H. Water, plants, and early human and time. Tectonics 33, 898–935 (2014).
habitats in eastern Africa. Proc. Natl Acad. Sci. 110, 1175–1180 (2013). 53. Bailey, G. N., Reynolds, S. C. & King, G. C. P. Landscapes of human evolution:
19. Lloyd, J. A review of aridity and groundwater. Hydrol. Process 1, 63–78 (1986). models and methods of tectonic geomorphology and the reconstruction of
20. Scanlon, B. R. et al. Global synthesis of groundwater recharge in semiarid and hominin landscapes. J. Hum. Evol. 60, 257–280 (2011).
arid regions. Hydrol. Process 20, 3335–3370 (2006). 54. Ertuna, C. Technical note: flow regime of the Mzima springs in Kenya. Proc.
21. Simmers, I. Understanding Water in a Dry Environment: IAH International Inst. Civil Eng. 67, 833–840 (1979).
Contributions to Hydrogeology 23 (Taylor & Francis, 2003). 55. Xu, Y. & Beekman, H. E. (eds). Groundwater Recharge Estimation in Southern
22. Cuthbert, M. O. & Ashley, G. M. A spring forward for hominin evolution in Africa (UNESCO, 2003).
East Africa. PLoS ONE 9, e107358 (2014). 56. Voss, C. I. & Soliman, S. M. The transboundary non-renewable Nubian Aquifer
23. Trauth, M. H. et al. Human evolution in a variable environment: the amplifier System of Chad, Egypt, Libya and Sudan: classical groundwater questions and
lakes of Eastern Africa. Quat. Sci. Rev. 29, 2981–2988 (2010). parsimonious hydrogeologic analysis and modeling. Hydrogeol. J. 22, 441–468
24. Bonabeau, E. Agent-based modeling: methods and techniques for simulating (2014).
human systems. Proc. Natl Acad. Sci. 99, 7280–7287 (2002). 57. Vasseur, G., Rousseau-Gueutin, P. & de Marsily, G. Time constant of hydraulic-
25. Strait, D. S. & Wood, B. A. Early hominid biogeography. Proc. Natl Acad. Sci. head response in aquifers subjected to sudden recharge change: application to
96, 9196–9200 (1999). large basins. Hydrogeol J. 23, 915–934 (2015).
26. Brunet, M. et al. A new hominid from the Upper Miocene of Chad, Central 58. Hargrove, W. W. & Westervelt, J. D. in Ecologist-Developed Spatially-Explicit
Africa. Nature 418, 145–151 (2002). Dynamic Landscape Models 211–222 (Springer, 2012).
27. Lorenzen, E., Heller, R. & Siegismund, H. R. Comparative phylogeography of 59. NetLogo (Wilensky, U., Center for Connected Learning and Computer-Based
African savannah ungulates. Mol. Ecol. 21, 3656–3670 (2012). Modeling, Northwestern University, Evanston, IL, 1999).
28. Zinner, D., Wertheimer, J., Liedigk, R., Groeneveld, L. F. & Roos, C. Baboon 60. Jenness, J. S. Calculating landscape surface area from digital elevation models.
phylogeny as inferred from complete mitochondrial genomes. Am. J. Phys. Wildl. Soc. Bull. 32, 829–839 (2004).
Anthropol. 150, 133–140 (2013). 61. Rees, W. Least-cost paths in mountainous terrain. Comput. Geosci. 30, 203–209
29. Calow, R. C. et al. Ground water security and drought in Africa: linking (2004).
availability, access, and demand. Ground Water 48, 246–256 (2010). 62. Popkin, B. M., D’Anci, K. E. & Rosenberg, I. H. Water, hydration, and health.
30. Tucker, J., MacDonald, A., Coulter, L. & Calow, R. C. Household water use, Nutr. Rev. 68, 439–458 (2010).
poverty and seasonality: wealth effects, labour constraints, and minimal 63. Tarawneh, M. S. Evaluation of pedestrian speed in Jordan with investigation of
consumption in Ethiopia. Water Resour. Rural Dev. 3, 27–47 (2014). some contributing factors. J. Saf. Res. 32, 229–236 (2001).
31. Beven, K. J. Rainfall-runoff Modelling: the Primer (John Wiley & Sons, 2011). 64. Browning, R. C., Baker, E. A., Herron, J. A. & Kram, R. Effects of obesity and
32. Wood, B. & K Boyle, E. Hominin taxic diversity: fact or fantasy? Am. J. Phys. sex on the energetic cost and preferred speed of walking. J. Appl. Physiol. 100,
Anthropol. 159, S37–S78 (2016). 390–398 (2006).
33. Hammer, M. F., Woerner, A. E., Mendez, F. L., Watkins, J. C. & Wall, J. D. 65. Levine, R. V. & Norenzayan, A. The pace of life in 31 countries. J. Cross-cultural
Genetic evidence for archaic admixture in Africa. Proc. Natl Acad. Sci. 108, Psychol. 30, 178–205 (1999).
15123–15128 (2011). 66. Minetti, A. E. Biomechanics: walking on other planets. Nature 409, 467–469
34. Prăvălie, R. Drylands extent and environmental issues. A global approach. (2001).
Earth-Sci. Rev. 161, 259–278 (2016). 67. Mohler, B. J., Thompson, W. B., Creem-Regehr, S. H., Pick, Jr H. L.
35. Taylor, R. G. et al. Ground water and climate change. Nat. Clim. Change 3, & Warren, Jr W. H. Visual flow influences gait transition speed and preferred
322–329 (2013). walking speed. Exp. Brain Res. 181, 221–228 (2007).
36. Holman, I. P., Allen, D. M., Cuthbert, M. O. & Goderniaux, P. Towards best 68. McGarigal, K., Cushman, S. A. & Ene, E. FRAGSTATS v4: Spatial Pattern
practice for assessing the impacts of climate change on groundwater. Hydrogeol. Analysis Program for Categorical and Continuous Maps. Computer software
J. 20, 1–4 (2012). program produced by the authors at the University of Massachusetts, Amherst.
37. GRDC, Global Runoff Data Centre in the Bundesanstalt fuer Gewaesserkunde. Available at: http://www.umass.edu/landeco/research/fragstats/fragstats.html
http://www.bafg.de/GRDC (56068 Koblenz, Germany, 2007). (2012).
38. Hayes, G. et al. Seismicity of the Earth 1900–2013 East African Rift. Report No. 69. Cuthbert, M. O. Data and Models for Groundwater Hydro-Refugia & Hominin
2331–1258 (US Geological Survey, 2014). Evolution in East Africa. Figshare https://doi.org/10.6084/m9.figshare.c.
39. Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from 3721141 (2017).
spaceborne elevation data. EOS, Trans. Am. Geophys. Union 89, 93–94 (2008).
40. Jarvis, A., Reuter, H., Nelson, A. & Guevara, E. Hole-filled Seamless SRTM Data
V4 (International Centre for tropical Agriculture (CIAT), 2008).
41. Gleeson, T., Moosdorf, N., Hartmann, J. & Beek, L. A glimpse beneath earth’s Acknowledgements
surface: GLobal HYdrogeology MaPS (GLHYMPS) of permeability and M.O.C. was financially supported by the European Community’s Seventh Framework
porosity. Geophys. Res. Lett. 41, 3891–3898 (2014). Programme (FP7/2007–2013) under the grant agreement no. 299091; T.G. was supported
42. MacDonald, A. M., Bonsor, H. C., Ó Dochartaigh, B. É. & Taylor, R. G. by the NSERC funding; we are grateful to Liz Flanary for GIS support and Weronika
Quantitative maps of groundwater resources in Africa. Environ. Res. Lett. 7, Wodzinska for providing hominin site location data. M.R.B. was supported by the
024009 (2012). UK Natural Environment Research Council (NE/H004211/1).

10 NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15696 ARTICLE

Author contributions Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in
M.O.C./T.G.: conceived the paper and designed the hydrological approach; C.J.M./M.O.C./T.G.: published maps and institutional affiliations.
spring mapping and linked GIS analysis; M.O.C./C.J.M.: hydrological modelling; S.C.R./M.R.B./
A.C.N.: design and implementation of the agent-based modelling and associated GIS prepara-
Open Access This article is licensed under a Creative Commons
tions; M.O.C.: drafted the paper with all authors providing data interpretation and comments.
Attribution 4.0 International License, which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give
Additional information appropriate credit to the original author(s) and the source, provide a link to the Creative
Supplementary Information accompanies this paper at http://www.nature.com/ Commons license, and indicate if changes were made. The images or other third party
naturecommunications material in this article are included in the article’s Creative Commons license, unless
indicated otherwise in a credit line to the material. If material is not included in the
Competing interests: The authors declare no competing financial interests. article’s Creative Commons license and your intended use is not permitted by statutory
Reprints and permission information is available online at http://npg.nature.com/ regulation or exceeds the permitted use, you will need to obtain permission directly from
reprintsandpermissions/ the copyright holder. To view a copy of this license, visit http://creativecommons.org/
licenses/by/4.0/
How to cite this article: Cuthbert, M. O. et al. Modelling the role of groundwater
hydro-refugia in East African hominin evolution and dispersal. Nat. Commun. 8, 15696
doi: 10.1038/ncomms15696 (2017). r The Author(s) 2017

NATURE COMMUNICATIONS | 8:15696 | DOI: 10.1038/ncomms15696 | www.nature.com/naturecommunications 11