Did a 3800-year-old Mw ~9.5 earthquake trigger major social disruption in the Atacama Desert?

Early inhabitants along the hyperarid coastal Atacama Desert in northern Chile developed resilience strategies over 12,000 years, allowing these communities to effectively adapt to this extreme environment, including the impact of giant earthquakes and tsunamis. Here, we provide geoarchaeological evidence revealing a major tsunamigenic earthquake that severely affected prehistoric hunter-gatherer-fisher communities ~3800 years ago, causing an exceptional social disruption reflected in contemporary changes in archaeological sites and triggering resilient strategies along these coasts. Together with tsunami modeling results, we suggest that this event resulted from a ~1000-km-long megathrust rupture along the subduction contact of the Nazca and South American plates, highlighting the possibility of Mw ~9.5 tsunamigenic earthquakes in northern Chile, one of the major seismic gaps of the planet. This emphasizes the necessity to account for long temporal scales to better understand the variability, social effects, and human responses favoring resilience to socionatural disasters.

INTRODUCTION

Earthquakes and tsunamis are among the most catastrophic events to affect human societies. They are marked and revealing moments in which human adaptation to an environment is challenged, and resilience capacities are manifested (). Coastal communities throughout the world are highly vulnerable in the face of these natural hazards, and the social disasters they produce will certainly recur in the coming decades, putting human lives, infrastructure, and socioenvironmental assets at risk (). As a result, seismic and tsunami risk assessments are becoming critical studies for the world’s growing population that is increasingly exposed to active tectonic boundaries (). Since the beginning of the present millennium, a cluster of giant earthquakes has occurred around the world with notable social consequences including the loss of more than 700,000 human lives (). These include the giant megathrust earthquakes of moment magnitude (Mw) ~9.1 Sumatra-Andaman (2004); Mw 8.8 Maule, central Chile (2010); and Mw 9.1 Tohoku-Oki (2011), which caused large tsunamis that severely affected adjacent coasts and produced transoceanic effects.

The unexpected characteristics of some of these events have significantly changed our views about their genesis and impact (). These findings have demonstrated that long temporal scales must be considered to adequately assess the magnitude, frequency, and sources of these events (–) and that understanding the different trajectories followed by human societies to face these socio-natural disasters may teach us how to address them in the future. The tectonic margin along northern Chile is a privileged location to study the occurrence and characteristics of tsunamigenic megathrust earthquakes given its high potential for these kind of events due to the adjacent subduction of the Nazca plate beneath the South American plate (Fig. 1A) (, ). However, a lack of interdisciplinary research has restricted our understanding of the occurrence, magnitude, and effects of prehistoric megathrust earthquakes and tsunamis to sparse geological records in the region (–), providing an insufficient time scale for analysis but lacking fundamental data to assess the social impacts of these events and the human responses that were triggered by them. Here, we present the first geoarchaeological evidence for a major earthquake and tsunami that struck coastal northern Chile ~3800 calendar years before present (cal yr B.P.), causing notable social disruption among ancient hunter-gatherer-fishers (HGFs) of the Atacama Desert and triggering resilience strategies. By resilience, we mean the capacity of human communities to absorb changes occurring after a socioenvironmental disturbance, allowing for their long-term adaptation (). In this sense, human societies follow different historical resilience trajectories, in contrast to the “return to the preshock state” that characterizes resilient behavior in noncultural communities ().

Fig. 1.

Tectonic plates setting, major historical tsunamigenic megathrusts, sites, and tsunami modeling results supporting a Mw ~9.5 earthquake scenario ~3800 cal yr B.P. in northern Chile.

(A) Proposed ca. 1000-km-long megathrust rupture for a ~3800 cal yr B.P. tsunamigenic earthquake, encompassing three major segments in front of Iquique, Taltal, and Atacama regions. (B) Sites of geological and archaeological observations for coastal uplift and paleotsunami deposits dated to ~3800 cal yr B.P. and tsunami modeling results. 1: Pabellón de Pica; 2: Paquica Sur; 3: Cobija; 4: Michilla; 5: Playa Grande de Mejillones; 6: Mejillones Peninsula; 7: Zapatero; 8: San Ramón-15 mine; 9: Los Bronces; 10: Hornos de Cal; 11: Playa San Pedro; 12: Bahía Cisne, Caldera. (C) Tsunami run-ups associated with the major historical tsunamigenic earthquakes of 1868, 1877, and 1922 CE in northern Chile (from National Hydrographic and Oceanographic Service). (D) Map showing tectonic plates coupling indicating a high potential for large megathrust earthquakes along northern Chile (). (E) Seismotectonic setting along northern Chile, with the subduction of the Nazca plate beneath the South American plate (, ) and historical megathrusts ruptures (, ).

The coastal Atacama Desert

The coastal Atacama Desert along northern Chile is an arid to hyperarid environment (rainfall of <1 mm/year) encompassing around 1000 km between Arica and Huasco (Fig. 1A), with low terrestrial biomass production, limited availability of dependable water sources, almost barren of vegetation except for fog-driven vegetation assemblages, and virtual absence of precipitation in many areas with only sporadic heavy rainfalls driving localized alluviation mostly during El Niño–Southern Oscillation episodes (). However, the offshore environment is one of the most productive marine ecosystems in the world because of the cold Humboldt Current system and southwesterly wind-driven upwelling cells that create one of the richest pelagic and benthic fisheries in the world (). Because of the stability, richness, and predictability of marine resources, there has been continuous human occupation for over 12,000 years along this coastline, showing an efficient adaptation to maritime conditions until traditional HGF communities were absorbed into the capitalist economy over the past few centuries (, ). Faunal, technological, and isotopic data show evidence for Early to Late Holocene continuity, suggesting a conservative way of life characterized by a traditional subsistence base and an improved specialized technology, well adapted to local ecosystems (). This millennial continuity in dietary practices and technology was favored by resilience strategies that allowed these communities to effectively adapt to this extreme environment despite environmental change and the effects of episodic catastrophic events common on the Pacific coast of South America, such as debris flows and especially tsunamigenic earthquakes ().

