Ancient Maya wetland fields revealed under tropical forest canopy from laser scanning and multiproxy evidence.

We report on a large area of ancient Maya wetland field systems in Belize, Central America, based on airborne lidar survey coupled with multiple proxies and radiocarbon dates that reveal ancient field uses and chronology. The lidar survey indicated four main areas of wetland complexes, including the Birds of Paradise wetland field complex that is five times larger than earlier remote and ground survey had indicated, and revealed a previously unknown wetland field complex that is even larger. The field systems date mainly to the Maya Late and Terminal Classic (∼1,400-1,000 y ago), but with evidence from as early as the Late Preclassic (∼1,800 y ago) and as late as the Early Postclassic (∼900 y ago). Previous study showed that these were polycultural systems that grew typical ancient Maya crops including maize, arrowroot, squash, avocado, and other fruits and harvested fauna. The wetland fields were active at a time of population expansion, landscape alteration, and droughts and could have been adaptations to all of these major shifts in Maya civilization. These wetland-farming systems add to the evidence for early and extensive human impacts on the global tropics. Broader evidence suggests a wide distribution of wetland agroecosystems across the Maya Lowlands and Americas, and we hypothesize the increase of atmospheric carbon dioxide and methane from burning, preparing, and maintaining these field systems contributed to the Early Anthropocene.

Precolumbian Maya civilization persisted from ∼3,000 to 1,000 or 500 calibrated years before present (BP), building vast numbers of cities, farms, roads, and reservoirs. Over the last decade, lidar (light detection and ranging) imagery has greatly expanded our estimate of ancient infrastructure that ground survey had laboriously identified for more than a century (1, 2). This massive built environment indicates large ancient populations over a wide area with complex economies that required large-scale and diverse subsistence strategies (1, 3). Many studies since the 1970s have used a wide array of tools to reconstruct food and farming systems, such as indigenous swidden or milpa systems and more intensive agroecosystems on terraces and in wetlands (4). Research in the 1970s (5, 6) heralded the growing evidence of ancient agricultural intensification, which many volumes examined in depth (7–9). Evidence has matured steadily since then for crop, soil, and water management systems (10, 11). Wetland fields and canals make up one such system, and similar systems are widespread across indigenous Mesoamerica and the neotropics (12–19), as well as in ancient China (20), Angkor Wat (21, 22), New Guinea (23), and in modern Africa (24). Wetland systems have been keys for human subsistence, and their sediments and stratigraphy can provide evidence for human responses to droughts, floods, and sea level change (21, 25–27), as well as global scale human-induced environmental change that may amount to an Early Anthropocene (28, 29). Nonetheless, research on Maya wetland farming that started with fervor in the 1970s (30–40) languished in controversy and doubt after the 1990s (27, 41). This decline had multiple causes, including uncertainty about wetland formation processes, chronology, extensiveness, and importance for Maya food production (13, 33–38). Wetlands under tropical forests are also inherently difficult to excavate, and widespread recent plowing, draining, and deforestation are destroying many relict field systems (28). Fortunately, tropical forest cover is still extensive in some areas such as Guatemala’s Petén, where a recent Science article estimated 67 km2 of seasonally wet, canalized fields based on lidar mapping with some ground validation beneath the forest canopy (1). These Petén field systems and canals will require further extensive validation, excavation, dating, and multiple proxy evidence to confirm their uses, chronologies, and extents over Maya history.

We report here the verification of widespread ancient Maya wetland fields in the Rio Bravo watershed of Belize based on lidar data and many excavations with multiple lines of evidence for cultivars, formation, and chronology (Fig. 1). These ancient wetland systems occur in four main areas within the watershed. We focus on the area with the most evidence thus far, the Birds of Paradise (BOP) fields (26, 27), to present a synthesis of formation, use, and chronology based on 42 radiocarbon ages, including 15 from six wetland excavations (Fig. 2 and ; these include all radiocarbon ages amassed for the BOP fields and canals for this project). The extensiveness of the field systems based on lidar and the multiple lines of evidence satisfies the doubts that led to the decline of research on Maya wetland agriculture. We provide clear evidence of ancient Maya construction, subsistence, and chronology along a river system connecting scores of ancient Maya centers and the Caribbean coast beyond. These newly verified wetland fields also add to the growing evidence for early and sometimes intensive human use of tropical forests (42), tropical wetland extensiveness (43), and temporal correlations with vital questions of Maya civilization. For example, the formation and use of the wetlands coincides in time with the Maya Late Preclassic (2300–1700 BP) droughts, Late Classic population expansion (1400–1120 BP), and Terminal Classic droughts (44–47) and the Terminal Classic abandonment and major population realignments (∼1000–500 BP) (48). Both the Late Preclassic and Terminal Classic are periods of major population and trade-network realignments, and both have growing evidence for significant drought periods. Our findings here also begin to place Maya agriculture into a framework of Late Holocene global change (29) based on the expansion of wetland fields from ∼1800 to 900 BP through burning and canalization and, thus, their greenhouse gas (CO2 and CH4) emissions.

