An integrative approach to estimating productivity in past societies using Seshat: Global History Databank

Abstract

This article reports the results of a collaborative effort to estimate agricultural productivities in past societies using Seshat: Global History Databank. We focus on 30 Natural Geographic Areas (NGAs) distributed over 10 major world regions (Europe, Africa, Southwest Asia, South Asia, Southeast Asia, East Asia, Central Eurasia, North America, South America, and Oceania). The conceptual framework that we use to obtain these estimates combines the influences of the production technologies (and how they change with time), climate change, and effects of artificial selection into a Relative Yield Coefficient, indicating how agricultural productivity changed over time in each NGA between the Neolithic and the 20th century. We then use estimates of historical yield in each NGA to translate the Relative Yield Coefficient into an Estimated Yield (tonnes per hectare per year) trajectory. We tested the proposed methodology in two ways. For eight NGAs, in which we had more than one historical yield estimate, we used the earliest estimate to anchor the trajectory and compared the ensuing trajectory to the remaining estimates. We also compared the end points of the estimated NGA trajectories to the earliest (the 1960s decade) FAO data on crop productivities in the modern countries encompassing Seshat NGAs. We discuss the benefits of this methodology over previous efforts to estimate agricultural productivities in world history.

Keywords

agricultural revolution agricultural technologies artificial selection crop productivity historical estimates methodology Seshat: Global History Databank

Introduction

Agricultural productivity and its variation across time and space plays a fundamental role in many theories of human social and cultural evolution (Currie et al., 2015). Because agricultural societies are characterized by higher population densities, compared to hunter-gatherer societies, the transition from foraging to agriculture fueled a series of population expansions world-wide, which spread the farmers’ culture, language, and genes with them (Bellwood, 2005; Diamond and Bellwood, 2003; Renfrew, 1992; Reich, 2018; Shennan, 2018). Adoption of a farming way of life is also central to many theories about how large-scale complex societies evolved over the past 10,000 years. According to these theories, agricultural systems allowed for “surplus” production, which enabled an extensive division of labor and the appearance of artisans, artists, specialized managers, and governing elites (Johnson and Earle, 2000; but see also Meller et al., 2018).

However, although the Neolithic Revolution was one of the most consequential transitions in human evolution, we lack systematically collected quantitative estimates of the productivity of past agricultural systems on a large enough scale. In particular, by how much did the adoption of agriculture increase human carrying capacity, the maximum population density a particular environment can sustain indefinitely, without depleting critical resources needed for survival and reproduction (for a review, see Cohen, 1995)? Even more importantly, how did agricultural productivity/carrying capacity among agrarian populations vary over time (as new technologies were invented and spread, and as a result of climate change) and in space (between geographic regions characterized by different soils, landscapes, and climate)? It is this latter question that we begin to address here.

The relationships between crop yields, weather, and climate have been the focus of a great deal of attention in the Earth system science literature. This is due to concerns about securing food supplies for our growing populations and the potential challenges that climate change poses (Oyebamiji et al., 2015). Most studies have been concerned with establishing the current relationships between climate and crop yields, or making projections about changes in crop yields due to future climate change rather than extending this approach back into the past (Tóth et al., 2012). Where historical information is used, it tends to be on a relatively recent time scale. For example, Ramankutty and Foley (1999) project geographic distribution of permanent croplands on a global scale from 1992 back to 1700 (but don’t offer estimates of cropland productivity).

More recently, researchers have attempted to infer the location and intensity of agricultural production during the Holocene on a global scale (Goldewijk, 2005; Kaplan et al., 2011). These approaches are ultimately derived from estimates of past population sizes (e.g. McEvedy and Jones, 1978) and make assumptions about how human populations use land for agriculture. Although such studies should be applauded for their ambitious scale, they have a number of features that make them less-than-ideal for our purposes. First, in order to test certain theories it is desirable to separate out achieved production and realized population density from potential production and carrying capacity (i.e. the population that could theoretically be supported). A number of interesting hypotheses about human social and political evolution invoke “population pressure” as a key variable in causing changes in human societies (i.e. how close actual population is to potential population, and the stresses induced when there is competition for land and resources). For example, demographic-structural theory (Goldstone, 1991; Turchin and Nefedov, 2009), argues that state instability and societal collapse is a result of the pressures on resources from population growth, which, in turn, leads to population decline. Boserupian models of agricultural change see agricultural innovations themselves as resulting from population pressure (Boserup, 1966, 1981). Second, this approach does not make full use of the historical and archaeological information about past agricultural systems that could potentially inform estimates of productivity. Finally, the data on past populations are often quite imprecise (and sometimes just guesswork) and are typically made at the coarse-grain level of a province or whole country (Boyle et al., 2011). Additionally, realized population densities can be a very poor guide for crop productivities, because they are affected by many other factors than agriculture. For example, the provision of security from violence by large-scale territorial states allows a greater proportion of land to be cultivated and can increase population densities by a factor of 3–4 (Bennett, 2015, 2020; Turchin, 2003).

More localized studies have achieved a greater degree of success with incorporating historical or archaeological information into crop productivity estimates. For example, Nicholas (1989) developed assessment of changes in agricultural productivity in pre-contact Oaxaca, Mexico, based on ecological features of the landscape and changes in productivity of maize due to artificial selection. Goodchild (2013) used GIS techniques to model the distribution of agricultural productivity in the Roman Republican period using historical information about yields of important crops in different environments. The approach taken in this paper attempts to merge these two kinds of approaches.

For a comprehensive empirical test of various theories of social and cultural evolution (see Turchin et al., 2015) we need, ideally, such estimates spanning all major world regions and time since the invention of farming. In a previous publication (Currie et al., 2015) we outlined an empirically informed framework for deriving such estimates for a set of past societies that together constitute the Seshat sample of NGAs (the World Sample-30; (Turchin et al., 2018)). Here we report the results of this multi-annual multi-investigator research project. We emphasize, however, that the estimates that we have produced so far should be thought of only as initial approximations, which need to be refined in future research by us and other investigators. Our goal in presenting these quite preliminary results is to serve two important purposes; firstly to accelerate the ultimate development of better models through review, criticism, and collaboration with the present framework, and secondly to demonstrate which observed effects can already be reproduced with a minimal model.