The coastal Atacama Desert has a high potential for the occurrence of megathrust earthquakes as it is characterized by a first-order tectonic plate boundary with a subduction convergence rate of the Nazca plate beneath the South American plate estimated at 65 mm/year (). Coastal northern Chile is considered a highly coupled mature seismic gap (Fig. 1B) (, ), since no major episodes have occurred here since the 1877 CE Mw ~8.8 Iquique event, a giant tsunamigenic megathrust earthquake that ruptured at least ~500 km of the subduction margin and produced tsunami waves of 12 to 20 m in the epicentral zone (, ). Significant tsunami run-ups were recorded for this event, as well as for the historical 1922 CE Atacama (Chañaral-Huasco) and 1868 CE Arica earthquakes (Fig. 1C) (, ). Holocene paleotsunami deposits along this extremely arid and rocky coast are more difficult to identify but have been recently recognized in association with archaeological sites in the Taltal region (). In addition, submarine slump deposits and an associated erosional basal contact preserved in laminated sediments in Mejillones Bay (23°S) suggest that a predecessor of the major 1877 CE event probably occurred around 1409 to 1449 CE (1429 ± 20 CE) (). This tentatively corresponds with a paleotsunami dated from the study of coastal boulders to around 1300 to 1600 CE in the Atacama region () and closely dates to a historically documented distant tsunami that struck the Japanese coasts on 7 September 1420 CE (, ). An 863 ± 199 CE episode described from Taltal () could also correspond to the distant tsunami of 18 September 799 CE reported from historical records in Japan (). This paper deals with a probable predecessor of these more recent events, occurring at ~3800 cal yr B.P. and affecting more than 1000 km of the Atacama coast and the ancient local population of HGFs inhabiting it.

RESULTS

Regional geological evidence for Mid-Holocene uplifted littoral deposits

Geological and archaeological evidence for a ~3800 cal yr B.P. megathrust earthquake (Fig. 1, D and E) derives from uplifted littoral beach (Fig. 2A) and paleotsunami deposits (Fig. 2B), erosional contacts, and fractures in archaeological sites (Fig. 2, B and C) (see figs. S1 to S7, S9 to S11, S16, S17, S19, S20, and S24 to S28). Geological evidence from these and previous observations (, ) shows conspicuous littoral deposits located up to 6 to 7 m above the present mean sea level (a.s.l.). Littoral deposits are characterized by coarse well-rounded gravel, some angular and boulder clasts embedded within a sandy matrix with abundant carbonate shell fragments and well-preserved mollusk shells, as in the case of those previously reported from Los Bronces and Hornos de Cal sites in the Taltal area (), or by sandy facies with low-angle cross-bedding with abundant shell fragments such as Playa Grande site in Mejillones (). Deposits with similar characteristics, which have been interpreted as uplifted beach-berm or shoreface facies, have been observed over more than 600 km of coastline (Figs. 1D and 2A).

Fig. 2.

Selected sites for geological and archaeological observations of Late Holocene coastal uplift and tsunami deposits (see Fig.1B for site locations).

(A) Uplifted littoral deposits are located up to 6 m a.s.l. at Pabellón de Pica (20.89°S), overlain by distinctive deposits interpreted as tsunami layers and debris flow. Littoral facies, showing stratification and abundant shell fragments interpreted as beach-berm deposits, are overlain by distinctive deposits interpreted as tsunami layers showing significant lateral extension, fragments, and some articulated bivalves and other faunal content typical from subtidal and intertidal environments (Protothaca thaca, Choromytilus chorus, and Megabalanus sp.), extensive angular fragmentation, and incipient imbrication toward the sea [e.g., ()]. (B) Tsunami deposit located at 14 m a.s.l. embedded in and reworking Mid-Holocene archaeological layers at Paquica Sur (21.91°S). (C) Archaeological strata at San Ramón-15 mine (25.38°S), affected by a seismically induced crack overlain by colluvial wedge and soil layer.

These deposits have yielded 17 radiocarbon ages from seven sites, with minimum values of 3145 ± 373 cal yr B.P. at Pabellón de Pica (20.89°S), 3441 ± 257 cal yr B.P. at Zapatero (24.93°S), 3345 ± 403 cal yr B.P. at Hornos de Cal (25.51°S), and 4048 ± 350 cal yr B.P. at Playa San Pedro (25.51°S), which are consistent with previous calibrated age ranges from other sites in the region (Fig. 3) (, , ).

Fig. 3.

Inferred ages for a tsunamigenic megathrust episode occurred during the Late Holocene, deduced from 12 sites located along northern Chile (see location in Fig. 1).

Horizontal bars correspond to the minimum and maximum ages for each calibrated age range. PP, Pabellón de Pica; PS, Paquica Sur archaeological site; CO, Cobija; Mi, Michilla; PG, Playa Grande de Mejillones; MP, Mejillones Peninsula; Za, Zapatero archaeological site; SR, San Ramón archaeological site; LB, Los Bronces archaeological site; HC, Hornos de Cal archaeological site; SP, Playa San Pedro; BC, Bahía Cisnes.

Taking into account the maximum inferred elevation of global sea level during the Early and Mid-Holocene (), and that it was higher than contemporary levels according to regional models (), our field observations provide evidence for coastal uplift along the major northern Chile seismic gap since the Late Holocene (Figs. 1D and 2A). This Holocene coastal uplifting was the product of megathrust episodes and reflects a process noted over the late Quaternary period for the Mejillones Peninsula (, , ).

Geoarchaeological evidence for a ~3800 cal yr B.P. tsunamigenic earthquake

Identification of distinctive deposits laid down by tsunamis in arid and hyperarid siliciclastic coastal environments of the Atacama Desert is challenged by a series of factors. Among them are the scarcity of depositional systems favoring the accumulation of sediments, limited accommodation space for tsunami deposition and preservation, and the sporadic occurrence of floods and debris flows, storms, or even wind deflation that can erode previous deposits (, ). However, paleotsunami deposits have been recently described from the hyperarid coastal Atacama Desert in protected areas associated with archaeological sites from the Taltal area (). These deposits exhibit characteristics typically described for paleotsunamis from onshore environments (), such as erosive basal contacts, seaward flame structures, landward thinning, rip-up clasts, massive texture or incipient lamination, imbrication, distinctive grain size distributions with respect to the underlying and overlying deposits, and the presence of small gastropods and other faunal remains typical of shoreline environments (). Palaeotsunami deposits from Los Bronces and Hornos de Cal sites were previously reported as dating to ~4000 cal yr B.P. ().