Fig. 1.

(Upper Left Inset) Location of study area. (Left) Sentinel-2 image overlayed with lidar-derived DEM. (Right) (E) Chan Cahal North shown using the lidar intensity image, Chan Cahal West shown using the lidar bare-earth model (D), a portion of BOP (C), and the central Rio Bravo floodplain (B) using the the lidar bare-earth model. We discuss Chawak But’o’ob zone (A) in the text.

Fig. 2.

AMS dates chart and conceptual model of wetland formation ( for AMS ages).

Research from the 1970s through 1990s on wetland fields 40–50 km north and closer to sea level (5 meters above sea level [masl] and below) than the current study area (∼8.1–35.3 masl) suggested chronologies from the Archaic period before ∼3000 BP to the Terminal Classic period ∼1000 BP. Conflicting results created still-unresolved debates about these chronologies, natural versus anthropogenic formation factors, and whether wetland systems were widespread (30–41). Research since 2001 on the Rio Bravo watershed has revealed different geographies, timing, and uses of ancient Maya wetland fields and canals (26–28). These relict field complexes cluster around ancient Maya nonelite residential groups that occupy small rises (20–40 masl) above the wetlands. Larger Maya cities occupy the higher adjacent escarpments (∼150 masl), and many show evidence of imports from mollusks, turtles, and fish, which we have found in abundance in ancient wetland fields and canals. Intensive excavation along with earlier, passive remote sensing reported evidence for chronology and past uses of a small fraction of these wetland field systems (26–28).

Rio Bravo Wetland Complexes

The ancient wetland farming systems occur in four main zones of the Rio Bravo watershed (Fig. 1). In the upper watershed (Fig. 1), investigations of the ancient Maya village of Chawak But’o’ob suggest that intensive terrace and wetland field systems arose and declined in the Late/Terminal Classic (1200–1050 BP) (49). Near the Rio Bravo’s mouth in groundwater fed wetlands around the Maya residential group of Chan Cahal (Fig. 1), maize agriculture started in the Archaic (∼4000 BP) (26) and intensive, wetland agriculture arose as early as the Late Preclassic (∼1800 BP) although most evidence here dates from the Late/Terminal Classic (50, 51). In the lower middle Rio Bravo floodplain, previous work indicated the BOP area (Fig. 1) mainly also arose during the Late/Terminal Classic and declined by the Terminal Classic and Early Postclassic (1300–900 BP), paralleling the later history of two adjacent ancient Maya centers (Fig. 1) (26, 27, 52). We also report a rediscovered wetland field complex (Fig. 1). All of the excavations discovered exclusively ancient Maya artifacts except at Chan Cahal (Fig. 1), where modern cultural items occurred in the plow zones of deforested pastures.

Results

Lidar Mapping and Areal Extent.

For decades, our research group has used many types of passive remote sensing with little success to discern wetland fields previously discovered from small planes and from ground surveys under forests of the Rio Bravo Conservation Management Area. Our objectives were to test the reasons that Maya wetland research had languished, including doubts about whether it was a widespread phenomenon. We hypothesized that lidar could help quantify the aerial extent of ancient Maya intensive agriculture. In June 2016, we engaged the National Center for Airborne Laser Mapping (NCALM) to collect airborne lidar data of the wetlands and surrounding areas (Fig. 1). NCALM employed a Teledyne Optech Titan multispectral lidar sensor (53) onboard a Piper Chieftain (PA-31-350) aircraft flying at 570 m above the ground to map an area of 274.6 km2. The sensor lasers were fired at a combined repetition rate of 525 kHz with a scan angle of ±30° and a scan frequency of 25 Hz. The processed point cloud yielded 11.065 billion returns from 6.516 billion laser pulses (∼1.7 returns per pulse), on average 23.7 pulses/m2 and 40.1 returns/m2 with a vertical accuracy estimated at of 4.4 cm root mean square error (RMSE). The classified point cloud was used to produce digital elevation models (DEM) as well as active multispectral images based on the return intensity at three different laser wavelengths (1,550, 1,064, and, 532 nm).

