The Mitla Landslide, an Event That Changed the Fate of a Mixteco/Zapoteco Civilization in Mesoamerica

Before the arrival of the Spanish conquerors, Mitla was the second most important city in the valleys of Oaxaca, México. However, the ruins that are visible today do not seem to match the size of a city of more than 10,000 inhabitants. Geological and geophysical studies suggest that part of the city was covered by the deposits of a dry landslide likely to have been caused by an earthquake with a magnitude that could vary between 6 and 7.This landslide is monolithological, which is why two geophysical methods were used in order to evaluate its geometrical characteristics and to suggest the possible existence of archeological remains under the landslide.


Introduction
Many hypotheses regarding the collapse of Mesoamerican cultures have been linked to wars, conquests, epidemics, or climate changes. However, Mesoamerica is located in a context presenting exceptional seismic, volcanic, and meteorological phenomena. Past and current events in México (Tehuantepec, Pinotepa Nacional or Morelos-Puebla earthquakes ∼M8, SSN) indicate that tsunamis, earthquakes, and volcanic eruptions countered the development of many important civilizations. For instance, the collapse of Maya culture is insistently related to a large drought, yet very few theories associate these human extinctions with exceptional events. Much of this disaster information was written in the codices. These documents contain data with place and time specifications, also providing information on natural disaster scenarios. For example, in the Telleriano-Remensis codex, tlacuilos recorded approximately 12 seismic events that include glyphs that suggest a scale of environmental effects that appears similar to ESI2007 [1].
The location of Mesoamerica compels us to think that those cultures witnessed very important seismic and volcanic events. For instance, the birth of the Tarascan empire is clearly linked to changes in the landscape due to an earthquake and to the birth of Capaxtiro volcano [2]. Besides the areas towered by volcanos, earthquakes caused important imbalances or prophecies, such as the fourth bad omen described by León Portilla [3], which warned the Aztecs about the arrival of the Spanish conquerors (The burning lake, El Lago en Llamas): "Texcoco Lake burned, its boiling waters swelled and swept the houses, killing many". This description is certainly very close to an account of the formation of a tsunami caused by an earthquake.
On the other hand, the ongoing seismic activity in the State of Oaxaca brought about a particular architectural practice: in order to resist violent earthquakes, the churches have very thick walls, huge buttresses, and short towers ( Figure 1).
The study of the effects of important earthquakes can be carried out according to different environmental or anthropogenic results. In present case the research focuses on the geophysical analysis of the Mitla Landslide. As in this case the collapse was caused by an earthquake, the techniques selected were ground-penetrating radar and electrical tomography and cross-correlations of seismic ambient noise records. The two last were those that yielded the best results. Eventually the use of the radar was not a good option due to the homogeneity of the landslide and for the little penetration of this, which in this case is monolithologic. Conversely, electrical tomography was an adequate method. The use of seismic noise and the advantage of a regional earthquake register allowed us to find a clear lithological contrast.