Palaeotsunami deposits with the above-mentioned characteristics and contemporary with the uplifted littoral beaches previously described have been recorded from eight geological and archaeological sites across the study area, suggesting minimum run-ups of up to 19 to 20 m a.s.l. (Figs. 1 and 2). Maximum ages from these paleotsunami records are provided by radiocarbon data from reworked material from underlying littoral deposits, as in the cases of Pabellón de Pica, Mejillones Peninsula, Los Bronces, and Bahía Cisne, or from underlying archaeological strata, as in the case of Paquica, Cobija, Zapatero, and Hornos de Cal (Fig. 3). In the case of Zapatero, the distribution of maximum ages from paleotsunami deposits most probably reflects a combination of both the reworking of material from the preexisting Mid-Holocene shell midden and from underlying littoral deposits. As expected, the 19 maximum ages from paleotsunami deposits yielded dates older than 3800 cal yr B.P., corresponding to reworked material (table S1). Minimal ages from paleotsunami deposits are scarce (six ages), but these correspond to charcoal from small hearths or shells from archaeological layers directly overlying these deposits (Fig. 3). However, several factors challenge precise radiocarbon age determinations for geoarchaeological records along these coasts. Among them is the uncertainty related to the contemporaneity of the analyzed material in a stratigraphic layer with the event to be dated (i.e., “old wood” effect) (); the possible inclusion of reworked material in distinctive layers; and, in the case of marine organisms, the regional radiocarbon reservoir effect, at present known only for a large regional scale encompassing all of northern Chile (). Consequently, it is not possible to assert whether the geoarchaeological evidence presented here resulted from a single, unique event or from more than one event occurring at close temporal intervals. Considering the partial recovery of the land that follows large subduction earthquakes (), the possibility that this evidence occurred as the result of a single event is supported by (i) the singularity of these contemporary well-preserved uplifted littoral marine deposits, with no other similar conspicuous Mid- to Late Holocene deposits widely recognized in the region; (ii) the relative contemporaneity with tsunami deposits distributed in several locations all along the northern Chilean coast (Fig. 3 and table S1); (iii) the occurrence of a conspicuous single-event deposit corresponding to a trans-Pacific tsunami identified in New Zealand and other islands and dated between 3500 and 4500 cal yr B.P. (); and (iv) the severe social transformations that occurred simultaneously and which have not been reported for any other time period (see below).

Considering this single-event scenario, we identified the maximum age deduced from the minimum age ranges from uplifted littoral deposits (3517 cal yr B.P., from 3145 ± 373 cal yr B.P. at Pabellón de Pica), together with the minimum age taken from the maximum age ranges from hearths and shells overlying paleotsunami deposits (4150 cal yr B.P., from 4279 ± 129 cal yr B.P. from a charcoal fragment at Hornos de Cal), to obtain a chronostratigraphically consistent age range of 3834 ± 317 cal yr B.P. for the tectonic event that would have uplifted littoral deposits all along the study region, generated a paleotsunami, and triggered social disruption at a regional scale (Fig. 3 and table S1).

Archaeological evidence for a ~3800 cal yr B.P. tsunamigenic earthquake

Evidence for paleotsunami-induced erosion and sedimentation contemporary with the ~3800 cal yr B.P. event has been identified at five archaeological sites along the coast of northern Chile (Paquica, Cobija, Zapatero, Los Bronces, and Hornos de Cal). The Zapatero site (Figs. 4 and 5 and fig. S8) is a 6000-m2, 2-m-deep, shell midden (, , ). The effects of strong currents that we associate with a Late Holocene paleotsunami are evident in the erosional surfaces of the upper part of the midden at different locations on the site (Fig. 5C and figs. S9 to S11, S16, and S17) and by the partial or total destruction of two Late Holocene stone-built architectonic features identified during excavations. The stones of the walls from one of these buildings (structure 3) are toppled seaward (Fig. 6A and figs. S12 and S14), indicating that the agent responsible was the strong currents consistent with tsunami backwash (Figs. 6C and 5E). Another structure, also destroyed during prehispanic times (structure 1), preserves evidence of landward imbrication of toppled stones (fig. S13). A minimal age of 4191 ± 108 cal yr B.P. from a hearth overlaying structure 3 and a maximum age of 4689 ± 136 obtained from in situ material under the collapsed stones of structure 1 are both consistent with the wide age range determined from uplifted littoral and paleotsunami records around 3834 ± 317 cal yr B.P. (Fig. 3 and table S1).

Fig. 4.

Geomorphological context for the Zapatero archaeological site.

(A) View of the site toward the northeast. (B) View toward the south. NW, northwest; W, west.

Fig. 5.

Geoarchaeological evidence for coastal uplift and tsunami records at the Zapatero archaeological site (site 7 in Fig. 1B).

(A) Stratigraphic column and associated results from microstratigraphic analyses (whole individual radioles per sample and algae fragments per sample were counted). (B) Geomorphological profile and stratigraphic column showing site location with respect to modern sea level. (C) Distinctive tsunami deposit with gravel and broken shells overlaying supralittoral aeolian sands. (D) Geomorphological context of Zapatero showing the location of main excavations. (E) Proposed model for the ~3800 cal yr B.P. event that uplifted littoral deposits and caused reworking of archaeological strata by tsunami inundation. VE, Vertical Exaggeration.

Fig. 6.

Toppled stones from structure 3 at Zapatero archaeological site.

(A) Plan view sketch of toppled stones from structure 3 after excavation. (B) Interpretation of the impact of backwash of tsunami waves leaving toppled stones. (C) Photo of toppled stones from structure 3 during archaeological excavation.