Lidar technology has been a remarkable boon to Maya archaeology (1, 2, 53–59), and to other areas of ancient civilizations in tropical forests, such as Angkor Wat (60) because of lidar’s ability to penetrate and to map the ground through gaps in the forest canopy. In tropical forests, lidar has helped to identify massive ancient infrastructure. This present study in the Maya region couples the quantification of large-scale, wetland canal complexes through lidar with multiple lines of verification for the types and chronology of these ancient farming systems. In addition, this study complements lidar’s spatial and spectral data successfully to document these agricultural features under different types of land covers.

Because of the topographic and spectral characteristics of these ancient agricultural features and their varied vegetation cover, we were able to use spatial (topographic) and spectral (intensity) lidar data together. In the study area, the key source reports tropical forest cover is ∼77%, savanna ∼6.1%, agriculture ∼16.3%, and wetlands ∼1.1% (61). In regions with tree cover, we were able to use spatial imagery derived from DEM to map the topographic signatures of the features. The DEM represents a bare-earth model of the surface without vegetation. In regions of savanna and agriculture, we were able to use imagery derived from the lidar multispectral intensity (LMI) at the 1,550-, 1,064-, and 532-nm wavelengths to map and identify the spectral signatures of the agricultural features in similar fashion of what has previously been done with satellite imagery. This represents a fortunate set of conditions given that all of the previous archaeological lidar research in the region has primarily used the spatial data of lidar to identify anthropogenic topographic signatures, with a few attempting to use intensity as an additional feature (62, 63) but without success documenting large-scale features as presented here.

The topographic and spectral lidar imagery indicated wetland agricultural systems at all four zones where ground survey had previously discovered them. DEMs, however, helped define many more wetland complexes at two of the sites with thick forest canopy: BOP and Chawak But’o’ob (Fig. 1 and ). DEMs also helped identify a zone, the largest, in the mostly unexplored central Rio Bravo floodplain (Fig. 1). In addition, the LMI facilitated defining more wetland complexes at Chan Cahal (Fig. 1). Here, deforestation and growth of dense tropical pasture grasses for the last 50 y (50) limited DEM fidelity.

BOP contains the largest verified area of Maya wetland fields visible on lidar (Fig. 3 and ), with a ∼5 km2 area of canals and fields including the original ∼1 km2 area of BOP complex discovered by aerial photography (Fig. 3). Many of the BOP canals under canopy are visible in DEM (Fig. 3) because laser penetration mapped the surface topography of canals and fields well. The BOP wetland complexes visible in the topographic imagery include ∼71 km of 3-m-wide canals and a median field size of ∼2,790 m2 (mean of ∼7,031 m2) (Figs. 1 and 3). The field and canal elevation ranges from ∼9.5 to 17 masl. Canal length may be longer and field size may be smaller than these lidar estimates because we counted the area between canals as fields although there may be undetected canals in some of the larger fields with denser vegetation (62). Ironically, the originally discovered wetland complexes (Fig. 3) (26) are less visible in DEM, although these are identifiable from ground-verified Geoeye imagery and the LMI. Here, in dense tropical marsh, the canals are visible spectrally although topographic expression is low, ∼0.5 m. North of the Rio Bravo’s riparian forest (Fig. 3), canals and the only Maya mounds yet discovered within these wetland systems were visible from pedestrian survey but not with DEM or LMI.

Fig. 3.

BOP field areas (Left, A–C) with the Maya site of Akab Muclil and Center of Gran Cacao. Left has color enhancement for elevation, and Right is a shaded relief map of the DEM.

The DEM imagery also indicated a small wetland field complex of ∼0.3 km2 at Chawak But’o’ob (Fig. 1) (49) and the largest complex of ∼7.7 km2 in the middle of the Rio Bravo valley (Fig. 1). At this latter field complex, ground-verified canals lie in close proximity to Maya architecture. These field complexes range from ∼33.5 masl to ∼13 masl in the upper valley. Most wetland fields and canals are clearly visible from both aerial photography and the lidar DEM at Chan Cahal (Fig. 1), a perched wetland fed by perennial springs, and at Chan Cahal East in the coastal plain floodplain (50). These wetland complexes together are ∼1 km2 and range from ∼17.6 to 22 masl and ∼8 to 12 masl, respectively. At Chan Cahal North, from ∼17.6 to 20 masl, the DEM did not reveal most of existing canals because of thick pasture grasses and low elevation difference of ∼0.5 m between canals and fields. The LMI again clearly shows the canals and fields, but as the difference in soil moisture in the canals versus the fields reflected in the visible and infrared signatures of the vegetation (Fig. 1 and ). Therefore, the LMI combination of intensity of the first returns of each lidar pulse distinguishes canals based on spectral characteristics of the vegetation in the canals rather than topographic differences. We also note the DEM allowed the first local identification of another form of intensive agriculture in terrace groups at Blue Creek (Fig. 1), the region’s largest intensively studied ancient Maya center between La Milpa and Lamanai.