This article is organized as follows. The Methods section presents a conceptual framework for estimating agricultural productivity in the past and provides information on various components that go into the estimates. The Results and Discussion section describes these estimates, and the diversity of intensification trajectories followed by societies in the Seshat World Sample-30, spanning all major world regions and going back to the origins of agriculture. The Conclusion, once again, emphasizes the limitations of the current estimates, highlights the gaps in our knowledge of past agriculture production practices, and outlines a research agenda for addressing these issues.

Methods

The general idea underlying the proposed integrative approach is to take account of multiple processes and factors that can affect crop yields, and then document how they change with time. The first set of factors we need to know about is the production technologies: the focal cultivar, the cropping system (most importantly, what proportion of arable land is under cultivation at any given time), and whether fertilization and/or irrigation is practiced. Second, crop yields are affected not only by production practices but also by cumulative effects of artificial selection that produced higher-yielding crop varieties over time. Third, crop yields fluctuate in response to climate, which has changed over the past millennia. Putting together the temporal sequence of crop yield improvements resulting from better technologies and estimated climate influences gives us a time-series of relative changes in overall crop yields. To translate this Relative Yield Coefficient (RYC) into absolute yields expressed in (metric) tonnes per hectare, we need to “anchor” the time-series in a historical estimate of crop yields for the focus area. The translation of relative yields into absolute ones is an important step, because we are interested not only in when and how much productivity increased in a particular spatial location, but also in comparing productivities across space (between different geographic locations). Further, absolute yields are key to translate productivity measures into estimates of potential carrying capacity. Finally, we use additional historical yield estimates, if available, and FAO (Food and Agriculture Organization of the United Nations) data on modern yields as a check of the method (see Figure 4 below for a flow chart summarizing the inputs and outputs of the procedure). The rest of the Methods section deals with each of these factors in turn and then explains how they were put together to estimate yields.

Production technologies

We took data on production technologies from the Seshat: Global History Databank, an interdisciplinary project which brings together data on archaeological and historical variables in order to allow comparative analysis of human social evolution (François et al., 2016; Turchin et al., 2015). The agricultural data in the databank currently features codified information on the past 10,000 years for thirty bounded geographic regions, termed Natural Geographic Areas (NGAs), selected to provide global geographic coverage as well as encompassing a range of historical societies, from small-scale societies to chiefdoms, states, and empires. The databank also contains coded data relating to more than 1000 variables, spanning themes including social complexity, agriculture, warfare, rituals, and institutions. Here we are concerned with the agricultural section of Seshat.

At the most basic level, the database records the main carbohydrate sources used by the population for subsistence. Our estimating procedure focuses on the most important – focal – crop variety, or staple. These focal crops include wheat, rice, maize, yams, sweet potatoes, and breadfruit.

The next kind of information we need is on the proportion of potentially arable land that is under cultivation at any given moment in time, which provides the basis for estimating the cropping coefficient. In continuously cultivated systems (no fallow) the cropping coefficient is 1. A two-field system, such as that broadly practiced in medieval Europe, in which crops are grown in one field, while the second lies fallow, has the cropping coefficient of 0.5. It is possible for the cropping coefficient to be greater than 1 under the conditions of multi-cropping, in which land produces more than one crop per year. Conversely, under the conditions of swidden (also known as shifting, or slash-and-burn cultivation) only a small proportion of arable area is utilized in any given year. Where data on the length of cultivation and fallow are available, we use them to estimate the cropping coefficient. Otherwise, we use a generic estimate, obtained as follows. A review of 330 studies of shifting cultivation by Merz (2002) found that most fallow periods fall in the interval between 5 and 30 years, with the modal period of 10 years. Assuming that each plot is cultivated for 3 years and is left fallow for the modal 10 years implies a generic cropping coefficient for swidden (proportion of arable land under cultivation) equal to 3/(10 + 3) = 0.23. This estimate falls in the middle of estimated fallowing coefficients used in the FAO-IIASA study of Global Agro-Ecological Zones (Tóth et al., 2012).

Two other agricultural technologies can have a large effect on crop productivities: irrigation and fertilization. The Seshat Databank records information about the appearance and implementation of these practices in each NGA. However, in order to quantify how these management practices improve crop yields, it is necessary to look to modern agricultural studies. Experimental studies from the field of agricultural science provide quantitative measures of how crop yields can be improved through the addition of fertilizer or water (Table 1).

Table 1.

Generic fertilization and irrigation coefficients used in this study. Coefficients calculated from agricultural experiments which compared treatment (fertilizing, irrigation) to a control that used traditional agricultural techniques.

Crop	Practice	Impact on yields	Reference
Wheat	Application of animal manure as fertilizer	65% improvement	Zhao et al. (2009)
Wheat	Irrigation	48% improvement	Chapagain and Good (2015)
Rice	Application of animal manure (bovine)	30% improvement (judged by tiller size)	Sudarsono et al. (2014)
Rice	Irrigation	30% improvement	Aggarwal et al. (2008)
Maize	Application of animal manure as fertilizer	78% improvement	Celik et al. (2010)
Maize	Irrigation	16% improvement	Meng et al. (2013)

Artificial selection

Wheat

Over the roughly 12,000 years since wheat, barley, and similar cereals were domesticated in the Fertile Crescent, the agricultural yields of these grains changed dramatically (Araus et al., 2007). Grain yields can be analyzed in terms of three primary yield components: number of plants per unit of cultivated area, number of kernels (grains) per plant, and mean kernel weight (Moragues et al., 2006). By measuring the dimensions of fossilized grains, archaeologists are able to reconstruct the changes over time in the last component, kernel weight. Fortunately, kernel weight is also one of the best predictors of the overall yield per unit area (Araus et al., 2014; Moragues et al., 2006), as shown in Figure 1. Although the relationship between kernel weight and yield is slightly curvilinear, in the historically relevant range of [0–30 mg] a linear approximation works as well, and we use the linear relationship for simplicity.