A third structure excavated at Zapatero (structure 2) does not present evidence of collapsed stones but shows a post-abandonment layer both inside and outside the structure, interpreted as a tsunami deposit (), because of, to the lower irregular and erosive contact with underlying layers, a small channel with an erosional basal contact filled with sand and scarce archaeological material and the presence of rip-up clasts (figs. S16 and S17). In areas of the Zapatero shell midden where no stone structures were built by prehispanic populations, we also identified tsunami deposits. This interpretation is based on field observations (lack of internal structure, randomly orientated artefactual and ecofactual remains, and abrupt erosional basal contacts) from at least four profiles located at different areas of the site (fig. S9 to S11C2 and S29 for location of profiles 1 to 7). A microstratigraphic analyses performed on one of these profiles showed a high frequency of marine organisms such as whole echinoderm spines and algal fragments in the inferred tsunami deposit, which are either absent or scarce in the lower sections of the midden (Fig. 5C and figs. S30 to S36) (). The textural and morphological characteristics of the calcareous algae suggested that it was most likely red algae (Rhodophyta), a marine organism that live in deep waters, sandy beaches, or on rocky substrate, indicating a marine origin for the sediments. Furthermore, a 30× to 60× increase in the concentration of echinoderm radioles in these layers further supports the marine origin (table S2). Microstratigraphic analyses also showed an increased amount of terrigenous clay mixed with nodules of secondary phosphates in the tsunami deposit and increase in the grain size of moderate sorted rock fragments, as opposed to the coarse sand and well-sorted rock fragments in underlying units (table S3).

Four other archaeological sites show similar and contemporary evidence to Zapatero. Hornos de Cal, at 10 m a.s.l., is located ~90 km south of Zapatero (); San Ramón 15 (SR15), 3 km from the current coast and at 125 m a.s.l., is located ~55 km south of Zapatero (); while Cobija and Paquica, at 19 and 15 m a.s.l., respectively, are located ~260 to 330 km to the north (Fig. 2 and fig. S39). The SR15 site is a prehispanic iron oxide mine with evidence for a near-vertical crack with an average lateral separation of 3 cm in the prehispanic infill of the mining trench. Two hearths located stratigraphically near the upper end of the crack were radiocarbon-dated to 4161 ± 37 cal yr B.P., indicating that the crack was formed immediately after these ages (). A radiocarbon result from a shell fragment obtained from the top of the crack provided a maximum age of 4087 cal yr B.P. (from an age of 3812 ± 276 cal yr B.P.), while a radiocarbon result from a shell fragment in a layer directly overlying the crack provided a minimum age of 3630 cal yr B.P. (from an age of 3914 ± 284 cal yr B.P.), yielding an age range of 3859 ± 229 cal yr B.P. (Fig. 3 and table S1). This is consistent with the estimated age for the tectonic event that regionally uplifted littoral deposits and produced minimum tsunami run-ups of 19 to 20 m a.s.l. at Zapatero, Hornos de Cal, Cobija, Paquica, and other sites (Fig. 1D).

Hornos de Cal (), Cobija, and Paquica show evidence of reworked material and erosional basal contacts overlaying archaeological strata. At Zapatero, Hornos de Cal, and Paquica, shells found within the tsunami deposits were radiocarbon-dated to ages older than the underlying in situ deposits, showing evidence for the reworking of the earlier middens. These latter sites correspond to Mid-Holocene shell middens documenting residential and semi-sedentary occupations, with evidence of paleotsunami deposits indicative of the region-wide ~3800 cal yr B.P. event.

DISCUSSION

Northern Chile superrupture dynamics hypothesis

Seismic ruptures associated with Mw ≥ 9.0 earthquakes have large asperities within a segment area in the order of 0.3 (asperity area/segment area); together with a high stress drop, this would imply that seismic barriers or structural boundaries would not be effective (). If a coherent seismic rupture nucleates in the region of a large slip-weakening distance (i.e., the region with the greatest seismic slip), then the possibility that it will grow into a larger size event increases ().

The ~1000-km-long supersegment that ruptured during the 1960 Mw 9.5 Valdivia earthquake would have resulted from four to five seismic patches of higher slip (). The subduction contact along northern Chile encompassing ~1000 km between Arica and Atacama can be divided into five highly coupled patches capable of generating large megathrust earthquakes (Mw > 8.0; Fig. 1B) (). The highly coupled segments (, ) have asperity lengths that are twice those of the intersegment patches, which are necessary to synchronize for generating an Mw ≥ 9.0 superearthquake (). The northern edge of this superrupture would be limited by the Arica bend, a regional geometric feature that marks the northwest shift of the Andean forearc and where seismic and geodetic research has indicated that this part of the subduction megathrust would be a zone with dominant aseismic slip (). The southern edge of the supersegment would encompass the area up to the north of La Serena (Fig. 1), a region defined as a major decoupling area, where the Challenger Fracture Zone enters into the subduction contact ().

We propose that the ~3800 cal yr B.P. event identified here ruptured three giant segments (Fig. 1E), encompassing the northernmost Chile (e.g., Iquique, 1877 CE), the southernmost Atacama region (e.g., 1922 CE), and the intermediate Taltal segment (no historical episodes), resulting in an Mw ~9.5 earthquake as suggested by our tsunami modeling based on field observations (Fig. 1, D and E). Numerical modeling for this tsunamigenic megathrust event considered scenarios for Mw 9.2 and Mw 9.5 earthquakes with results showing reasonable agreement with inferred run-up heights for both (Fig. 1D and fig. S44). However, evidence from modern tsunami deposits indicates that a recognizable coarse sediment layer may reach only 56 to 90% of the total inundation distance depending on local relief (). Accordingly, estimated minimum inundation distances were extended by assuming that deposits represented a conservative 75% of the total, resulting in a closer match with Mw 9.5 scenarios (Fig. 1D). A Mw 9.5 scenario matches well with all our observations from sites located between Pabellón de Pica (20.89°S) and Caldera (27.25°S), as well as with numerical models for tsunami inundation in New Zealand based on an Mw 9.5 earthquake in the northern Chile-Peru region ().

Therefore, according to our archaeological and geological evidence along coastal northern Chile, together with the results of our tsunami modeling, we propose that three large segments (Iquique, Taltal, and Atacama segments) formed a supersegment, referred to as the “northern Chile superseismic gap” (see Fig. 1E), which released the greater part of the total stored elastic strain into one single seismic superrupture triggering an Mw ~9.5 earthquake (model 4; fig. S44) ~3800 years ago. The latter caused the total activation of the large segments present along coastal northern Chile, provided that the return period of these segments was synchronized, leading them to an imminent stage of multiseismic rupture.

Major social disruption in HGF communities

The occurrence of the ~3800 cal yr B.P. tsunamigenic earthquake in coastal northern Chile is contemporary with major changes in the settlement system of local hunter-gatherer communities, which had inhabited this coastline for nearly 10 millennia before this event. Although there are limitations in some of the available datasets, changes are consistent at different spatial scales, from individual sites to the entire hyperarid coast of the Antofagasta and Tarapacá regions, considering independent lines of evidence.