The collective area of wetland complexes is 14.08 km2 within one medium-sized watershed with a dry season discharge of ∼6–7 m3·S−1. All of these wetland areas are floodplains or groundwater seeps subject to aggradation (64) and, thus, a large area of canals may be buried and invisible to lidar.

Field Chronology and Uses.

To understand how and when the ancient Maya used these ∼14 km2 of canal and field complexes requires evidence of timing, use, and processes of formation. Ancient wetland canal and field systems are diverse across Mesoamerica (13) but greatly understudied, and a key distinction is whether they are perennially or seasonally wet (34). Possible canals and fields in upland bajos of Guatemala and Mexico (1, 54) lie mostly within seasonally wet-dry areas. The BOP fields, however, are perennial with the water table near the surface during the wet season and within the capillary fringe in the dry season. Previous studies found that BOP wetland complexes formed through ancient Maya canal excavation and field raising to manage water for water quality and quantity. Deposition from human-induced erosion, natural flooding, and natural gypsum precipitation also buried the wetland systems (26, 27, 64).

Evidence for the BOP Late/Terminal Classic wetland agroecosystems came from excavations of more than 17 trenches from 2005 to 2018 across fields and canals, which we mapped with aerial and ground survey (26, 27) after first identifying the canals from a small airplane in 2004. The BOP wetland field complex preserves its stratigraphic records because of sediment accumulation in the Rio Bravo floodplain that buried the wetland fields and canals (Fig. 2). The surrounding BOP floodplain has been aggrading from 0.82 to 1.92 mm·yr−1 at nearby sites over the last 3,484–3,718 y (Fig. 2), but, locally, fields experienced faster rates in the Late/Terminal Classic period, accelerated by ancient Maya field raising from canal excavation (49). Gypsum precipitation, from groundwater with high concentrations of dissolved Ca+2 and SO2-4, added to aggradation (65).

Since lidar imagery indicated five times more wetland field area than what was previously estimated from pedestrian and aerial photo surveys, we studied areas A and C (Fig. 3) with six new trenches and 15 additional AMS (accelerator mass spectrometry) ages (Fig. 2 and ). Based on these, the stratigraphy of the clay-rich floodplain soils of the wetland fields provides evidence of ancient Maya field building (Fig. 4). Here, O and A horizons lie over Cy (gypsum rich), Cg (gleyed), and Ab (buried topsoil) horizons above earlier Cg horizons. The buried soils occur in field zones, between the canals, and lie below widespread ash layers from burn events preceding field and canal building. The burn layers are largely continuous, except where canal digging erased them, and 1 to 3 cm thick with pale colors of their major minerals (CaCO3, CaSO4, and silicate clays) where burning of the dark, organic matter occurred. AMS ages date to the Late to Terminal Classic within three separate ash layers (1309–1083 BP), sandwiched by three Terminal to Postclassic (1170-800 BP) dates 12 to 20 cm above the ash and two Late to Terminal Classic (1300–1060 BP) dates within 15 cm below the ash (). One AMS age, which lies 180 cm below the surface and 10 cm below a burn layer, dates to the Late Preclassic, which is the earliest date for field burning (Fig. 4 and ). We interpret the zones 20 cm below and above the ash layers as buried topsoils, Ab horizons, composed of the preburn paleosol, the ash layer, and sediment added above the ash (Figs. 2 and 4). The depths of this whole zone vary from 72 to 180 cm below the surface (). These field zones also had ceramics, faunal remains, greater melanism, charcoal, and geochemical evidence (26). Above the Ab horizon, most AMS ages become younger in the aggrading floodplain over the last millennium, whereas below the Ab, most radiocarbon ages become older.

Fig. 4.

Photo of the BOP 3 field excavation profile and photos of the fields and canals excavation profiles for BOP 11 and 16 (Upper) and the model of formation for these wetland complexes (Lower) (see for all AMS ages).

