Significance Understanding agricultural subsistence is vital for understanding past complex societies. Lidar data are indicating widespread ancient Maya infrastructure. Wetland agriculture was crucial to ancient cultures, but no previous study coupled lidar with multiproxy evidence to demonstrate the extent and uses of Maya wetland fields. We conducted a lidar survey around wetlands that multiple use proxies established were ancient Maya polycultural systems. Lidar indicated the Birds of Paradise (BOP) wetland field complex was five times larger than we had previously mapped and identified an even larger wetland agroecosystem. We ground-verified the BOP fields through excavations and dating, creating a study to couple these multiproxy data with lidar, thereby demonstrating widespread ancient Maya wetland agroecosystems.

Abstract 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 SI Appendix, Table S1; 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 (CO 2 and CH 4 ) 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 (SI Appendix, Table S1 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. 1A), 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. 1D), 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. 1C) 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. 1B). All of the excavations discovered exclusively ancient Maya artifacts except at Chan Cahal (Fig. 1D), where modern cultural items occurred in the plow zones of deforested pastures.

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 (SI Appendix, Table S1). 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 CO 2 from burning to build the wetland fields and CH 4 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 CH 4 over the last 5,000 y. The largest premodern increase of CH 4 (∼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 (SI Appendix, Tables S2 and S3). 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 SI Appendix, Table S1). 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 C 4 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 CO 2 , and the various masses of CO 2 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).

Acknowledgments We thank the Belize Institute of Archaeology, the Programme for Belize Rio Bravo Conservation Management Area, and the Blue Creek Community, Belize; the National Center for Airborne Laser Mapping; The Northwestern Belize LiDAR Consortium; N. Brokaw, M. Canuto, K. Cox, N. Dunning, T. Garrison, S. Houston, T. Inomata, J. McNeill, S. Ward, G. Wells, the anonymous reviewers, and our many collaborators in Belize. This work was supported by National Science Foundation Grants 0924501, 0924510, 1114947, and 1550204; National Geographic Society Committee for Research and Exploration Grants 7861-05 and 7506-03; The C.B. Smith Sr. Centennial Chair, College of Liberal Arts, and Center for Archaeological and Tropical Studies, University of Texas (UT) Austin; the Maya Research Program, UT Tyler; and the Cinco Hermanos Chair, Georgetown University.

Footnotes Author contributions: T.B. and S.L.-B. designed research; T.B., S.L.-B., S.K., T.G., F.V., J.C.F.-D., S.E., and C.D. performed research; T.B., S.L.-B., S.K., J.C.F.-D., S.E., and C.D. analyzed data; and T.B., S.L.-B., S.K., J.C.F.-D., S.E., and C.D. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1910553116/-/DCSupplemental.