At a spatial mesoscale (102 to 104 km2), significant changes can be identified in the frequency and characteristic of archaeological sites corresponding to Archaic IV (~5700 to 4000 cal B.P.), Archaic V (~4000 to 2700 cal B.P.), and Archaic VI (~2700 to 1300 cal B.P.) periods (Fig. 7 and table S4). In the Taltal-Paposo area, where the most systematic archaeological surveys have been carried out, covering more than 100 km of coastline, there is a ~65% decrease in the frequency of sites after 3800 cal yr B.P. (Fig. 7B and table S5). All Archaic IV sites mentioned in the previous section were almost entirely abandoned following the tectonic event and corresponding paleotsunami inundation. We have only found archaeological evidence for ephemeral occupation consisting of small hearths and a lower frequency of artefactual and ecofactual remains overlying the fallen stone walls and inundation surfaces at Zapatero (Fig. 7) and the paleotsunami deposits at Hornos de Cal. Sedimentological evidence from the deposits overlying the crack at SR15 suggests extensive mining ceased, and the site was finally abandoned around 3900 cal yr B.P. (Fig. 2C) (). Other sites—such as Paso Malo Arcaico in the Taltal-Paposo area, Cobija 1N, Mamilla 7, and Punta Guanillos in the Tocopilla area (400 km north of Taltal)—show stratigraphic evidence for aeolian and/or colluvial deposition after 3800 cal yr B.P., consistent with the after effects of a large earthquake and with scarce archaeological materials (figs. S40 to S42).

Fig. 7.

Changes in the archaeological record after the 3800 cal yr B.P. event at regional and local scales.

(A) Artefactual and ecofactual density from occupations at the Zapatero site. (B) Artefactual and ecofactual density in residential sites of the Taltal-Paposo area. (C) Total site counts for residential occupations and cemeteries from the Taltal-Paposo area. (D) SPDs of available radiocarbon ages for coastal Antofagasta and Tarapacá regions. (E to G) Change in location of cemeteries between the Archaic IV and Archaic VI periods along the coast of the Antofagasta region.

A similar pattern emerges on a spatial macroscale along the hyperarid Atacama coast (104 to 107 km2), since Archaic V sites are also fewer in number, while Archaic VI sites have been reported as showing an increase in frequency along the coast of Tarapacá and Antofagasta regions (table S6) (, ). Unlike preceding and subsequent occupations, no evidence has been reported throughout the region for agglutinated cemeteries between 3800 and 2500 cal yr B.P. (only few sites show evidence for individual inhumation in rockshelters). Of the 19 Archaic IV sites with architecture in the region (most of them with burials), thought to be residential bases (, , ) and/or social aggregation sites () (fig. S39), none of them seems to show occupations after 3800 cal yr B.P. The youngest ages available for this kind of site come from Caleta Huelén 42 (3780 ± 90 B.P., wood) and Punta Guasilla 1 (3490 ± 290, charcoal). Other Archaic IV sites from the Antofagasta region with evidence of habitational occupation, but without architecture (i.e., Mamilla 7, Cobija 1 N), also show abandonment of the shell middens around 4000 to 3800 cal yr B.P.. In these cases, the sites were not reoccupied until several millennia later (figs. S40 and S42). Also, at a spatial macroscale, a summed probability density analysis for available radiocarbon ages in the region show an evident drop after 3800 cal yr B.P. (Fig. 7A).

Sites throughout the region are not only fewer in number during Archaic V (i.e., after the ~3800 cal yr B.P. event) but also smaller in terms of size and artefactual density. In all the sites in the Taltal-Paposo area showing open-air domestic occupations during Archaic V period (Zapatero, Morro Colorado, and Punta Morada), the horizontal extension of these deposits is smaller than those from the Archaic III and Archaic IV periods. The few Archaic V sites reported in previously uninhabited areas (), as well as the archaeological evidence overlying the paleotsunami deposits in previously occupied sites, also show occupations with comparatively lower frequencies of artefactual and ecofactual remains, indicative of ephemeral and infrequent residence by small groups. Zapatero is an ideal site to assess these differences at the scale of individual sites because it shows little evidence of looting. It is probably the single most systematically studied site of the Antofagasta coast and one of the few sites with comparable excavated areas corresponding to the Archaic IV, Archaic V, and Archaic VI periods. Extensive excavations at different locations on the site have shown that Archaic V deposits were the most poorly represented in terms of horizontal extension. Moreover, Archaic V materials were identified in yellow sand deposits with scattered charcoals and fine lenses of small shell debris (fig. S37), showing low frequency and low diversity in ecofactual categories when compared to material from underlying Archaic III-IV and overlying Archaic VI deposits (Fig. 7D and tables S7 to S9). The faunal assemblage of the Archaic III to IV period shows the highest abundance, and richness (number of species), together with a relatively even distribution of the most predominant taxa. Archaic V, on the other hand, shows a marked decrease in terms of abundance, richness, and diversity and declining evenness values for this assemblage. In the Archaic VI period, all parameters increase in respect to the previous period, with the highest diversity and evenness values of all of the samples studied (fig. S38). Kruskal-Wallis tests showed that different comparative density parameters in faunal and lithic assemblages show statistically significant differences between Archaic III to IV and Archaic V and between Archaic V and Archaic VI periods, both in terms of frequency and the density of remains (Fig. 7D and table S7 to S10). The only exception was the values for tetrapods, which show no statistical difference between periods. This result is expected, given that tetrapods present low frequencies throughout the prehistoric sequence and have been considered a marginal foodstuff in local prehistoric diets (, , ).

Smaller human occupations in the Archaic V period in relation to the Archaic IV and Archaic VI periods can also be inferred when comparing between different residential sites in terms of the frequency and density of cultural materials. Considering data produced in a previous work (), we compared densities of artifacts and ecofacts between residential sites corresponding to Archaic IV (sites 228 to 230 and Punta Morada), Archaic V (site 183), and Archaic VI (sites Plaza de Indios Norte and Paso Malo Alfarero) periods. Results are consistent with our data from Zapatero, inasmuch as the Archaic V residential site shows the lowest frequencies and densities in virtually all cultural and ecofactual categories considered (Fig. 7C and table S11) ().