Figure 1.

(a) Relationship between kernel weight and cereal yield for northern Mesopotamia and western Mediterranean (data sources: Araus et al., 2014; Moragues et al., 2006). Linear regression (solid line): Y = 0.038 W, R² = 0.62, where Y is yield in t/ha/year and W is kernel weight in mg. Curvilinear regression (dashed curve): Y = 0.016 W + 0.0008 W², R² = 0.74. (b) Change in the kernel weight with time since domestication, c. 10,000 BCE (data from Araus et al., 2007, 2014; Ferrio et al., 2006). “Wheat” includes several Triticum species (einkorn, emmer, and naked wheat). The stepped line traces out the relationship between kernel weight and time that we used in estimating past wheat productivities.

Figure 1b shows that kernel weight increased after domestication (c. 10,000 BCE) until roughly 6000 BCE. Subsequently, between 6000 BCE and 1500 CE, kernel weights showed no evidence of systematic change, fluctuating around the level of 20 mg. Finally, the post-1500 period was another time of rapid artificial selection. We approximated this evolution of the domesticated cereals by a piece-wise linear function (see Figure 1b).

Maize

A highly informative archaeological proxy for past maize yields is the mean corn cob length (Kirkby, 1973). Figure 2a shows that this relationship is curvilinear. Kirkby (1973) also estimated how the mean corn cob length changed since domestication. We fitted this relationship with a second-degree polynomial (Figure 2b).

Figure 2.

(a) Relationship between corn cob length and maize agricultural yield for Oaxaca, Mexico (data source: Kirkby, 1973). Regression: Y = -0.03 L + 0.013 L², R² = 0.80, where Y is maize yield in t/ha/year and L is corn cob length (cm). (b) Change in the mean corn cob length over time (Kirkby, 1973). Regression: L = 21.2 + 0.0031 (T–5000) + 9.9 × 10⁻⁸ (T − 5000)², R² = 0.99, where T is year BCE/CE.

We estimated the effect of artificial selection on maize yields by combining the two fitted relationships. First we estimated the mean corn cob length (L) for each century after 5000 BCE, and then substituted the estimated L into the regression model for Y to calculate the artificial selection coefficient.

Rice

As with wheat and maize, we use rice grain size as a proxy for the effects of artificial selection. Broadly speaking, archaeobotanists recognize two main rice sub-species, Japonica (short/medium grain) and Indica (long grain), which reflect different histories of domestication and hybridization. Grain breadth in both species has a positive relationship with time, with similar slopes. Here we focus on Japonica, because the data on grain size that we have so far been able to locate for Indica, has poor temporal coverage. Figure 3 shows the evolution of rice grain breadth between 6000 BCE, when the earliest field systems for rice cultivation appeared in China (Fuller et al., 2014), and the present. The pattern of change is similar to what we observe for wheat: initial increase, a pre-modern plateau, and then a second increase during the modern period.

Figure 3.

Change in the mean grain breadth in Japonica rice (data sources: Fujita et al., 1984; Fuller et al., 2014; Kim et al., 2013). The stepped line traces the evolution of grain breadth over time, which we used in estimating past rice productivities.

Climate influences

We used a dynamic global vegetation and crop model, the Lund-Potsdam-Jena managed land model (LPJmL) (Bondeau et al., 2007), to estimate how agricultural productivity changed over time in different parts of the world in response to fluctuations of climate. The LPJmL was developed to predict future global vegetation patterns under possible climate change scenarios, but Oyebamiji et al. (2015) used this model to build an emulator of future crop yields, which we have simplified and reconfigured here as a palaeoemulator to look into the past. The palaeoemulator uses the statistical relationship between climate variables and crop yields to produce global maps of crop yields for 1000-year time slices going back in time from 2000 CE to 8000 BCE. The Oyebamiji et al. model can be configured to produce results for both rainfed and irrigated agricultural systems. We ran the palaeoemulator for each timestep and extracted the “low input” rainfed results for each NGA. A more detailed description of the approach and illustration with two case study NGAs (Latium and the Oaxaca Valley) are in Collins et al. (2020). Notably, the palaeoemulator does not currently distinguish between soil types (e.g. loess vs clay soils), which may have been a factor affecting cultivar choice and cropping practices among early agrarian populations (Bakels, 2013; Shennan, 2018); this is a ripe area for future expansion of the approach described here.

The crop yield emulator requires detailed information on past climate as input. We produced an estimate of climate through the Holocene using an emulator of the PLASIM-ENTS climate model (Holden et al., 2014). The paleoclimate emulator (Holden et al., 2015) uses dimensional reduction to derive global climate fields, and was constructed from an ensemble of simulations that were forced with constant atmospheric CO₂ but variable orbital forcing. Holocene climate time slices were calculated by applying the time-varying orbital configuration (Berger, 1978) to this emulator. To convert climate output from the PLASIM-ENTS paleoemulator (approximately 5°) to the high spatial resolution used by the LPJmL emulator (0.5°), the output was filtered through the ClimGen tool, which applies anomalies to observed high-resolution climatology (Osborn et al., 2016).

Historical yield estimates

We derived the historical yields used in this analysis from secondary source literature. In most literature, historical yields are estimated using historical documentation or archaeological information (Goodchild, 2013; Slicher Van Bath, 1963) or are reconstructed through the modeling of environmental factors (Bakels, 2013; Lee et al., 2006). Other methodologies include reconstruction through archaeological proxies (Kirkby, 1973) and experimental archaeology (Reynolds, 1994). The relevant literature lacked estimates only for the Finger Lakes region. Our estimate for Finger Lakes is based on data for the Cahokia region (Mt.Pleasant, 2015) but accounts for differences in climate between the two locations. Some of the underlying data were collected as part of previous research on carrying capacity by Seshat: Global History Databank team (Currie et al., 2015), expanded in the ensuing years through our ongoing work.