Canal sections (Figs. 2 and 4) provide evidence of their last uses and deposition after disuse. The canal surfaces are today 20 to 50 cm lower than field surfaces but were ∼100 cm lower when they functioned, having filled faster after abandonment due to higher trap efficiency. Many canals still flow in the wet season. The canal soils usually have O horizons or A horizons built atop layered organics, clays, and silt or sand layers with precipitated gypsum. Canal zone excavations often reveal layers of running and slack water deposits, lying above mélanges of the last Maya uses mixed with subsoils (Cg horizons left behind from soil profile truncation by ancient Maya canal digging). A series of radiocarbon samples across this boundary shows Terminal to Postclassic ages through the depositional canal sediments (Figs. 2 and 4). The buried fields and canal bottoms produced the only artifacts (ceramics and lithics); the few diagnostic ones date to the Late or Terminal Classic. Some wood artifacts appear in the Postclassic-aged sediments above the canal bottom zone, and the nearby ancient Maya centers of Akab Muclil and Gran Cacao have some Postclassic artifacts (52).

Other lines of evidence for the uses of these field complexes includes pollen, phytolith, charcoal, macrobotanical, and faunal. Previous investigations from Maya period regional wetlands indicate multiple ancient food species such as maize (Zea mays), avocado (Persea americana) and other fruits, squash (Cucurbita spp.), cassava (Manihot esculenta), and arrowroot (Maranta arundinacea) as well as Gossypium spp (26, 27, 65–67). Multiproxy sediment cores obtained within 1,000 m of the fields also produced Classic through Postclassic maize pollen (52). Faunal remains also indicate widespread protein harvesting. Burned animal bones and shells, especially the mollusk genera Pachychilus and Pomacea, were common in excavations. Another line of evidence for how the Maya used these field systems are carbon isotope ratios (δ13C) from dated profiles of the soil’s humin fraction that provide geochemical evidence of vegetation because plants using the C3 photosynthetic pathway dominate the area today while the main C4 plants here are maize and weedy species associated with human disturbance. Synthesized findings from four canals and two field zones show the δ13C rising from the C3 plant signatures (∼−28‰) at the surface to as much as 64% derived from C4 plants (∼−12‰) in the Late Classic sediments, which also have maize pollen and macrobotanical evidence, before returning back to C3 signatures in older, lower sediments (26–28, 49).

The watershed’s three other wetland field zones near the nonelite centers at Chan Cahal, Akab Muclil, and Chawak But’o’ob provide nearby comparisons (Fig. 1). At Chan Cahal, near the mouth of the watershed, phytoliths and macrobotanicals from a Late Classic structure (∼1200 BP) indicated squash (Cucurbita spp.), maize cobs (Zea mays), beans (Phaseolus spp.), cacao (Theobroma cacao), and sapodilla (Manilkara zapota) fruit (50). Wetland agriculture may have started in the Late Preclassic based on maize pollen from this time (2140–1870 BP (26, 50)) and persisted into the Postclassic (after 1000 BP) at Chan Cahal (50) and Akab Muclil (52), but most agricultural evidence from the canals and fields dates to the Late and Terminal Classic (26–28, 50). At Chawak But’o’ob, in the upper end of the watershed, evidence for wetland farming in the Late and Terminal Classic included the surface expression and stratigraphy of canals and fields, copious charcoal, and carbon isotopic evidence of C4 species. Here, the δ13C from soil humin rose through the Classic sediments by 6‰, an increase of ∼37% in C4 plants in a region today virtually without C4 plants (49).

Taken together across the watershed, these multiple lines of evidence, coupled with the lidar imagery, indicate large-scale, ancient, polycultural wetland agricultural systems, especially in the Late, Terminal, and Early Postclassic.

Conclusions and Discussion

We confirmed our hypothesis that lidar data could identify more wetland field complexes under canopy, estimating ∼14 km2 in a small watershed. Moreover, our original dating from before lidar holds well for the new areas that we excavated and dated based on lidar (). Most of the BOP and other three wetland field zones have dating consistent with the Late and Terminal Classic with some Postclassic evidence (26–28, 50). We presented evidence of Late Preclassic burning in the BOP complex, which chronologically coincides with evidence at Chan Cahal (50). For the Rio Bravo watershed, at least, we show that wetland agroecosystems were extensive and persistent, which counters the doubts that had caused research to languish in the 1990s. Our findings parallel evidence from 40 km farther north in the Belize coastal plain that these systems lasted late in ancient Maya history (41).