Tectonic Setting of Oaxaca
Oaxaca is one of the major earthquake zones of México, due to the Cocos plate subduction under the North American plate and the seismic activity related to intraplate structure that causes shallow quakes, such as the Xalapa earthquake of 1920. Oaxaca has been affected by a large number of seismic events with magnitudes above 7 (Table 1, data provided by the Servicio Sismológico Nacional (SSN), UNAM 2018), so much so that Fray Juan de Córdova's Vocabulario en lengua zapoteca (1578) contains words to designate the tremor of the earth: xoo, but also a name for the deity of earthquakes: Pitao xoo. Beside subduction earthquakes, Oaxaca presents intraplate faulting related to the limits of tectonic terranes expressed in the Central and Tehuacan valleys. The Central valleys of Oaxaca are part of the Oaxaca fault system, which is formed at the borders between the Zapoteco and Cuicateco terranes (Figure 2(a)), and besides its normal component has always been thought to have a horizontal component in arguable direction (Figure 2(c), [4].
The seismicity reported by the Servicio Sismológico Nacional (SSN) and presented in Figure 2(c) corresponds to the intraplate events zone in the south-east zone of México named by Zúñiga et al. [5] as NAM and matches the intraplate seismicity of the southeastern zone of México that is not related to the volcanism of the Trans Mexican Volcanic Belt; its maximum depth has been measured at 20 km.
The focal mechanisms presented in Figure 2(b) correspond to quakes of magnitude over 7 and were taken from the IRIS (Incorporated Research Institution for Seismology) consortium. Most of these focal mechanisms are associated with the zone of strong coupling between the Cocos and North American tectonic plates (interplate, reverse faulting) and others to rupture processes of the subducted Cocos plate (normal faulting).
Recently, in 2017 and 2018, major events occurred again in Oaxaca, the first on September 7 with magnitude 8.2 Mw and epicenter located in the Gulf of Tehuantepec, 137 km southeast of Pijijiapan (Chiapas), and 46 km deep. This event is mainly related to normal faulting (79 ∘ dip) that caused the complete rupture of the lithosphere in a zone related to the seismic gap of Tehuantepec, in which no event of this magnitude had taken place since 1787, when an earthquake of Mw≈8.6 occurred [6].
The second event happened on Friday February 16, 2018, of magnitude 7.2 Mw and epicenter 11 km from Pinotepa Nacional, Oaxaca. Both events caused material and human losses.
It is important to mention that regarding recurrence periods for the central zone of the State of Oaxaca, Zúñiga et al. [5] report periods of 6 years for quakes with magnitudes above 5 occurring in the NAM zone (shallow events with depth under 20 km), periods of 109 years for magnitudes above 7 taking place within the subducted Cocos plate (intraplate events with depths of 40 to 180 km), and periods of 37 years for quakes of magnitudes above 7.5 occurring in the coupling zone of the Cocos and North American tectonic plates (subduction area).
On the other hand, previous refraction and reflection studies across the Oaxaca fault [7,8] did not confirm that it runs as deep as presumed. Supporting this result, the regional magnetotelluric study by Jödicke [9] did not detect the expected large electric contrast across the Oaxaca fault at crustal depths along a transect passing south of Oaxaca City. By analyzing the electrical impedance along a regional transect, Jording et al. [10] concluded that the major structural change at depth in this region is displaced northeastwards of the Mitla Valley, about 30 km east from the generally accepted Oaxaca-Juárez terrane boundary. A more detailed MT study of the region [11] concluded that the Oaxaca fault system does not appear to penetrate deep into the crust but remains relatively shallow (5 km). Other results [12] also suggest that the contact between the Oaxaca and Juarez terranes at crustal depths occurs along a SW-dipping interface south of Oaxaca City [13].

