Bayesian analyses question the role of climate in Chulmun demography

We investigate the relationship between climatic and demographic events in Korea during the Chulmun period (10,000–3,500 cal. BP) by analyzing paleoenvironmental proxies and 14C dates. We focus on testing whether a cooling climate, and its potential negative impact on millet productivity around the mid 5th-millennium cal. BP, triggered the population decline suggested by the archaeological record. We employ a Bayesian approach that estimates the temporal relationship between climatic events and change-points in the rate of growth in human population as inferred from radiocarbon time frequency data. Our results do not support the climate-induced population decline hypothesis for three reasons. First, our Bayesian analyses suggest that the cooling event occurred after the start of the population decline inferred from the radiocarbon time–frequency record. Second, we did not find evidence showing a significant reduction of millet-associated dates occurring during the cooling climate. Third, we detected different magnitudes of decline in the radiocarbon time–frequency data in the inland and coastal regions, indicating that the even if cooling episodes were ultimately responsible of these population ‘busts’, their impact was most likely distinct between these regions. We discuss our results highlighting the long tradition of mobility-based subsistence strategy in coastal regions as a potential factor contributing to the regional differences we were able to detect.


Extended Dataset Description
We gathered a total of 894 14 C dates from samples collected from Chulmun sites in South Korea. We built our 14 C date database using the previous compilation of the Chulmun dates made available to the scholarly communities by Chulmun specialists 1,2 .
To ensure the scientific reproducibility and validity of the data used for analysis, we further applied the following filters to the collected dates. First, we eliminated samples of which we could not clearly identify its published primary source (i.e. site excavation reports). Second, we eliminated dates which were deemed too early or late to establish a reliable cultural association to the Chulmun period. More specifically we removed samples whose uncalibrated 14 C age were earlier than 6,200 BP or later than 2,800 BP. Third, we eliminated shell and animal bone dates, which contain potential marine signature, to control for the marine reservoir effect potentially impacting the result of our analysis. Finally for multi-component sites spanning multiple chronological periods, we only selected the dates with clear contextual association with a Chulmun feature to eliminate chronological ambiguity in dates. After these filtering, the total number of 14 C dates used for our analysis was 682. Of these, 640 (93.8%) were wooden charcoal, 29 (4.3%) various annual grains, and 13 (1.9%) uncharred wood. We further distinguished the dates that have contextual association with millet. We establish millet association to a date when the date itself was from a millet grain or was found in the same stratigraphic layer with confirmed presence of millet either by archaeobotanical or pottery impression evidence.
We then divide the 14 C datasets into two spatial groups: inland and coastal. We deemed a date to be coastal, if it was recovered from a site located within 2 km of the current coastal line. Otherwise, the date was placed into the inland group. To account for the sea-level change that occurred since the Chumun period, an exception was made if an inland site is known to have been located on the coast during the Chulmun occupations. The sites applicable to this exception were Bibongri, Sugari, Gadong, and Jukrimdong shell mounds sites in the Southeastern region of Korea 3,4,5,6 . These now-inland sites yielded artifacts and ecofacts of clear marine origins such as marine mammals, fish, and shells, thus the dates recovered from these sites are regarded as coastal.

Data Preparation for Summed Probability Distribution and Bayesian Analysis
We generated summed probability distributions of calibrated radiocarbon dates (hereafter SPD) of the coastal and inland regions between 7,000 and 3,000 cal. BP, using the rcarbon package 7 (v.1.4.1) of the R statistical computing language 8 . To account for potential biases arising from idiosyncratic definitions of archaeological sites (which often do not necessarily equate to past settlements) and inter-site variations in the sampling intensity of radiocarbon dates, we first grouped sites into spatial clusters using the DBSCAN algorithm 9 , using a distance threshold ε equal to 1km and setting the minimum number of points per cluster to 1, and subsequently generated temporal bins 10 using the binPrep function in rcarbon and setting a temporal distance (h) to 100 years. We thus define each spatio-temporal bin to be a group of radiocarbon dates that are close to each other in space and time. Both the visual display of the regional SPDs ( Fig. 1, main text) and the mark-permutation test for assessing the difference in their shapes (Fig. S1, see below) were carried out by using as a basic unit of analyses these spatio-temporal bins. In both cases, radiocarbon dates from the same bin were calibrated using the IntCal20 calibration curve 11 , and their resulting vector of probabilities summed and divided to the number of dates in the same bin (i.e. normalised to sum to unity).
The Bayesian model fitting procedure had a slightly different workflow. Firstly, to account for potential edge effects we excluded from our sample all dates with a calibrated cumulative probability within the time range of analyses (i.e. 7,000 -3,000 cal. BP) below 0.5. Then we selected from each spatio-temporal bin the radiocarbon date with the smallest 14 C error (randomly selecting one in case multiple dates had the same 14 C error). Table S1 shows the number radiocarbon dates, bins, site clusters, and the sample used for Bayesian analyses for each region.

Coastal vs Inland Site Dates
Temporal density of 14 C dates associated with coastal and inland sites were compared using mark permutation test 12 with 1,000 Monte-Carlo iterations to compute the local and global P-values. Results (Fig. S1) show significant departures from the null hypothesis of the two SPDs being samples from the same statistical population, with a Global P-value < 0.001.