This large area of intensive farming has several implications. First, the role of wetland farming in Maya subsistence requires greater attention, especially since other studies may show large areas of canals and fields (1, 54). This article on northwestern Belize provides strong evidence for chronology and crops and suggests integration with populations and trade routes, but we need these lines of evidence from many more field complexes across the Maya world. Second, the wetland fields expand in the Late Classic, a time of regional urban population increase (48, 68, 69). Along slopes at this time, the Maya were also building other forms of intensification—agricultural terraces (28). The large size of the BOP and surrounding fields may reflect local use or longer-range trade. The fields sprawl around two ancient Maya sites including the extensive but little-studied Gran Cacao (70) and other nearby centers (52). The BOP fields had access to both human tumplines and canoe transport through a complex network of rivers and canals. By canoe, the BOP system linked to the Rio Hondo, Chetumal Bay, and the Caribbean, known trading routes of agricultural products (66). This large-scale agricultural expansion may have been a response to a demand for food production or perhaps for commodities and specialty items. The BOP wetlands produced several crops as well as mollusks, but we have yet to find cacao, a high-status crop that would confer some possible elite connections to the wetland fields. Excavations at Chan Cahal, 10 km north of BOP, however, recovered evidence of cacao and more than 200 pieces of jade, although of lesser quality than at the elite sites (66).

Third, evidence suggests most field and canal construction was coincident with two major environmental changes: water table rise associated with sea level rise (26, 27, 64) over the Maya Preclassic and Classic periods and droughts during the Late/Terminal Classic and Early Postclassic (46) (Fig. 2). Especially during the Maya Classic population expansions, water tables inundated these field systems, to which the Maya adapted by building wetland field complexes. Later, the Terminal and Postclassic were times of widespread population abandonment in the central Maya lowlands (69, 71) but the wetland fields persist, perhaps providing subsistence oases as uplands far above the water table became riskier. We note that most canals start to infill with no dredging in the Terminal Classic and Early Postclassic, but some show continued use while the Precolumbian centers of Akab Muclil and Gran Cacao along the margins of BOP did have Postclassic evidence (72). We also note that the large Maya center of Lamanai, just 25 km east, had a long Postclassic occupation (73). The wetland field environment, during regionally dry conditions, may have reattracted or maintained less intensive land uses like hunting, fishing, gathering, and milpa farming.

Lastly, the vast agricultural infrastructure of Maya civilization that lidar indicates in this study and possibly in others (1, 2) adds to the growing evidence for earlier, more intensive, and more wide-ranging anthropogenic impacts on tropical forests (42, 74) and wetlands (20). We also hypothesize that these activities added atmospheric CO2 from burning to build the wetland fields and CH4 from wetland creation through canal and wetland expansion (29). Studies (20, 75), for example, mapped the evolution of paddy rice agriculture in China, which coincides with the rise of 70 ppb of CH4 over the last 5,000 y. The largest premodern increase of CH4 (∼60 ppb) from 2000 to 1000 BP coincides with the rise of Maya wetland agriculture and likely a much larger area of neotropical wetland canals and fields (24, 76). The discernible area from the lidar imagery of this study is ∼14 km2 in one drainage basin, and the area is probably larger because modern plowing, aggradation, and draining have masked other canals. Other recent studies have identified additional, possible wetland field complexes in northern Belize (77), Guatemala’s Petén (1, 63), and Mexico’s Yucatan, which may be much larger than those in this study (54). However, multiple proxy research will need to confirm their areal extents, quantify their greenhouse gases, and identify their uses over time to explore their possible global roles in the Early Anthropocene.

Materials and Methods

Lidar Data Collection, Processing, and Analysis.

For this research, NCALM collected Lidar data between July second and fourth of 2016 for an area covering ∼275 km2 (Fig. 1) employing the Teledyne Optech Titan MW multispectral lidar (53). The collection was performed from a height of 570 m above the ground with each laser channel firing at 175 kHz (total 525 kHz) and the scanner oscillating at ±30° and 25 Hz. The processed point cloud yielded 11.065 billion returns from 6.516 billion laser pulses (roughly 1.7 returns per pulse), on average 23.7 pulses/m2 and 40.1 returns/m2. A vertical accuracy assessment based on 2,238 checkpoints collected via kinematic GPS along a stretch of open road within the project area yielded a RMSE of 4.4 cm. The point cloud was classified and a variety of rasters were generated based on the point cloud including half-meter resolution first surface (digital surface model [DSM]) and bare-earth (digital elevation model [DEM]) models as well as multispectral Lidar intensity images.