The estimation procedure

The procedure for estimating the yield trajectory in each NGA works as follows (see Figure 4). First, the program reads from the file the identity of the focal crop. Using Seshat data, the program constructs a sequence of cropping coefficients for each century between 10,000 BCE and 2000 CE in the NGA. For example, if Seshat data indicates that between 1245 and 1634 CE agriculture in the NGA was based on a two-field cropping system, then years 1300, 1400, and 1600 are assigned values of 0.5. For periods when agriculture was not practiced in the NGA, we assign the value of 0.

Figure 4.

Flowchart of inputs and calculations needed to estimate the agricultural yields in each NGA (Natural Geographic Area) as they change over time (at 100-year steps).

The program next constructs fertilizer, irrigation, and artificial selection coefficients in a similar way, using the values appropriate for the NGA’s focal crop (see Table 1). Because palaeoemulator data are resolved at 1000-year steps, we need to interpolate them for the century step used in our procedure (we used linear interpolation). For each century mark in the NGA sequence these coefficients are multiplied together to produce the Relative Yield Coefficient (RYC).

The final step is to scale the RYC using the historical yield data for the NGA coded in Seshat. The scaling coefficient is estimated by fitting RYC on historical yields using a linear regression with no intercept (because RYC = 0 should map onto historical yield = 0). The estimation procedure was implemented in R (version 3.6.3) and the script, as well as the data files it uses, are published as Supplemental Online Materials (https://osf.io/kjw8c/).

Comparisons between estimates and independent data

The estimation procedure requires at least one estimate of historical yield to translate RYCs into estimated yield in tons/ha. For the NGAs where we have historical yield estimates for more than one time period, we used the earliest one to scale the RYC, and the subsequent ones (not used in generating estimates) as tests of the predicted yield trajectory. An additional comparison is provided by comparing the end point of the Estimated Yield trajectory (at year 1900) to the FAO data on the modern yields in the country that encompasses the NGA.

Results and discussion

Testing the approach with historical yield estimates

As we stated in the Methods, to test the approach we used one (the earliest) historical yield estimate in an NGA to translate the Relative Yield Coefficient into absolute annual yields in t/ha. Additional yield estimates, reported by archaeologists or historians, then serve as an empirical check of the estimation procedure. Figure 5 shows trajectories for eight NGAs for which we have such multiple independent estimates. We now discuss each of these trajectories in some detail.

Figure 5.

Trajectories of estimated agricultural yields (black curves) in a selection of NGAs, compared to historical yield estimates from the literature (red symbols: downward and upward pointing triangles indicate the upper and lower ranges of reported productivities, crosses indicate the median, and the horizontal line is the extent of the temporal period to which the estimate applies). The blue diamonds are productivity estimates calculated from the FAO data for the modern countries encompassing the NGAs (an average of three periods: 1960–1970, 1995–2005, and 2010–2016).

Latium

We are lucky to have multiple estimates for this well-studied region, covering the development of Italian agriculture over the past 3000 years. The early Iron Age estimate of c. 0.5 metric tons of wheat per hectare was used to anchor the trajectory. Historical estimates from different authors for the Roman Empire range over a broad interval between 0.5 and 2 t/ha (Goodchild, 2007: Appendix V). Our estimate (1 t/ha) underestimates the median (1.5 t/ha). Conversely, the estimates of the productivity of Medieval agriculture in Italy from the literature (just under 0.5 t/ha in c.1000 CE) are below our estimates. The trajectory also somewhat overestimates Early Modern yields, and then gets on track with the data during the 19th century. The trajectory endpoint of 1.75 t/ha in 1900 is just below the average wheat productivity in Italy during the 1960s, as suggested by the FAO data. The decline in historical productivity estimates between the Roman Optimum and Early Modern Periods is probably driven, at least in part, by the worsening of the climate (Harper, 2017; McCormick et al., 2012). Unfortunately, the temporal resolution of the climatic data we used (at 1000 years) is insufficient to capture this effect.

Upper Egypt

The estimated trajectory appears to capture accurately the two transitions in the history of Egyptian agriculture: (1) the doubling of yields from the early Neolithic to the Pharaonic period and (2) another doubling during the Modern period. The 1900 end point is right on top of the average productivity of Egyptian agriculture estimated from the FAO data for 1965. The remarkable stability of Egyptian agriculture for more than 4000 years in between these two transitions is probably due to the effect of the Nile, whose annual floods provide both irrigation and fertilization (Allen, 1997).

Paris Basin

As with Italy, we are lucky to have multiple literature estimates covering the four millennia between the early Neolithic and the present. Again, our trajectory overestimates the productivity of the Medieval French agriculture, but then gets on track by the Early Modern Period. As we noted above for Latium, this decline could have resulted from the worsening climate. However, another consideration is the switch to the three-field rotation system, which took place at about this time (Bakels, 2013; Wickham, 2005). Lower wheat yields, thus, could be compensated by pulses – a hypothesis that would be worth exploring in future research. The 1900 endpoint, on the other hand, is very close to the FAO-1965 estimate, which preceded the dramatic productivity increase in modern French agriculture during the second half of the 20th century.

Susiana

The estimated trajectory is anchored by the estimate of 0.6–1.0 t/ha for winter barley during the Ur III period. The next estimate, 1.33–1.46 t/ha, is by Robert McC. Adams for “primitive agricultural technologies in the Diyala region in the 1950s” (Adams, 1965: 17). Although Susiana today is part of Iran, we compare this NGA’s trajectory to the FAO data for Iraq, since the source of our estimates is derived from Mesopotamia. The trajectory shows a very good fit with this second estimate, but we really need more historical data to test this important region more properly.

Middle Yellow River Valley

The trajectory for China is anchored in the estimate of 0.26 t/ha for millet during the early Neolithic period (note the broad temporal range indicated by the horizontal red line). The estimated trajectory is very close to the next estimate (1 t/ha for the early Han Dynasty). Historical yield estimates for the early Modern and Modern periods exhibit a great degree of variance, with our trajectory passing through the higher end. Similarly, the 1900 endpoint is slightly above the 1965 FAO estimate.

Cambodian Basin

The estimated yield trajectory is anchored by a literature estimate for the peak of Angkor Empire, 1.5 tons of rice per hectare. The next historical estimate, 2.5 t/ha in 1900 is right on the trajectory. The 1965 FAO estimate, however, is much lower, perhaps reflecting the observation that the Cambodian Basin is a high productivity area within Cambodia, due to plentiful water provided by the seasonal inundation of Tonle Sap lake.

Valley of Oaxaca

The trajectory is anchored by an estimated yield of maize agriculture, 0.5–0.6 t/ha, during the early Monte Alban I period. The subsequent continuous improvement of maize yields, due to its remarkable response to artificial selection (see Figure 2) is well-captured. However, the 1900 endpoint overpredicts the FAO 1965 estimate, and the overall productivity of Mexican agriculture catches up with the predicted level only by 2000.

Kansai

Japan’s trajectory, anchored by an estimate of 1.3 t/ha for 1100–1300 CE, does a good job for the Early Modern period, but underpredicts the rapid improvement of Japanese agriculture during the modern period.

Comparison of historical and modern yields

Overall, these eight case-studies covering all major world regions (Europe, Asia, Africa, North and South America, and Oceania) for the three major world crops (wheat, rice, maize) suggest that the procedure that we describe in this article results in reasonable estimates. Historical estimates predating the Green Revolution are particularly useful in anchoring and validating the Estimated Yield trajectories, because the timing of the onset, and of the improvement rate associated with the Green Revolution is highly variable across modern countries. Nevertheless, it is instructive to compare the most recent premodern Seshat estimate (for 1500) to the productivities of the modern countries, reported by FAO (Figure 6). This comparison was done on all NGAs.

Figure 6.

Relationship between estimated yields in Seshat NGAs in 1500 CE (calculated in this study) and the average yields in the modern country containing the NGA in 2000 CE (FAO data). Dashed line is the linear regression (R² = 0.69). Two-letter codes refer to ISO abbreviations for modern countries, with the exception of HI (Hawaii), MO (Missouri for Cahokia), NY (New York for Finger Lakes), and YK (Yakutia).

On one hand, we observe a strong correlation between Seshat-1500 and FAO-2000 values (R² = 0.70, corresponding to a correlation coefficient of 0.84). Modern yields are much higher than the Seshat estimates (on average, 3.5 times higher), but this is as expected, given the dramatic improvements in agricultural technology in the second half of the 20th century. On the other hand, there is enough scatter to indicate that modern yields could be poor indicators, even in relative terms, of pre-modern productivities in some world regions. Consider the cluster of three NGAs/countries, for which the pre-modern Seshat estimates lie well above the fitted line: Oaxaca (MX), Cambodian Basin (KH), and Sogdiana (UZ). These three contemporary countries are all relatively under-developed and over-populated. Underdevelopment means that they do not use the most productive high yielding (but also high input) agro-technologies. Overpopulation means that crops are often grown in marginal areas, which depresses average yields for the whole country. Compare, for example, Cambodian and Paris Basins. Around 1500 the productivities in both areas (supported by historical yield estimates) were c.1.5 t/ha/year (in rice and wheat, respectively). Yet by 2000 French agriculture, which is one of the most advanced in the world, generated average yields of 7 t/ha, while the average yield in Cambodia was still 2 t/ha (improving to 3 t/ha by 2015). A 20th century estimate of rice productivity for the flood-recession area is 2.5 t/ha, reflecting better suitability of such land for agriculture, compared to the rest of Cambodia. Our overall conclusion from this comparison, then, is that modern productivities can be a poor guide to the productivities in the past, even when used as relative indicators. This underscores the importance of developing independent, dynamic models of past productivities using inputs (ecological, climatological, technological, and societal) from historical societies themselves, rather than relying on mechanical projection from present circumstances to the past.

Relative strength of input factors

In this section we ask the following question: how important are the relative contributions of the various inputs that we have quantified (climate, artificial selection, irrigation, fertilizing, and cropping coefficients) to the model output, estimated yield? We addressed this question in two ways. First, we regressed the estimated yields on input factors. Table 2 (the column Estimate) reports the standardized regression coefficients associated with each factor (standardized coefficients, or beta weights, are the estimates resulting from a regression analysis in which the underlying data have been standardized so that the variances of dependent and independent variables are equal to 1.) This analysis indicates that the greatest effects are due to artificial selection and cropping coefficients. In contrast, paleoclimate and fertilizer coefficients have a much smaller effect on the estimated yields. The effect of irrigation is very close to zero.

Table 2.

Relative contributions of various factors to estimated yields. Estimate: the standardized regression coefficients. Difference (%): mean absolute deviation between the estimates using all factors and estimates omitting one factor at a time, expressed in percent of all-factor estimates.

	Estimate	Difference (%)
Paleoclimate	0.225	9.98
Artselect	0.700	214.88
Irrigation	-0.024	10.03
Fertilizer	0.160	14.79
Cropping	0.439	144.12

The second approach is to calculate the estimated yields by omitting one of the input factors and determine how much this omission affects the estimates. Table 2 (Difference %) suggests, again, that artificial selection and cropping system are the most important input factors – omitting either of these processes from estimation would result in gross errors of more than 100%. In contrast, the omission of any of the other three factors results in errors of only 10%–15%; in other words, an order of magnitude difference.

This analysis, thus, suggests that, in order to further refine yield estimates for past societies, most effort should go in obtaining better data on the effects of artificial selection and cropping systems on crop yields. However, one important caveat is that the inclusion of shorter timescale variability would be expected to alter the relative importance of different factors, and even for the long-timescale changes considered here, our estimates, in particular of climate influence, require much more work, as we discuss in the next subsection.

Comparison to other climate models

The climate emulator used here is derived from PLASIM-ENTS (Holden et al., 2014), a 3D model of the atmosphere coupled to slab ocean, slab sea-ice, and dynamic vegetation models. This climate model was used because earlier work has already addressed the PLASIM-ENTSem-ClimGEN-LPJmLem coupling (Warren et al., 2019), which converts seasonal climate outputs from an emulator of PLASIM-ENTS to the high spatial resolution and monthly variables (temperature, precipitation, cloud cover, and wet-day frequency) required by the crop emulator (Oyebamiji et al., 2015). The principal weaknesses of this particular climate model emulator relative to some more recently published paleo-climate model products (Fordham et al., 2017; Holden et al., 2019) are the simplified ocean and sea-ice dynamics and boundary conditions (we here consider only orbital forcing) and the neglect of sub-millennial variability (Fordham et al., 2017).

We validate aspects of our model in Figure 7, which compares the emulated spatial pattern of mid-Holocene precipitation in northern summer with (i) PALEO-PGEM (Holden et al., 2019), an emulator of the intermediate complexity coupled atmosphere-ocean GCM PLASIM-GENIE (Holden et al., 2018), (ii) Paleoview (Fordham et al., 2017), derived from a transient simulation of the Community Climate System Model version 3 (Liu et al., 2009), and (iii) the mean and variance of the Paleoclimate Model Inter-comparison Project PMIP2 OAV (coupled ocean-atmosphere-vegetation) ensemble (Braconnot et al., 2007). We choose the mid-Holocene (c. 4000 BCE) as it is a well-studied warm interval, with climate quite different from today due to the different orbital forcing.

Figure 7.

Comparison of June to August averaged precipitation between the climate model used in this study (Seshat) with other climate models (see the text for explanation).

As Figure 7 shows, our climate emulator captures the large-scale dynamics of the more complex models, most notably the strengthening of the Asian monsoon system and the increased precipitation in the Sahel. The spatial patterns of change lie within the uncertainty envelope (standard deviation) displayed by the PMIP2 ensemble. Given that our results are relatively insensitive to climate forcing, we expect that using an existing higher-complexity climate model would not greatly affect our predicted yield trajectories, which are dominated by artificial selection and cropping. However, this hypothesis needs to be tested directly, and we plan to do this in future work using a more complex climate model and crop emulators targeted specifically at paleoclimate states.

Conclusion: Limitations of this study and prospects for the future

Estimating agricultural productivity on a global scale over long periods of time, going back all the way to the beginnings of agriculture, is not an easy task. Here we proposed and implemented a procedure that systematically incorporates various factors that affect agricultural yields in different world regions and different time periods. Any such scheme must be based on an explicit model, which means making simplifying assumptions. In particular, we assumed that we can use a single agricultural crop, a major carbohydrate source for populations inhabiting each NGA. Previous studies focusing on maize agriculture in Mexico and wheat in Italy suggest that the resulting approximations are not too far from reality (Nicholas, 1989; Scheidel, 2001; Turchin and Nefedov, 2009). However, this assumption needs to be further tested in additional world regions.

Next, due to limitations of historical data on how various cultivation practices (fallowing, irrigation, and fertilization) affected yields in the past, our approach relies heavily on generic coefficients derived from modern studies. As we saw in Results and Discussion, implementation of these practices has a substantial effect on the resulting trajectory of estimated yields. While using generic coefficients is a reasonable first step (and certainly better than not taking these processes into account), what is needed next is an investigation that will allow us to tailor this approach to individual Seshat NGAs, with their different climates and physical environments.

Yet another limitation of the present study regards the way in which climate effects are currently handled. A millennium-long time step is far too long for many of the questions that we want to ask about the evolution of complex societies and, particularly, for the periodic breakdowns that they go through. A time step of 1 century, or, better, of a decade, would be much more informative. This is a high priority research direction for the Seshat project.

Finally, the estimation procedure used in this article calculates the Relative Yield Coefficient by combining the input coefficients multiplicatively. This approach ignores possible interaction effects between different processes that affect crop yields. For example, irrigation can ameliorate the effect of climate drying, while not being necessary as climate becomes wetter. Certain technological advances can alter the productivity of areas, for instance the introduction of the heavy plow in Europe allowed farmers to bring previously marginal clay soils into productive cultivation. This is another avenue for future investigation.

Despite these caveats, we argue that a quantitative approach that explicitly incorporates various influences shaping agricultural yields in the past represents an important step forward in enabling us to understand some of the key processes underlying the rise and fall of historical societies. A significant advantage of the proposed quantitative procedure is that it is modular. This means that different modules (processes affecting yields, such as climate, artificial selection, and cultivating practices) can be improved one by one, resulting in better overall estimates (notwithstanding the caveat above on interactions). For example, once we have higher resolution climate data, we can simply substitute our crude interpolated series with a more detailed one, which immediately generates better estimates. Finer-grain climatic data will enable researchers to investigate resilience of agricultural systems to major disturbances, a topic of great interest, but which was beyond the scope of the present article.

The R-scripts together with the data on which they operate have been made publicly available at this URL: https://osf.io/kjw8c/. We invite other investigators to join us in exploring the implications of the dynamics of agricultural productivities in the past, by bringing additional data and refining the overall estimation procedure, as well as by analyzing how fluctuations in their productive bases affected the rise and demise of past civilizations.

Footnotes

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by a John Templeton Foundation grant to the Evolution Institute, entitled “Axial-Age Religions and the Z-Curve of Human Egalitarianism,” a Tricoastal Foundation grant to the Evolution Institute, entitled “The Deep Roots of the Modern World: The Cultural Evolution of Economic Growth and Political Stability,” an ESRC Large Grant to the University of Oxford, entitled “Ritual, Community, and Conflict” (REF RES-060-25-0085), a grant from the European Union Horizon 2020 research and innovation programme (grant agreement No 644055 [ALIGNED, www.aligned-project.eu]), an European Research Council Advanced Grant to the University of Oxford, entitled “Ritual Modes: Divergent modes of ritual, social cohesion, prosociality, and conflict,” and the program “Complexity Science,” which is supported by the Austrian Research Promotion Agency FFG under grant #873927. We gratefully acknowledge the contributions of our team of research assistants, post-doctoral researchers, consultants, and experts. Additionally, we have received invaluable assistance from our collaborators. Please see the Seshat website () for a comprehensive list of private donors, partners, experts, and consultants and their respective areas of expertise.

ORCID iDs

Peter Turchin

Neil R Edwards

Supplemental material

Supplemental material for this article is available online at the Open Science Framework’s website:

References

Adams

(1965) Land Behind Baghdad. Chicago: University of Chicago Press.

Aggarwal

Hebbar

Venugopalan

, et al. (2008) Quantification of yield gaps in rain-fed rice, wheat, cotton and mustard in India. Global Theme on Agroecosystems Report no. 43. Patancheru 502 324, Andhra Pradesh, India: International Crops Research Institute for the Semi-Arid Tropics. p.36.

Allen

(1997) Agriculture and the origins of the state in ancient Egypt. Explorations in Economic History 34: 135–154.

Araus

Ferrio

Buxó

, et al. (2007) The historical perspective of dryland agriculture: Lessons learned from 10 000 years of wheat cultivation. Journal of Experimental Botany 58: 131–145.

Araus

Ferrio

Voltas

, et al. (2014) Agronomic conditions and crop evolution in ancient near east agriculture. Nature Communications 5: 3953.

Bakels

(2013) The Western European Loess Belt: Agrarian History, 5300 BC–AD 1000. Dordrecht, NL: Springer.

Bellwood

(2005) First Farmers: The Origins of Agricultural Societies. Oxford: Blackwell.

Bennett

(2015) Modeling the large-scale demographic changes of the Old World. Cliodynamics 6: 57–76.

Bennett

(2020) Retrodicting the rise, spread, and fall of large-scale states in the ld world. SocArXiv Preprint.

10.

Berger

(1978) Long term variations of caloric insolation resulting from the Earth’s orbital elements. Quaternary Research 9: 139–167.

11.

Bondeau

Smith

Zaehle

, et al. (2007) Modelling the role of agriculture for the 20th century global terrestrial carbon balance. Global Change Biology 13: 679–706.

12.

Boserup

(1966) The Conditions of Agricultural Growth: The Economics of Agrarian Change under Population Pressure. Chicago, IL: Aldine.

13.

Boserup

(1981) Population and Technological Change: A Study of Long-term Trends. Chicago: University of Chicago Press.

14.

Boyle

Gaillard

M-J

Kaplan

, et al. (2011) Modelling prehistoric land use and carbon budgets: A critical review. The Holocene 21: 715–722.

15.

Braconnot

Otto-Bliesner

Harrison

, et al. (2007) Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum – Part 1: Experiments and large-scale features. Climate of the Past 3: 261–277.

16.

Celik

Gunal

Budak

, et al. (2010) Effects of long-term organic and mineral fertilizers on bulk density and penetration resistance in semi-arid Mediterranean soil conditions. Geoderma 160: 236–243.

17.

Chapagain

Good

(2015) Yield and production gaps in rainfed wheat, barley, and canola in Alberta. Frontiers in Plant Science 6: 1–10.

18.

Cohen

(1995) How Many People Can the Earth Support?. New York, NY: W. W. Norton.

19.

Collins

Currie

Oyebamiji

, et al. (2020) Estimating agricultural productivity and carrying capacity in past societies: A method combining palaeoclimatic crop models and historical/archaeological information. SocArixv Preprint. Available at: https://osf.io/preprints/socarxiv/ky9wt.

20.

Currie

Bogaard

Cesaretti

, et al. (2015) Agricultural productivity in past societies: Toward an empirically informed model for testing cultural evolutionary hypotheses. Cliodynamics 6: 24–56.

21.

Diamond

Bellwood

(2003) Farmers and their languages: The first expansions. Science 300: 597–603.

22.

Ferrio

Alonso

Voltas

, et al. (2006) Grain weight changes over time in ancient cereal crops: Potential roles of climate and genetic improvement. Journal of Cereal Science 44: 323–332.

23.

Fordham

Saltré

Haythorne

, et al. (2017) PaleoView: A tool for generating continuous climate projections spanning the last 21 000 years at regional and global scales. Ecography 40: 1348–1358.

24.

François

Manning

Whitehouse

, et al. (2016) A macroscope for global history: Seshat global history databank, a methodological overview. Digital Humanities Quarterly 10: 4.

25.

Fujita

Coronel

Yoshida

(1984) Grain-filling characteristics of rice varities (Oryza sativa l.) differing in grain size under controlled environmental conditions. Soil Science and Plant Nutrition 30: 445–454.

26.

Fuller

Denham

Arroyo-Kalin

, et al. (2014) Convergent evolution and parallelism in plant domestication revealed by an expanding archaeological record. PNAS 111: 6147–6152.

27.

Goldewijk

(2005) Three centuries of global population growth: A spatial referenced population (density) database for 1700–2000. Population and Environment 26: 343–367.

28.

Goldstone

(1991) Revolution and Rebellion in the Early Modern World. Berkeley, CA: University of California Press.

29.

Goodchild

(2007) Modelling Roman agricultural production in the middle Tiber Valley, Central Italy. PhD Dissertation. Birmingham: Institute of Archaeology and Antiquity, School of Historical Studies, University of Birmingham.

30.

Goodchild

(2013) GIS models of Roman agricultural production. In: Bowman

Wilson

(eds.) Oxford Studies on the Roman Economy, vol. 3: The Roman Agricultural Economy. Oxford: Oxford University Press.

31.

Harper

(2017) The Fate of Rome: Climate, Disease, and the End of an Empire. Princeton: Princeton University Press.

32.

Holden

Edwards

Garthwaite

, et al. (2014) PLASIM-ENTSem v1.0: A spatio-temporal emulator of future climate change for impacts assessment. Geoscientific Model Development 7: 433–451.

33.

Holden

Edwards

Garthwaite

, et al. (2015) Emulation and interpretation of high-dimensional climate model outputs. Journal of Applied Statistics 42: 2038–2055.

34.

Holden

Edwards

Rangel

, et al. (2019) PALEO-PGEM v1.0: A statistical emulator of pliocene–pleistocene climate. Geoscientific Model Development 12: 5137–5155.

35.

Holden

Edwards

Ridgwell

, et al. (2018) Climate–carbon cycle uncertainties and the Paris Agreement. Nature Climate Change 8: 609–613.

36.

Johnson

Earle

(2000) The Evolution of Human Societies: From Foraging Group to Agrarian State. 2nd ed. Stanford, CA: Stanford University Press.

37.

Kaplan

Krumhardt

Ellis

, et al. (2011) Holocene carbon emissions as a result of anthropogenic land cover change. Holocene 21: 775–791.

38.

Kim

Ahn

S-M

Jeong

(2013) Rice (Oryza sativa L.): Seed–size comparison and cultivation in ancient Korea. Economic Botany 67: 378–386.

39.

Kirkby

AVT

(1973) The Use of Land and Water Resources in the Past and Present Valley of Oaxaca, Mexico. Ann Arbor: University of Michigan.

40.

Lee

Tuljapurkar

Vitousek

(2006) Risky business: Temporal and spatial variation in preindustrial dryland agriculture. Human Ecology 34: 739–763.

41.

Liu

Otto-Bliesner

, et al. (2009) Transient simulation of last deglaciation with a new mechanism for Bølling–Allerød warming. Science 325: 310–314.

42.

McCormick

Büntgen

Cane

, et al. (2012) Climate change during and after the Roman Empire: Reconstructing the past from scientific and historical evidence. Journal of Interdisciplinary History 43: 169–220.

43.

McEvedy

Jones

(1978) Atlas of World Population History, New York, NY: Facts on File.

44.

Meller

Gronenborn

Risch

(2018) Surplus Without the State: Political Forms in Prehistory. Halle: Tagungen des Landesmuseums Vorgeschichte.

45.

Meng

Hou

, et al. (2013) Understanding production potentials and yield gaps in intensive maize production in China. Field Crops Research 143: 91–97.

46.

Mertz

(2002) The relationship between length of fallow and crop yields in shifting cultivation: a rethinking. Agroforestry Systems 55: 149–159.

47.

Moragues

Moral

LFGd

Moralejo

, et al. (2006) Yield formation strategies of durum wheat landraces with distinct pattern of dispersal within the Mediterranean basin I: Yield components. Field Crops Research 95: 194–205.

48.

Mt.Pleasant

(2015) A new paradigm for pre-columbian agriculture in North America. Early American Studies 13: 374–412.

49.

Nicholas

(1989) Land use in prehispanic Oaxaca. In: Kowalewski

Feinman

Finsten

, et al. (eds.) Monte-Alban’s Hinterland. Part II: The Prehispanic Settlement Patterns in Tlacolula, Etla and Ocotlan in the Valley of Oaxaca, Mexico. Univ. Mich. Mus. Anthropol. Mem. No. 23., pp.449–505.

50.

Osborn

Wallace

Harris

, et al. (2016) Pattern scaling using ClimGen: Monthly-resolution future climate scenarios including changes in the variability of precipitation. Climatic Change 134: 353–369.

51.

Oyebamiji

Edwards

Holden

, et al. (2015) Emulating global climate change impacts on crop yields. Statistical Modelling 15: 499–525.

52.

Ramankutty

Foley

(1999) Estimating historical changes in global land cover: Croplands from 1700 to 1992. Global Biogeochemical Cycles 13: 997–1027.

53.

Reich

(2018) Who We Are and How We Got Here: Ancient DNA and the New Science of the Human Past. New York, NY: Pantheon.

54.

Renfrew

(1992) Archaeology, Genetics, and Lingusitic Diversity. Man 27: 445–478.

55.

Reynolds

(1994) Experimental Archaeology: A Perspective for the Future. The Reuvens Lecture 5. Stichting voor de Nederlandse Archeologie.

56.

Scheidel

(2001) Progress and problems in Roman demography. In: Scheidel

(ed.) Debating Roman Demography. Leiden, Netherlands: Brill, pp.1–81.

57.

Shennan

(2018) The First Farmers of Europe: An Evolutionary Perspective. Cambridge: Cambridge University Press.

58.

Slicher Van Bath

(1963) The Agrarian History of Western Europe. New York, NY: St. Martin’s Press.

59.

Sudarsono

Melati

Aziz

(2014) Growth and yield of organic rice with cow manure application in the first cropping season. Agrivita 36: 334.

60.

Tóth

Kozlowski

Prieler

, et al. (2012) Global Agro-Ecological Zones (GAEZ v3.0): User’s Guide.

61.

Turchin

(2003) Historical Dynamics: Why States Rise and Fall. Princeton, NJ: Princeton University Press.

62.

Turchin

Brennan

Currie

, et al. (2015) Seshat: The Global History Databank. Cliodynamics 6: 77–107.

63.

Turchin

Currie

Whitehouse

, et al. (2018) Quantitative historical analysis uncovers a single dimension of complexity that structures global variation in human social organization. Proceedings of the National Academy of Sciences of the U S A 115: e144–e151.

64.

Turchin

Nefedov

(2009) Secular Cycles. Princeton, NJ: Princeton University Press.

65.

Warren

Edwards

Babonneau

, et al. (2019) Producing policy-relevant science by enhancing robustness and model integration for the assessment of global environmental change. Environmental Modelling & Software 111: 248–258.

66.

Wickham

(2005) Framing the Early Middle Ages. Oxford: Oxford University Press.

67.

Zhao

Wang

, et al. (2009) The effects of two organic manures on soil properties and crop yields on a temperate calcareous soil under a wheat–maize cropping system. European Journal of Agronomy 31: 36–42.