Figure S1
Mark permutation test between coastal and inland site.

Millet Dates
To assess the relationship between foxtail (Setaria Italica) and broomcorn millet (Panicum miliaceum) cultivation on the one hand and demography on the other, we extracted all radiocarbon dates associated with either crops from both regions (Fig. S2), generated an SPD (Fig. S3) and compared its shape to an SPD based on all radiocarbon dates via a mark permutation test 12 as implemented in the rcarbon R package 7 . The technique consists of comparing the observed SPD (in this case of millet dates, Fig. S1) against an envelope of simulated SPDs generated by randomly shuffling the marks associated with each date (in this case millet vs non-millet dates). Thus in practice, in this case the procedure is equivalent to generating SPDs by randomly sampling 45 radiocarbon dates from the pool of all radiocarbon dates available in the window of analysis. The null hypothesis is that there is no difference in the shape of the SPD based on millets vs the SPD based on all dates, which is equivalent to a stationary proportion of dates associated with the crop. SPD in this case were generated without spatio-temporal binning and the global P-value obtained via 1,000 permutations. The results (Fig. 4, main text) yielded a p-value of 0.69, suggesting that we cannot determine whether the rise and fall in the millet SPD is different from the overall fluctuations in the density of radiocarbon dates, and hence possibly human population dynamics.

Figure S2
Calibrates dates associated with foxtail and broomcorn millets from coastal and inland sites (n=45).

Figure S3
SPD of radiocarbon dates associated with either foxtail or broomcorn millet from coastal and inland regions. Vertical bars represent the median calibrated date of each specimen.

Age-Depth Modelling
We used a compound Poisson-Gamma model as implemented in the Bchron R package 13 to re-fit an agedepth model for the SSDP-102 14 (Fig. S4), GY-1 15 (Fig. S5), and Pomaeho 16 (Fig.S6) sediment cores using the recently published IntCal20 and Marine20 calibration curves 11,17 . For each dataset we used 50,000 MCMC iterations with 10,000 burn-in steps to infer the posterior dates of each layer of the sediments. The Poisson-Gamma model provides also an estimate of outlier probability for each radiocarbon date 13 , and thus we ran the model iteratively excluding all dates with an outlier probability above 0.1 at each iteration.  Table S4. Radiocarbon dates for the GY-1 sediments. Figure S4. Age-depth model of the SSDP-102 sediment core.  For each dataset we have identified the start and end point of the largest cooling event within the temporal window of analysis (Fig. S7), namely between layer 1050cm and 1035cm in SSDP-102 (d1 and d2), between 1020 and 1017cm at Pomaeho (f1 and f2), and between 890 and 887 cm at GY-1 (g1 and g2), and extracted posterior samples (Fig. S8) to determine their chronology. Figure S7. Climatic proxies and cooling events: alkenone-derived sea surface temperature reconstruction from the SSDP-102 sediment core (top) 14 , ratio of arboreal pollen to total pollen from Pomaeho site (middle) 16 , ratio of arboreal pollen to total pollen from the GY-1 sediment core (bottom) 15 Figure S8. Posterior distribution of the start and end point of abrupt cooling events between 7,000 and 3,000 at SSDP-102 (d1 and d2), Pomaeho (f1 and f2), and GY-1 (g1 and g2).

Bayesian Modelling of Growth Rates and Change-points
We used the nimble 21,22 and nimbleCarbon 23 R packages to fit the following Bayesian model: where is the calendar date of the sample , ( ) is the corresponding 14 C age based on the IntCal20 calibration curve, is the root of the sum of the squares of the 's 14 C age error and the corresponding error in the calibration curve, and is the observed 14 C age of . The BoundedDoubleExponential model is generalised bernoulli distribution where the probability of observing a sample from the calendar year (in BP) is given by the probability mass function: and with = 0 when > and when < . Here we are interested in estimating from the BoundedDoubleExponential the growth rates 1 and 2 , and their change point . To achieve this we used the following constants and priors: = 7000 = 3000 1~( = 0, = 0.0004) 2~( = 0, = 0.0004) ~( = 3000, = 7000, = 5000, = 1000) which ensures a wide range of realistic patterns (see Fig. S9) with growth rates comparable to those observed from other SPDs and a weakly informative prior for reducing the probability of the shift in growth rate at the extremes of the window of analysis. Posterior samples for the coastal, inland, and combined datasets were obtained using nimble's Metropolis-Hastings adaptive random-walk sampler, with three chains of 100,000 iterations, a burn-in of 10,000 steps, and a thinning interval of 6 iterates (see Fig. S10 for trace plots). The 90% highest posterior density interval, Gelman-Rubin convergence diagnostic (̂), and the effective sample size of the three parameters for the three models are shown on table S5, while the fitted model and marginal posterior distributions are shown on Fig. S11 and S12.    To assess the goodness-of-fit of the growth model we carried out a posterior predictive check using the postPredSPD() function in the nimbleCarbon R package. The function generates an envelope of SPDs according to the fitted model and visualises time-intervals where the observed SPD showed a higher or lower density of radiocarbon dates. We generated the fitted model envelope of the three datasets using a 95% percentile interval computed from 500 simulated SPDs using randomly drawn parameter combinations from the joint posterior distribution (Fig. S13). Figure S13. Posterior predictive checks of the fitted models