The lidar sensor, a Teledyne Optech Titan MW (S/N 14SEN/CON340), operates three lidar channels with different laser wavelengths and look angles. Channel 2 resembles a traditional single-channel nadir-looking lidar system, and it is based on a 1,064-nm laser (near infrared) scanned with an oscillating mirror on an arc of ±30° from nadir. Channel 1 is based on a 1,550-nm laser (near infrared), and it is scanned through the same mirror and scan angle of channel 2, but it is pointed 3.5° forward of nadir. Channel 3 operates a 532-nm laser (green) scanned with the same mirror and angle but oriented 7° forward of nadir. Each channel can operate at pulse repetition frequencies (PRF) of 50–300 kHz, for a total combined (and synchronized) PRF of 150–900 kHz (53). For this project, two GPS stations were temporarily installed while the aircraft was flying, one station was located in the outskirts of Belize City (17.513388563° N, 88.224689048° W −1.419 m) and a second one was located at the Maya Research Program camp near Blue Creek (17.888648859° N, 88.923225625° W, 157.322 m).

Lidar flight planning was based on a nominal laser shot density of 24 laser pulses per m2 with a configuration that ensures high penetration through the forest canopy based on experiments performed by NCALM in similar environments. The nominal flight parameters consist of a flying height of 570 m above ground level (AGL) and a ground speed of 70 m/s and a lateral swath overlap of 50% (edge of swath overlaps the center of the adjacent swath). The lidar instrument was configured with a PRF of 175 kHz per channel for a total combined PRF of 525 kHz. The scanning mirror was set at full motion of ±30° at a frequency of 25 Hz. While the instrument configuration is usually kept constant during a survey, the flying parameters vary due to topographic relief and winds conditions, and this will cause local and global variations on the actual laser pulse density with respect to the nominal plan.

Once the point cloud has been cleaned and classified different kinds of rasters are generated at 50-cm node spacing based on both the elevation and the spectral (intensity) data of the lidar returns. The elevation rasters include the first surface or DSM and the bare-earth or DEM according to the American usage of the term (the equivalent European usage is digital terrain model or DTM). The spectral or intensity rasters are derived from the first returns of each lidar channel (1,550, 1,064, and 532 nm) and are four different rasters, three individual gray scale images for each channel and a fourth raster where the individual intensity images are combined into a three-band false color raster (red 1,550 nm, green 1,064 nm, and blue 532 nm).

The point cloud interpolation source files are gridded tile by tile using Surfer by Golden Software, which is run in batch processing mode using a scripter. Both DEM and DSM tiles are generated by interpolating the source point clouds using the Kriging algorithm. The only parameter value difference used to generate DSM and DEM is the Kriging search radius. For dense DSM source data, the search radius is set to 5 m, while for the sparse ground return source data used to build the DEM, the search radius is set to 20 m. For the production of the intensity rasters, three different sets of source files were created with three attributes: the X and Y coordinates and the return’s lidar intensity for each first return, but separated per tile and per channel. The intensity rasters are interpolated using the inverse distance weighting algorithm based on a quadratic power and a search radius of 3 m. The individual surfer grid tiles are then mosaicked into larger raster datasets using the “Grid Mosaic” tool in Surfer, while at the same time they are exported into the ArcGIS raster float format (*.flt). Further information on NCALM lidar data processing workflows and considerations can be found in Fernandez-Diaz et al. (78).

Standard shaded relief images (azimuth 315°, elevation 45°, Z factor 1) were generated for the DSM and DEM in ArcGIS using the Hillshade tool in the Spatial Analyst Toolbox. Many of the wetland agricultural features were directly visible in the shaded relief images; in some other areas they were not as easily identifiable. Advanced lidar data visualization techniques were tried employing the Relief Visualization Toolbox (RVT) developed by the Research Centre of the Slovenian Academy of Sciences and Arts. We also used Simple Local Relief Models (SLRM) to enhance visualization of the agricultural features.

The geodetic elevations for the lidar products are referred to the WGS84 ellipsoid. Because a few geoid models exist for the study area, to convert the relative elevations of the produced DEM to elevation above mean sea level, we selected the Earth Gravitational Model 2008–WGS84 version (EGM2008) (79). The EGM2008 provides geoid undulation values (N) relative to the WGS84 ellipsoid at 2.5-min resolution for the globe. These values can be used to convert the ellipsoidal height (h) to orthometric height (H) in meters above mean sea level with the simple formula H = h − N. First, the EGM2008 was extracted for just the area covering the lidar swath and resampled to 0.5-m resolution to match the DEM. With both datasets in GCS WGS84, we applied the formula to produce a DEM in meters above mean sea level. For estimating the elevation ranges of the fields, we also removed extreme values using the “Fill” tool in ArcMap 12.2 to fill the sinks in the data.

Field and Canal Metrics.