The drop in summed probability distribution (SPD) of radiocarbon ages, the regional and local pattern of fewer, smaller and more ephemeral residential sites, and the lack of agglutinated cemeteries during Archaic V period all suggest a population decline following the ~3800 cal yr B.P. event. Population decline could be an effect of high death rates as a result of the tsunami, coastal abandonment by HGF communities after the event, or a combination of both. Although it is difficult to demonstrate either hypothesis with the data currently available for coastal Atacama, high death rates due to tsunamis have been documented both historically and archaeologically (). It is estimated that the 1896 AD Sanriku tsunami killed up to 83% of the population of some villages (). Several historically and archaeologically documented cases in Japan, New Zealand, North America, and the Pacific also attest that coastal areas were depopulated following massive tsunamis as the population moved inland and uphill (, , ). In the area potentially struck by the ~3800 cal yr B.P. event reported here, an agglutinated cemetery of coastal people appears after 3700 B.P. in the Tiliviche oasis, 40 km from the coast (), corresponding to a human group with a maritime economy but living far from the coastline. Coastal communities could have also chosen to move even farther from the coast and settle more permanently, for example, in the lower part of the Tarapacá quebrada (). Alternatively, some groups may have moved southward along the coast, since, in the La Serena area, a demographic increase has been suggested for the Punta Teatinos phase (after 4000 cal yr B.P.) ().

Future work is required to assess the hypothesis regarding high death rates and/or coastal abandonment, which may have followed the 3800 cal yr B.P. event. Available data show that HGF groups that chose to remain in the coastal area struck by the tsunami and occupied the smaller Archaic V sites show no evident change in technology, except for minor morphological changes in certain lithic tools or the replacement of specific raw materials (, ). Isotopic data also show evidence for continuity in local diet from before to after 3800 cal yr B.P. (), while faunal data show that the same species were exploited over time, albeit with a decrease in the diversity of the assemblages following the event. Nearly one millennium after the 3800 cal yr B.P. event, an increase in the frequency of residential sites has been reported for the Taltal-Paposo area (), while a reappearance of agglutinated cemeteries has been detected throughout the entire region (). However, after 3500 cal yr B.P., most residential sites were located farther from the shoreline and at higher altitudes than before 3800 cal yr B.P. (). Furthermore, we recorded the altitude and distance from the shore for a random sample of ~30 agglutinated cemeteries dating from 2500 to 1000 cal yr B.P. Results show a statistical difference when compared to pre–3800 cal yr B.P. burial sites (Fig. 7, E to G, and table S12), demonstrating that during the millennia that followed the event, both residential sites and agglutinated cemeteries were intentionally located at higher altitudes, usually over 20 m a.s.l., close to the minimum modeled run-up limit of the ~3800 cal yr B.P. paleotsunami.

The abandonment of previously occupied areas and changes in the mobility patterns and spatial arrangements of settlements and cemeteries were probably resilience strategies developed by hunter-gatherer societies to face the challenges of the ~3800 cal yr B.P. event and seem to have been effective for local and regional long-term adaptation. However, knowledge of these giant events and their consequences seems to wane over the passage of time, a common occurrence throughout the Pacific region (, ). By about 1000 cal yr B.P., residential sites along the coast of northern Chile were again located near the shoreline, and various agglutinated cemeteries were placed at altitudes of ~10 m a.s.l. Some of the residential sites abandoned after 3800 cal yr B.P. were again occupied as residential bases after 1500 cal yr B.P..

Implications for regional and Pacific Ocean basin-scale hazard assessment

Over millennia, repeated giant subduction zone events have provided revealing details of human adaptation to active seismic environments, and the challenges posed to the resilience capacities of social groups. Our results show that the largest tsunamigenic megathrust earthquake yet reported for northern Chile occurred ~3800 cal yr B.P. and is contemporaneous with a trans–Pacific Ocean paleotsunami that affected central Chile (), southern Chile (), and the coasts of Australia, New Zealand, Vanuatu, and possibly Japan and Russia ().

This event likely occurred as a product of spatiotemporal synchronization of large segments (), with seismic asperities large enough to produce a high stress drop outpacing barriers or structural boundaries (), which can result in superearthquakes at centennial to millennial time scales (). In subduction zones, the global occurrence of magnitude 9.0 earthquakes would be on the order of one to three per century () or even five (), but in northern Chile, a giant earthquake similar to the ~3800 cal yr B.P. event would be once every 250 to 10,000 years (, ). The latter implies a broad time span, where any of the subduction systems, not only northern Chile, may be capable of generating earthquakes much larger than known or expected today (), which causes large uncertainties in the seismic hazard assessments for a worst-case scenario for most of the subduction zones around the world.

Human communities have always developed and transmitted local knowledge about their environment, allowing them to adapt in situations of disaster risk (). Therefore, the major social disruption caused by the ~3800 cal yr B.P. event is unexpected, given that communities of HGF had inhabited the coastal Atacama Desert since at least 12,000 cal B.P. and probably faced previous tsunamigenic earthquakes. Although the specific conditions under which this event affected prehispanic communities cannot be inferred from available geoarchaeological records, one possible explanation could be that resilient responses based on local knowledge had not anticipated the possibility for a major (i.e., Mw ~9.5) tsunamigenic earthquake, precisely because these occur in the area only over large temporal intervals. Furthermore, modern research demonstrates that communities affected by socionatural disasters sometimes show limited social memory regarding these events, thereby increasing their vulnerability in the face of future hazards (, , , ).

Therefore, interdisciplinary scientific research on the occurrence of megathrust earthquakes and tsunamis over broad temporal scales, as well as the likely possibility for an Mw ~9.5 scenario in the mature northern Chile superseismic gap, needs to be urgently considered in regional and Pacific basin seismic and tsunami hazard assessment. The communication of results of such research through formal and informal means may aid in reintroducing the social memory of these events into local communities, favoring their resilience (). These results should also be used to recalibrate current hazard assessment policies, which, for coastal northern Chile, are mostly based on the historical tsunamigenic earthquakes. Our results suggest the need for increasing scientific interdisciplinary seismic and tsunami research to improve and foster hazard and disaster management for the future.

MATERIALS AND METHODS

Samples and radiocarbon dates