Study Zone and Landslide
The town of Mitla is located south of Oaxaca City (16 ∘    where sometimes the N-S drainage system is modified by NW-SE structures (Figure 2(c)).
The word Mitla or Mictlán is nahuátl and means "Place of the Dead" or "Underworld". In Zapotec it is called Lyobaa which means: "place of graves" or "place of burials", from the etymology Lio: Place and Baa: grave or burial [14]; in mexica it remained Mictlán, "Place of the Dead", or "place of many corpses", and was hispanicized to Mitla. Contrary to Monte Albán, the heyday of Mitla occurred between 950 and 1521 AD. The settlements of Mitla have been estimated to be over 7,000 km 2 large, with a population of around 10,000 inhabitants. One of its main characteristics are constructions decorated with elaborate mosaic fretwork (grecas) featuring variations of the same geometric design, and the crossshaped tombs that have been found under the palaces, in which important people and priests were probably buried. The archeological complex includes several structures, two of which were covered with Christian buildings, the most important of which was constructed in 1590. The San Pablo church was erected in one of the prehispanic courts of Mitla using stones from the temples themselves. The basic architectural form of the few pyramids shows paraseismic techniques, visible both in the Las Columnas site as in the Grupo Arroyo, where notches appear in the frames of all the entrances to the rooms.
During our field work the avalanche deposit was identified for the first time in all the surroundings of the archeological remains of San Pablo Villa de Mitla. Previous works mapped the avalanche as part of pyroclasts deposits [15]. Its distal facies seem to have covered part of the pyramids, as shown hereunder and on

Methodology
The methodology of this study focused on 5 main guidelines: (a) the first step was to find previous research on Mitla in the fields of Archeology and Geosciences; (b) subsequently, a 50 cm-definition terrain model was mapped with a drone; (c) International Journal of Geophysics   International Journal of Geophysics geology and all the data regarding geophysical, morphostructural; (d) stratigraphical characterization, structural geology in the Sierra Calaveras and (e) the geomechanical study of the landslide were based upon this model. Various columns were extracted from the landslide to study its geology. Once the lithology was understood, the area was mapped in detail.
Several methodologies were applied regarding geophysics, such as ground-penetrating radar (Mala) and electrical and seismic tomography. The first method was not adequate due to insufficient penetration and scarce lithological contrast. The seismic and electrical tomographies, however, yielded valuable results.
Reconnaissance methods, which mainly include remotesensing and aerial techniques, geological and geomorphological mapping, and geophysical and geotechnical techniques, had to be adapted to the characteristics of the landslide. According to McCann and Foster [16], a geotechnical appraisal of the landslide's stability has to consider the three following issues: (1) the definition of the 3D geometry of the landslide with particular emphasis on failure surfaces, (2) the definition of the hydrogeological regime, and (3) the detection and characterization of the movement.
The method for seismicity was to obtain velocity profiles using horizontal-to-vertical spectral ratios (H/V ratios) of microtremors [17,18] named Nakamura method, and inversion with the Spatial Auto-correlation Method [19] and with the H/V [20] based on ambient noise measurements carried out with 3 long-period Trillium Compact 120s sensors with frequency responses of 120s to 100 Hz, and their respective Ref Tek 130S three channel and 24 bits-resolution third generation recorders. Regarding microtremors, it was taken into account that the energy sources are multidirectional [18]; thus, the direction of maximum movement is unknown. Moreover, the use of this kind of technique does not cause environmental impact and yields precise results.
Microtremors, also called ambient vibration noise, seismic microtremors, or background seismic noise, are random vibrations generated by natural or artificial sources. The aim of the Nakamura method and the SPAC technique was to determine velocity profiles in the northern region of the archeological zone of Mitla in order to determine the depth of the base of the landslide.
The H/V technique was developed by [18] and is used for obtaining direct estimates of site-specific response during an earthquake (e.g., [21,22]), with relating geological and geophysical properties of Wells together with seismic records analysis on the various geological site conditions. It was hypothesized that the vertical component of the ambient noise at the ground surface keeps the characteristics of basement ground, is relatively influenced by Rayleigh wave on the sediments, and can therefore be used to remove both the source and the Rayleigh wave effects from the horizontal components. It is effective in identifying the fundamental resonance frequency of a sedimentary layer, where the depth of a soft soil layer with resonance frequency and shear wave velocity Vs is given by H = Vs/(4 * F).
The H/V spectral ratio of microtremors measured anywhere confirmed its ability to estimate the predominant frequency and the amplification factor. The result of the estimation is stable for the measured time and season (Nakamura, 1989).
The seismic noise records were visually inspected to select signals that are complete and that had a good signal-to-noise ratio. Two software packages were used for data processing: SAC (Seismic Analysis Code) and Geopsy (Geophysical Signal Database for Noise Array Processing) (SESAME WP05) [23]. The results of these programs were compared because they have a different smoothing function in obtaining H/V [24].
H/V spectral ratios were calculated for the three components of each record. This was done automatically when estimating the H/V ratio with the Geopsy software. For the calculation of H/V for each selected window, a smoothing function was applied with a bandwidth coefficient of b = 40 and a 5% cosine taper-window. This type of smoothing function employs a different number of points at low and high frequency; its use is strongly recommended for frequency analysis [25].
Using the SAC (Seismic Analysis Code) software, spectral ratios were determined in 81.92 seconds windows and a smoothing with a Von Hann window [26]. The results of both programs did not show significant differences.
The method known as SPAC (spatial autocorrelation method) was first introduced in 1957 by [27]. The SPAC method considers the dispersive characteristics of surface waves in a stratified medium [28]. When only the vertical component of seismic noise record is used, we can easily assume that the signal has a high content of Rayleigh waves.
In this method, it is proposed that from simultaneous microtremor measurements in some geometric arrangements of seismic stations, it is possible to obtain the dispersion curve by calculating the coefficient of spatial autocorrelation [29].
We represent harmonic wave circular microtremor frequency by u (0, 0, , t) and u (r, , , t), which are observed in the center C (0,0) of the array and point X (r, ) on the circle of radius r. Then spatial autocorrelation function is defined as ( , , ) = (0, 0, , ) ( , , , ) where u(t) is the average value in the time domain. The spatial autocorrelation coefficient for one angular frequency , p ( , r), or simply autocorrelation coefficient is defined as the power spectrum of one station within a spatial arrangement (middle circle): Therefore, the spatial autocorrelation coefficient at a frequency is related to the phase velocity c(f) through the Bessel function of first kind and zero order o . The phase velocity calculated from the argument of the Bessel function generated (2). For the SPAC technique, we used Geopsy software [23]. Prior to making rings, correlation graphs were generated, and the dispersion curve was then adjusted. Seismic noise measurements were taken during two days on the archeological site of Mitla. For each day, a triangular International Journal of Geophysics   array with three seismic stations and with seismic records with durations of 10 hours (array 1) and 6 hours (array 2) was built ( Figure 4).

Geoelectrical Tomography
Geoelectrical imaging is geophysical technique that was developed on the basis of physical methods that are helpful to reveal the presence or absence of buried bodies and structures that cannot be seen with the naked eye, but that can be detected because of their distinctive physical properties of being resistive to their surroundings (e.g., [30][31][32]. The aim of geoelectrical mappings is to determine the distribution of resistivity in the subsoil through measurements carried out on the surface. The true resistivity of the subsoil can be estimated from these measurements. In this 8 International Journal of Geophysics sense, the resistivity of the earth is related to geological parameters according to its mineral content, fluid saturation, and rock porosity. Geoelectrical studies have been used for some decades now for hydrological and mining research and for geotechnical research. However, more recently, these methods have become more sophisticated in mapping the subsoils, and they have increasingly been used to locate structures and characterize the subsoil.
The aim of this study is to obtain a 2D model of true resistivities from the pseudosection of apparent resistivities obtained from the field. This is why an inversion program is necessary.
On the other hand, samples of the landslide sediments were submitted to granulometric analysis and evaluation of their petrophysical properties. Unfortunately, neither organic soils, nor pieces of carbon were found to date the landslide.

Geology of the Mitla Area.
The stratigraphical column of the Mitla zone is composed of a basement of Precambrian rocks from the Zapoteco terrane covered by Mesozoic carbonated units veiled by pyroclastic flows that constitute Sierra La Calavera. On these sequences and on the hanging wall, huge Pleistocene lacustrine basins developed as part of the large valleys of the study zone, which are actually the hanging wall of the normal NNW-SSE fault carved by other faults following the same direction and which could be responsible for the Mitla Landslide (Figure 3).
The main structures are oriented NNW-SSE and clearly present the geometry of normal faults with strik-slip sinistral component. This fault, which we henceforth call Calaveras fault, is part of the large structures that are the boundaries of the tectonic terranes, especially those that limit the Zapoteco and Cuicateco terranes (Figures 2 and 3).
The ignimbrite sequence of rhyolitic composition is formed of at least two large packages; the inferior package corresponds to the hanging wall of the normal Calavera or Mitla fault outcropping at La Fortaleza and Mitla (Figure 3). The second is the package of Sierra La Calavera and is made of at least 4 thick ignimbrite packages, some of which are welded and others crumblier, but all with Q, PL, and mica crystals; some levels are tightly welded at different points of the column. These were dated between 14.4 ± 0.4 My and 15.48 ± 0.2 My, Miocene [15].
Above the ignimbrites of the hanging wall of the La Calavera fault is a thick landslide deposit formed of different levels but nonetheless of monolithologic composition, i.e., formed only of block and rhyolitic rocks matrix. The blocks range from gravel to blocks of more than 5 m in diameter, some of which are angular and others semirounded. More data regarding this landslide will be given hereunder.
On top of the landslide and the ignimbrites are lahar deposits, varying from blocks and matrix to hyperconcentrated lahars.
Fluviolacustrine deposits were observed in río Mitla formed by conglomerates, sands and clays, all of which heavily faulted. Large Pleistocene vertebrates were reportedly discovered in these deposits [15].

Geometry of the Collapses and Description of the Deposits.
In a detailed elevation model, two large landslide flows can be seen opening towards the south in the Sierra La Calavera, which is formed by tabular deposits of ignimbrites tilted towards the north. Because of its provenance, the most visible of the flows is the Mitla Landslide. This wide scarp and deposit bundle is clearly observable and delimited. The main scarp has a maximum diameter of 1350 m and a slope relief of 390 m. The landslide ran over more than 2500 m from the basis of the main scarp to the toe. The second main scarp, located NNO of Mitla, is less regular, both in its main scarp and in its body ( Figure 5). Both main scarps are located within the footwall of the normal fault that caused the Mitla depression.

Geomechanics and
Granulometry. An average per size was obtained, of samples from the avalanche front, from the granulometric data using the sizes of the wire mesh cloth of the sieves used to observe the behavior of one granulometry ( Figure 6). As in the case of individual strata, the line does not break, which reveals a fair distribution of the material and very little space among particles.

Φ Value.
With an elevation model presenting curve details every 50 cm, we obtained the value of the H/L relation for a length (L) of 1.633 km, and a drop height of 0.640 km, yielding Φ=0.39, which means its value is characteristic of dry landslides, i.e., of 0.5.

Slope Mobility.
Reference [33] introduced an indicator regarding the mobility of slopes called excessive travel distance. This parameter (Le) expresses the horizontal distance traveled in excess by the event, superior to what could be expected if the movement was that of a nonlubrified event of rigid mass, which would move across a tilted plane with a normal coefficient of friction of tan (32 ∘ ).
Le = 1033 As can be observed, the distance covered in excess indicates the movement of the body in the specific conditions mentioned above, and comparing it with the true distance of the slope we can observe that the difference is small, indicating scarcity of water. In cases wherein, water is presented as triggering factor; the difference can be four times higher.

Geophysics.
In order to know the geometry of the landslide, three 2D tomography lines were acquired with a= 4 m ( Figure 8). The vertical resolution obtained was of 2 m and the horizontal resolution of 4 m. Lines 1 and 2 were 108 m long with NE-SW orientation. Line 3 was 156 m long and oriented SE-NW (Figure 7). Data obtained so far reveal a landslide body of low resistivities (11)(12)(13)(14)(15)(16)(17) in blue and in the southern part of the profiles, between 4 and 5 meters from the ground, there is a body of average resistivity (22) colored green and yellow, interpreted as part of the ignimbrite basement or possible remains of pyramids constructed with ignimbrites.  6.6. Seismic Noise Records. On the assumption that the derived dispersion function represents the first mode of the Rayleigh waves, we derived a best-fitting velocity model that is consistent with the estimated velocity values. Figures 8  and 9 show the Rayleigh waves dispersion function for each seismic array. Figure 10 presents the H/V spectral ratio for stations 1 and 3 for seismic array # 1. Here, the H/V rotate tool was used to obtain the H/V in the horizontal plane, i.e., as azimuth function, from any type of 3D vibration signals (ambient vibrations, earthquakes, etc.). The black curve represents H/V geometrically averaged over all individual colored H/V curves. The two dashed lines represent the H/V standard deviation. The grey area represents the average peak frequency and its standard deviation. The frequency value is at the limit between the dark grey and light grey areas. The parameters were frequency sampling from 1 Hz to 15 Hz, global time range from T0 to the end of the time series and time windows to 81.92 s.

Inversion of the H/V Spectral Ratio.
The inversion of the spectral ratios allows us to obtain a velocity profile for each measurement site. In Table 2 we show the fundamental frequencies for each seismological station and the classification of the type of soil associated with the site based on the NEHRP (National Earthquake Hazards Reduction Program) classification [34].
It was determined that the depth of the base of the landslide is of 50 meters below array 1 (region south of the

Discussion
The Mitla avalanche opens a great discussion in the field of Geoarchaeology. The geological investigations give evidence that this event occurred in historical time that is during the Postclassic and before the arrival of the Spaniards. The important city of Mitla was found in decline during the arrival of the Spaniards.
The geophysical exploration data show important aspects, as the electrical tomography studies propose the existence of an important anomalous structure under the avalanche that could be part of the remains of the prehispanic Mitla. The elements that compose it were found between 5 and less of 20 m deep. On the other hand, with the earthquakes data and environmental noise, an avalanche thickness of 60 m was estimated, which allowed a correct calculation of the area and volume of the material removed (area of 2.79 km 2 and volume of 0.1674 km 3 ). In this same field of the relation between the volume or area removed by the collapse and the magnitude of the earthquake, two published works have been analyzed; one was by [35] and the other by [36]. Both make a great contribution to obtain the relationship between the destabilized areas and the Ms or Mw. However, in all these works they analyzed zones of Central America where the climate and the soils make up very different scenarios from Mitla. Consequently, it is likely that the magnitude proposed in our work does not correspond exactly to the real one, because the avalanche was in dry condition during the collapse, the geometries of the strata bedrock were not favorable to the destabilization, and the local climatic conditions have not allowed a great soils development. Anyway, here there is the morphological evidence that the Calaveras fault is part of an active fault system, as shown from the drainage networks that give the NNW-SSE structures a normal left-lateral component (Figure 2(c)). Despite these differences, we have considered keeping Keefer's equation.
In this paper we have only estimated the material volume of the Mitla avalanche. However, in later works we will make an estimate of the two avalanches that could have been generated by the same earthquake (Garduño-Monroy et al., in press) Starting from a maximum volume of 0.1674 km 3 for the landslide, the earthquake that caused the ground shaking up to the slope collapse could have had a maximum magnitude of 7.3 according to the relationship of [35,36].
If we consider the length of the Calaveras fault (with a minimum length of 10 km and maximum of up to 40 km, Figure 3), it would be potential to generate an earthquake of magnitude in the range of 6.2 to 6.9 Mw, based on the relations proposed by [37] or [38], respectively.
So then, two possible sources were feasible in the generation of this avalanche. The presence of the geodynamic framework of Oaxaca with the occurrence of earthquakes between 7 and 8 Mw could suggest that a comparable event could have generated the studied avalanche. On the other hand, if we consider also the morphological evidence of the Calaveras fault, an earthquake between 6 and 7 Mw could have easily generated an avalanche of this type if the fault was active in the period of the landslide event. Both scenarios are plausible, but now the most important thing is to see the possibility of starting new geoarchaeological explorations and consider that these segments of the Oaxaca fault system are active.

Conclusions
This research indicates that the Mitla Landslide is a collapsed body formed of ignimbrite blocks and matrix from the Sierra La Calavera. According to its morphology, the geotechnical characteristics, and the geophysical data interpretation, the landslide was provoked by an earthquake of a magnitude in the range from 6.2 to 7.3 Mw. The observations on the field suggest that the earthquake could have been caused by the La Calavera fault, which is a normal fault with small left component and is part of the Oaxaca fault systemoriented NNW-SSE to N-S. But we also do not rule out that this avalanche was instead generated by the frequent seismic activity of the Cocos plate. These geological facts (from regional to local tectonics), and the age of the landslide, could introduce a new scenario regarding seismic hazard in Oaxaca, because beside the intense seismic activity caused by the Cocos plate, more shallow earthquakes could also occur. On the other hand, the fact that the Mitla fault (Calaveras) is currently active indicates that many of the fault segments are favorably oriented to the current stress field.
There is no precise age established for the landslide occurrence; however, the event presumably damaged the pyramids of Mitla, and a large part of the pyramids are probably located under the avalanche deposit as evidenced by the outcomes of this preliminary investigation. Considering that if the population reached 10,000 individuals, the extension must have been larger than what can be implied by the current remains.
From the point of view of geophysical exploration, the method of electrical tomography and the study of earthquakes and environmental noise were good and promising tools in the Geoarchaeology field of research, even if the procedure can be better calibrated in synergy with other survey techniques. It was not the case of the use of Georadar that it failed to provide great information in the highdetail characterization of the local subsoil and materials discrimination in depth as necessary for this type of research. For future work it is recommended to make a thinner network with electrical tomography in parallel with seismic profiles.

Data Availability
All the geological and geophysical data have been generated by the authors of this work.

Conflicts of Interest
The authors declare that they have no conflicts of interest.