We delineated wetland canals and field extents on a geographic information system using a combination of visualization techniques of the lidar data (). First, we overlaid a hillshade model at 80% transparency on top of the DEM with a color palette, which provides elevation data, but with the realistic topographic texture from the hillshade. Although this visualization is very useful, we knew from ground surveys that not all of the canals were visible with this technique. Next, we derived a SLRM using the RVT (80). The SLRM accentuates convex versus concave features in the landscape, more clearly revealing many of the canals that were not clear in the hillshade and DEM combination, such as at Chawak But’o’ob. We tested other products from the RVT commonly used in archaeological identification such as sky-view factor and PCA (62, 63), but these proved less useful than SLRM for identifying canals. In areas with tall (2 + m), dense grasses such as some of the savanna, the DEM is less well-defined and the surface appears rough because of the difficulty of separating vegetation from ground returns. Thus, canals under this type of vegetation are not visible in the DEM-derived visualizations. Lastly, the intensity data from the different wavelength lasers were viewed individually as grayscale images, and together as an RGB composite. The red and near infrared wavelengths most clearly defined ancient Maya canals that were not apparent in any of the other visualization methods, such as at Chan Cahal North and Chan Cahal East. Although the intensity images clearly show these canals that are not visible in the DEM, they only reveal those that are not under dense canopy but in the open modern agricultural fields and savanna. The different data products (topography and spectral intensity) are, therefore, complementary under different environmental conditions. The multiple visualization methods allowed for a comprehensive inventory of canal systems that multiple analysts were able to digitize. The linear vectors were buffered to 3 m wide to reflect average canal width based on field measurements and the lidar-derived DEM. These vectors provide the total canal length for each of the wetland field complexes. The vector product was then converted to a raster in order to classify distinct fields and canals. Lastly, the raster was converted back to a vector in order to calculate individual field dimensions and summary statistics for the systems.

Dating.

We report all 42 radiocarbon ages, including 15 new and 27 previously reported AMS radiocarbon analyses from samples of burned and not burned organic material in discrete strata collected from our field and canal excavation units in the BOP complex (Fig. 1 and ). BETA Analytics provided ages for 18 of these samples, International Chemical Analysis for 13 samples, and The University of Arizona’s AMS Lab provided ages for 11 more samples. Conventional ages are in years BP (BP = Before Present, 1950 AD), and the AMS laboratories calibrated ages using INTACAL13 and corrected all ages for fractionation. Here we report all ages using 2-sigma calibration (95% probability).

Carbon Isotopes.

Geoarchaeological studies in Central America have most frequently measured δ13C in the humin fraction of the soil organic matter, as it is the most stable part of the soil carbon (81, 82). Because of the stability, this fraction is the most likely to represent the vegetation of the past and has more potential to reveal the C4 signature of ancient maize cultivation (82). The humin fraction of the soil was extracted using a process used in previous studies of the region (81) in the Soils and Geoarchaeology Laboratory at the University of Texas at Austin. The process begins by weighing 2 g of sample homogenized to pass through a 100-mesh sieve (149 μm) into a 100-mL test tube and adding 1 M HCl until effervescence stops. The tubes are then placed in a 70 °C water bath during the reaction to assure calcium and magnesium carbonates are removed. Next, the samples are transferred to 50-mL Oakridge Centrifuge tubes and centrifuged at 10,000 × g for 30 min and the supernatant discarded. The samples are then twice rinsed overnight with deionized water and centrifuged. We then use alkaline pyrophosphate extraction (0.1 M sodium hydroxide and 0.1 M sodium pyrophosphate) to remove the humic and fulvic acid fractions from the samples. About 25 mL of the solution is mixed into the centrifuge tube with the sample, capped with a septum cap, and the headspace is flushed with nitrogen gas to remove oxygen and prevent oxidation of the humic acids. The samples are shaken in the alkaline pyrophosphate solution on a mechanical shaker overnight and centrifuged at 30,000 × g for 2 h. This step of adding alkaline pyrophosphate, shaking, and centrifuging is repeated 3 times total, discarding the supernatant each time. We then rinse the samples with water, with 0.05 M phosphoric acid to remove the added alkalinity, and then one more time with water. Finally, the samples are oven dried at 105 °C and the pellets are crushed with a mortar and pestle and are packed into tin capsules. We analyzed the samples using a Leco Elemental Analyzer and Isotope Ratio Mass Spectrometer in the Jackson School of Geosciences at the University of Texas at Austin. The samples are combusted in the analyzer and converted to CO2, and the various masses of CO2 are measured to determine the ratio of 13C to 12C, as well as the total organic carbon. The results are reported as δ13C in ‰ as compared to the standard PDB (Pee Dee Belemnite).