A total of 107 radiocarbon ages were obtained from different sites along coastal Atacama Desert of northern Chile, from Pisagua to Bahía Cisne: 64 correspond to marine fauna, 40 correspond to charcoal, 2 correspond to human bone, and 1 corresponds to organic sediments (table S1). The ages obtained from samples of terrestrial and marine origin were calibrated using the SHCal13 and Marine13 curves (OxCal software) (). The reservoir marine effect for northern Chile () was used to calibrate these dates (see table S1 for specific values). For human bone samples, the Marine13 curve was also used following the recommendations of previous works from stable isotope studies (), which show the high consumption of marine resources by prehispanic inhabitants of coastal northern Chile.

SPDs of radiocarbon dates from archaeological sites were calculated using a compiled database for all ages known from published material (, , ) and our own results for the coast of the Tarapacá and Antofagasta regions for the period 6500 to 1500 cal yr B.P. (n = 221). Radiocarbon ages were tested to determine statistical differences between them within individual sites, so as to avoid overrepresentation of specific time moments. Eighteen ages were excluded. The analysis was run with a total of 203 ages, just above the lower threshold recommended to produce reliable SPDs of calibrated radiocarbon dates. The resulting SPD constitutes an independent line of evidence, consistent with other archaeological data discussed here (site counts, site size, and density indexes).

Geological and archaeological sites description

We described 12 geological and archaeological sites located above present mean sea level (figs. S1 to S37, S41, and S42), with littoral outcrops and layers corresponding to uplifted and interpreted paleotsunami deposits with a minimum age close to ~3800 cal yr B.P., allowing the possibility for the analysis of multiple source scenarios from model results (fig. S44). Given that many archaeological sites in coastal northern Chile have been partially looted, we cleaned profiles left by different illegal excavations at the archaeological sites, giving us the opportunity to make detailed stratigraphic studies with a low investment in excavation. We also considered naturally exposed profiles, such as at the Paquica site. After cleaning several profiles in a site, whether naturally exposed or due to illegal excavations, we made detailed profile descriptions considering the texture, color, sedimentary structures, geometry, sorting, and archaeological remains of stratigraphic units. When available, we focused the chronological determinations on samples extracted immediately below and above stratigraphic contacts in deposits related to tsunami and/or uplifted littoral deposits. The determination of distinctive layers interpreted as tsunami deposits was based on stratigraphic and textural criteria from field observations. The complete description of these deposits from Hornos de Cal and Los Bronces () is included in this work. Similarly, the determination of uplifted littoral deposits was done following the previous descriptions from the Michilla and Playa Grande sites (), which are also included in this work.

Excavations and analyses at the Zapatero archaeological site

Between 2011 and 2017, eight naturally exposed profiles caused by alluvial erosion were cleaned at the Zapatero site (fig. S8). Stratigraphical observations were made, and samples for radiocarbon dating were obtained. Two 2 m by 2 m units (unit 1 and unit 2) were excavated at the central and deepest part of the site (fig. S29), although unit 2 was not excavated to the base of the deposit, but only as deep as a collapsed structure (structure 3). Last, two other structures were excavated at the southern end of the shell midden (structures 1 and 2). Unit 1, unit 2, and profiles 1 to 4 were excavated in the main mound (fig. S29), while structures 1 and 2 correspond to excavations on the southern edge of the site, where post–5700 cal yr B.P. structures were identifiable on the surface (fig. S29). Excavations from unit 1 and structure 2 have been detailed in previous work (, ). A backhoe trench and pits were also excavated from the current coastline up to the western limits of the site to assess the geological and geomorphological context (figs. S18 and S19).

Commonplace methods and techniques were applied in each archaeological excavation unit, following stratigraphic layers during excavation. Artificial spits of 5 cm were also used. Sediments were systematically dry-sieved using coarse (4 mm) and fine meshes (3 mm), with all artifacts and faunal remains recovered. Column samples (50 cm by 50 cm) were excavated from each unit (with the exception of structure 1) to characterize faunal remains, especially shells, and were subjected to water flotation and later sorted in the laboratory.

A total of 34 radiocarbon ages have been obtained from the site, making it the most intensively dated site along the northern Chilean coast. Samples for radiocarbon dating were obtained from different stratigraphic profiles, considering the main archaeostratigraphic units. A total of seven radiocarbon ages are available from unit 1, three for unit 2, one for structure 1, and five for structure 2. The remaining ages come from the cleaned profiles at different locations on the site.

All recovered archaeofaunal remains from the excavation units (shells, fish, and tetrapods) were quantified and subjected to anatomical and taxonomical identification to the lowest taxon identifiable, most often to species level, using reference collections from the study area and osteological manuals (–). A total number of identified specimens per taxa (NISP) of 32,830 faunal remains were analyzed for the stratigraphical units reliably dated to the Archaic III-IV, Archaic V, and Archaic VI periods in Zapatero (see below). We identified a minimum number of individuals (MNI) of 2661 specimens that included 68 families and species of marine invertebrates, marine vertebrates, and tetrapods, all presently found along the northern coast of Chile. Detailed descriptions of faunal groups and their representation have been reported in previous work (, , ). NISP and MNI were used to quantify faunal abundance (–). Lithics recovered at Zapatero were analyzed on the basis of macroscopic, petrographic, and technomorphological analyses of individual pieces, considering specific attributes for detached pieces and tools (). The quantification considered the number of artefactual specimens (NAS) and number of flake specimens (NFS), whether complete or as fragments (). To provide a more exact representation of the lithic production on-site, we estimated the minimum number of flakes (MNF), taking into account complete flakes and proximal fragments that contain the striking platform (, , ). Quantification of the minimum number of tools (MNT) considered complete specimens and nonconjoinable fragments. Last, the minimum number of artifacts (MNA) was obtained by adding MNF + MNT ().

We made preliminary comparisons between time periods at Zapatero in terms of frequency and density of artifacts and ecofacts. These comparisons were limited by several factors. Among them, the fact that no exclusively Archaic IV deposit has been excavated at the site, since these deposits include reworked Archaic III materials (). Furthermore, there are differences in the size and depth of excavated areas at Zapatero. Last, deposition rates and taphonomical processes have not been sufficiently determined for the site. Considering these limitations, to compare density of remains for each time period at the site, we standardized frequencies of lithic and archaeofaunal categories by the volume of excavated sediments (MNI/m3) (–). In the case of the lithic assemblages, we calculated the lithic volumetric density, i.e., total lithics recovered per cubic meter of sediment excavated (). Stratigraphic profiles, field notes, diagnostic artifacts, and radiocarbon ages were studied to correlate stratigraphic layers identified during excavation with specific cultural periods. Despite having more than 30 radiocarbon ages for the site, there are still several incompletely dated stratigraphic sequences. Hence, we used a proxy measure of deposition rate to estimate density of occupation for each time period: NISP and MNI per unit of time in the case of the archaeofaunal assemblage and NAS and MNA per unit of time in the case of lithic remains. Unit of time was estimated in terms of the time span corresponding to each cultural period.

Despite these limitations, results were consistent across all independent categories of data: NISP per time span, MNI per time span, NAS per time span, MNA per time span, NISP/m3, MNI/m3, NAS/m3, and MNA/m3. Considering MNI/m3 results were also consistent in terms of abundance, diversity, richness, and evenness indexes for archaeofaunal remains, since all of these parameters show their lowest values during the Archaic V Period (tables S7 to S11 and fig. S38). A nonparametric Kruskal-Wallis test was used to compare all these categories of data between time periods (Archaic IV-V and Archaic V-VI), supporting statistically significant differences between Archaic V and previous and subsequent periods in all cases (P < 0.005), except for tetrapods. Statistical analyses were performed using the Minitab 19 software. To compare between periods in terms of abundance, richness, evenness, and Shannon’s diversity index, we used the paleontological statistics (PAST) software ().

Comparison between residential sites in the Taltal-Paposo area was based on a density index (MNI/m3) and the frequency of lithic instruments per cubic meter (), since raw data for these sites were not available (table S11). The consistency of these results warrants our claims and behavioral inferences derived from the comparative data despite the abovementioned methodological limitations.

Results of a microarcheological and micromorphological study of unit 1 are available in previous work (). Here, we provide a summary of this study, in addition to previously unpublished results from quantifications of whole individual echinoid radioles in the thin sections for micromorphology and in bulk sediment samples collected from the micromorphology blocks (table S2). A total of six undisturbed blocks of sediment were collected for micromorphological analysis from the east profile of excavation unit 1 at Zapatero (figs. S30 to S35 and table S3). The location of the micromorphology samples tried to capture the contact between the archaeofacies. The samples were dried and impregnated with a mixture of resin, diluent, and catalyst at Spectrum Petrographics (United States). Standard thin sections (7 cm by 4.5 cm) were made with an average thickness of 30 μm. Samples were analyzed following previous guidelines (). Observations were made at magnifications from ×1.0 to ×20 using plane-polarized light and cross-polarized light with a Leica DM 2700 P petrographic microscope and Leica S9i stereomicroscope. Analyses were done at the Microarchaeology Laboratory of the University of São Paulo.

Excavations and analyses of other archaeological sites

Between 2008 and 2013, four units were excavated inside the prehistoric mining trench of SR-15 and in the exterior tailings, to obtain diagnostic and dateable archaeological contexts and to understand the stratigraphic variability of the sedimentary deposits. Excavation followed natural and artificial layers (10 cm) within each natural stratum. All excavated sediments were sieved (3 cm). In unit 1, located at the western end of the trench, 16 squares were excavated (each square is 1 m by 1 m). The excavation exposed the original mining trench reaching the bottom of the prehistoric exploitation at ~3 m from the present day surface. In unit 2, located on the central part of the trench, nine squares were excavated (1 m by 1 m). Our excavation has reached more than 6 m from today’s surface, but still the bottom of the prehistoric mine has not been reached. Two further units were excavated outside the main trench, where evidence of prehispanic mining was found under at least 60 cm of natural colluvium. Unit 3 was 2 m by 2 m and reached more than 2.8 m in depth, but the bottom of the prehispanic mining operation could not be found. Unit 4 was 0.5 m by 2 m and 1.5 m deep, but the base of the original mining excavation could also not be reached.

Precise position of uplifted littoral and tsunami deposits

The precise position of the GPS base was defined for the post-processing of the base line. We used the TBC software by Trimble for the base line processing. Two points were established at Zapatero Bay to determine the topographic map of the study area.

For the purpose of obtaining the mean sea level, we used tide tables and procedures given by the Chilean Hydrographic Navy Service. Specifically, we obtained the sea level height on the day of observation at Antofagasta port. Using the “Tablas de Mareas de La Costa de Chile, PUB. 3009,” we first applied corrections to the times and the heights series by using the methodology for secondary ports, called “Puerto patron” (fig. S43).

We designed a procedure specifically adapted for taking direct observations of the height of the sea using a prism pole. The corrections were applied to the secondary ports, which, in this case, correspond to “Rada de Paposo”.

Once we had corrected the time and the height of the secondary ports, we calculated the Sea Level Zero using the technique called “tide clock,” and we used the total station to obtain the reference to the mean sea level (table S13).

Tsunami modeling

The tsunami modeling is divided in three stages: seismic source, initial condition, and tsunami wave propagation. Seismic sources were generated with a stochastic method based on a power law in the wave number domain, namely, k−2 (). Given an area and a target magnitude, it was possible to create a plausible and physically consistent finite fault model. The studied area was gridded following the Slab 1.0 model (), which statistically approximates the interface of the subduction zones (, ).

The finite fault model was used as input data to calculate the coseismic displacement of the seabed. Because of the stochastic nature of the sources, the kinematic effects of the source process can be ignored, and the coseismic displacement was used directly as the initial condition for the numerical tsunami model. This surface was computed using Okada’s equation (), which calculates the elastic deformation of a half-space due to a shear fault. The horizontal advections contribute to the final vertical displacement, as pointed out in previous work ().

Last, the selected model for the tsunami propagation and run-up is NEOWAVE (, ). This model numerically solves the nonlinear shallow water equations in a staggered finite difference scheme. It also included a term accounting for the nonhydrostatic pressure. The numerical model requires two main inputs: the initial conditions, which are described above, and the bathymetry data. In this case, bathymetry was set at one level only of 30 arc sec from the General Bathymetric Chart of the Oceans website. The interaction with the solid boundaries (coastlines) was numerically simulated as a full reflection boundary condition, whereas in the open boundaries, a radiation condition was set. We selected four models from the total number of simulations run that best fitted with our geological observations (fig. S44).