2.1 High-resolution data relating to the 774 and 993 CE cosmic events

Miyake et al. (Reference Miyake, Nagaya, Masuda and Nakamura2012, Reference Miyake, Masuda and Nakamura2013) discovered two large and rapid increases in atmospheric ¹⁴C concentration (Δ¹⁴C) from Japanese tree rings dated to 1176/5 cal BP (774/5 CE) and 957/6 cal BP (993/4 CE). Corresponding anomalous ¹⁰Be and ³⁶Cl events recorded in ice cores (Mekhaldi et al. Reference Mekhaldi, Muscheler, Adolphi, Aldahan, Beer, McConnell, Possnert, Sigl, Svensson, Synal, Welten and Woodruff2015) suggest extreme fluxes from high-energy solar particles were responsible. The NH expression of these events has been detailed in several studies (see Büntgen et al. Reference Büntgen, Wacker, Galván, Arnold, Arseneault, Baillie, Beer, Bernabei, Bleicher and Boswijk2018 and references therein). Güttler et al. (Reference Güttler, Adolphi, Beer, Bleicher, Boswijk, Christl, Hogg, Palmer, Vockenhuber, Wacker and Wunder2015) measured Δ¹⁴C levels in NZ kauri, showing their global nature. Annual or biennial SH Δ¹⁴C data from seven tree-ring series, six relating to the 774 CE event (four from NZ and one each from Chile and Tasmania) and one to the 993 CE event (NZ)—Table 1 are available from the Büntgen et al. (2018) study (ETH/AMS facility, identifier ETH, alpha-cellulose). They calculated offsets to contemporaneous NH data of ~30 ¹⁴C yrs (in terms of Δ¹⁴C, these offsets equate to 4.0 ± 0.4‰ in 774 CE and 3.5 ± 0.7‰ in 993 CE). These data are available from the SHCal database (http://intcal.org/shcal20/— datasets 9-1, 9-2, 9-3, 9-4, 9-5, 9-6). The high-resolution SH data compare well with existing SH decadal data (Figures 1A and 1B).

Table 1 Summary of new high-resolution SH tree-ring chronologies/data sets including the time span, block interval (number of yrs [i.e. rings] per analysis) and total number of analyses (N). From Büntgen et al. (Reference Büntgen, Wacker, Galván, Arnold, Arseneault, Baillie, Beer, Bernabei, Bleicher and Boswijk2018).

^a Series numbers 9-1 to 9-6 correspond with “Set No.” and “Division” in SHCal database (http://intcal.org/shcal20).

Figure 1A Southern Hemisphere tree-ring ¹⁴C data relating to the 774 CE event. Tree-ring series include: 3-1 (Wk, Manoao colensoi, Oroko Swamp, NZ and Libocedrus bidwillii, Takapari Forest Park, NZ); 3-6 (Wk, Lagarostrobos franklinii, Stanley River, Tasmania, Australia); 9-1, 9-2, 9-3 (ETH, Agathis australis, Dargaville, NZ); 9-4 (ETH, Manoao colensoi, Moana, NZ); 9-5 (ETH, Austrocedrus chilensis, El Asiento, Chile); 9-6 (ETH, Lagarostrobos franklinii, Stanley River, Tasmania, Australia).

Figure 1B Southern Hemisphere tree-ring ¹⁴C data relating to the 993 CE event. Tree-ring series include: 3-1 (Wk, Manoao colensoi, Oroko Swamp, NZ and Libocedrus bidwillii, Takapari Forest Park, NZ); 3-6 (Wk, Lagarostrobos franklinii, Stanley River, Tasmania, Australia); 4-1 (UB, Manoao colensoi, Oroko Swamp, NZ); 6-2 (CAMS, Lagarostrobos franklinii, Stanley River, Tasmania, Australia); 9-1 (ETH, Agathis australis, Dargaville, NZ).

2.2 New 450–0 cal BP data sets

Four new data sets in the interval 450–0 cal BP have been included in SHCal20 (see Table 2). Data set 3-7 contains 22 Wk decadal solvent-extracted alpha-cellulose high precision radiometric analyses spanning the interval 1725–1935 CE and first given in McCormac et al. (Reference McCormac, Hogg, Higham, Baillie, Palmer, Xiong, Pilcher, Brown and Hoper1998a). This group of analyses was accidentally omitted from SHCal04 and SHCal13 but is included in SHCal20.

Table 2 Summary of new 450–0 cal BP SH tree-ring chronologies/data sets showing the time span, block interval (number of yrs [i.e. rings] per analysis) and total number of analyses (N). Dating methods: HPLSC—high precision liquid scintillation counting; AMS—accelerator mass spectrometry.

^ 1915–1945 decades given preliminary multiple solvent (acetone) extractions.

¶ duplicate analyses.

# 1931–1950 decades from Takapari FP (40°04’S, 175°59’E).

+ solv. extr. = solvent (acetone) extracted.

^‡ kauri derived from the upper North Island of NZ but exact geographic location unknown.

Data sets 3-8 and 3-9 are composed from Wk decadal alpha-cellulose AMS analyses spanning the interval 1705–1945 CE and published in Turney et al. (Reference Turney, Palmer, Hogg, Fogwill, Jones, Bronk Ramsey, Fenwick, Grierson, Wilmshurst, O’Donnell and Thomas2016). This publication also contained complementary data from Lake Tay (Western Australia) and Campbell Island (NZ). The Lake Tay data set was not considered for SHCal20 because of suspicions that the very high resin content characteristic of the Callitris columellaris trees may have resulted in translocation of ¹⁴C across ring boundaries. The Campbell Island Dracophyllum longifolium data, also omitted, was systematically older than contemporaneous measurements from mid-latitude SH locations, probably as a result of its high latitude Southern Ocean location. See Hogg et al. (Reference Hogg, Heaton, Bronk Ramsey, Boswijk, Palmer, Turney, Southon and Gumbley2019) for a discussion of both the Lake Tay and Campbell Island data sets. Data set 3-10 contains duplicated Wk five-ring alpha-cellulose AMS analyses spanning the interval 1652–1827 CE and published in Hogg et al. (Reference Hogg, Heaton, Bronk Ramsey, Boswijk, Palmer, Turney, Southon and Gumbley2019).

Hogg et al. (Reference Hogg, Heaton, Bronk Ramsey, Boswijk, Palmer, Turney, Southon and Gumbley2019) presented a new Bayesian spline method for calibration curve construction (see below) and tested it on 11 SH data sets. The four new data sets presented here had relative offsets from the mean of ~ +6 to –6 ¹⁴C yrs indicating their suitability for inclusion in SHCal20.

2.3 New 3608–3453 cal BP data sets

Two new single tree ring data sets in the interval 3608–3453 cal BP have been included in SHCal20 (see Table 3). Data set 7-1 contains 16 University of Arizona AMS Facility analyses (identifier AA) spanning the interval 3608–3590 cal BP and 65 AA analyses spanning 3520–3453 cal BP. Pearson et al. (Reference Pearson, Wacker, Bayliss, Brown, Salzer, Brewer, Bollhalder, Boswijk and Hodgins2020 in this issue) compare a sub-set of these data with single year measurements (AA) from Irish oak and North American bristlecone pine, and provide additional fine scale information on the N-S offset during this period.

The two new AA data sets are plotted with SHCal13 in Figure 2.

Table 3 Summary of new 3608–3453 cal BP SH tree-ring chronologies/data sets showing the time span, block interval (number of yrs [i.e. rings] per analysis) and total number of analyses (N).

Figure 2 New single-ring ¹⁴C data sets (7-1) in the interval 3608–3453 cal BP. Series 7-1 (AA, Agathis australis, Dargaville NZ).

2.4 New Younger Dryas-age measurements

Hogg et al. (Reference Hogg, Southon, Turney, Palmer, Ramsey, Fenwick, Boswijk, Büntgen, Friedrich, Helle and Hughen2016a, Reference Hogg, Southon, Turney, Palmer, Bronk Ramsey, Fenwick, Boswijk, Friedrich, Helle, Hughen and Jones2016b) presented a floating ¹⁴C chronology containing 1022 measurements from NZ kauri trees spanning the interval ~11,285–9990 ¹⁴C yrs BP (mean 1σ uncertainty = ±28 ¹⁴C yrs). The kauri trees were located in Northland, at Towai (series name 8-1, CRW003) and Dargaville (series name 8-2, FIN11). The Towai decadal samples (~11,285–10,070 ¹⁴C BP) were derived from a well replicated and securely cross-matched tree-ring chronology (1451 rings, 91 radial strips derived from 37 trees with an average cross-correlation coefficient between all series of 0.71, Palmer et al. Reference Palmer, Turney, Cook, Fenwick, Thomas, Helle, Jones, Clement, Hogg, Southon and Bronk Ramsey2016). The 778 analyses were determined principally by University of California, Irvine (identifier UCI, holo-cellulose, AMS), Wk (alpha-cellulose, high precision liquid scintillation spectroscopy) and Oxford University (identifier OxA, alpha-cellulose, AMS) with a few additional analyses by ETH (alpha-cellulose, AMS) and University of Heidelberg (identifier Hd, holo-cellulose, Gas Proportional Counting - GPC).

The Dargaville decadal samples (~10,305–9990 ¹⁴C yrs BP) were obtained from a single tree (FIN11, two measured radii, estimated 533 rings), with the 244 analyses undertaken by UCI (holo-cellulose, AMS) and Wk (alpha-cellulose, predominantly high precision liquid scintillation spectroscopy).

The Towai and Dargaville ¹⁴C data sets were wiggle-matched against contemporaneous NH data contained within IntCal13, being careful to avoid aberrant sections of the curve (see Hogg et al. Reference Hogg, Southon, Turney, Palmer, Ramsey, Fenwick, Boswijk, Büntgen, Friedrich, Helle and Hughen2016a for details), returning calendar age estimates of ~13,134–11,695 cal BP (Towai) and ~11,870–11,367 cal BP (Dargaville). More recent unpublished higher resolution NH pine (Pinus sylvestris) and NZ kauri measurements by ETH have indicated that the Towai kauri YD ¹⁴C data set should be shifted towards older ages by 10 ± 3 cal yrs resulting in a final calendar age range of ~13,144–11,705 ± 3 cal BP (Sookdeo et al. Reference Sookdeo, Kromer, Buentgen, Friedrich, Friedrich, Helle, Pauly, Nievergelt, Reinig, Treydte, Synal and Wacker2019 in this issue, 2020). This shift also affects the Dargaville series, with an adjusted calendar age range of ~11,880–11,377 ± 5 cal BP. Both Towai and Dargaville data sets have been included in SHCal20 on the basis of a recent decision made by the IntCal Working Group to accept wiggle-matched data from pre-Holocene floating trees.

The four YD-age Huon pine samples analysed by Hua et al. (Reference Hua, Barbetti, Fink, Kaiser, Friedrich, Kromer, Levchenko, Zoppi, Smith and Bertuch2009) were compared with the Towai kauri series by Hogg et al. (Reference Hogg, Southon, Turney, Palmer, Ramsey, Fenwick, Boswijk, Büntgen, Friedrich, Helle and Hughen2016a). The four trees could not be cross-matched by dendrochronology alone and required additional linking by ¹⁴C wiggle-matching (Hua et al., Reference Hua, Barbetti, Fink, Kaiser, Friedrich, Kromer, Levchenko, Zoppi, Smith and Bertuch2009); of the four, only SRT-779 agreed closely with the well-replicated Towai chronology (n = 24; Acomb = 126.4%; An = 14.4%. ${\chi ^2}$-test: df = 23 T = 15.018 (5% 35.173) - see Hogg et al. Reference Hogg, Southon, Turney, Palmer, Ramsey, Fenwick, Boswijk, Büntgen, Friedrich, Helle and Hughen2016a for model details) and is retained for SHCal20. For more information on the chronological issues associated with some of the YD-age Huon pine ¹⁴C series and ultimately their exclusion from SHCal20, please refer to Hogg et al. (Reference Hogg, Southon, Turney, Palmer, Ramsey, Fenwick, Boswijk, Büntgen, Friedrich, Helle and Hughen2016a, Reference Hogg, Southon, Turney, Palmer, Bronk Ramsey, Fenwick, Boswijk, Friedrich, Helle, Hughen and Jones2016b).

The YD-age data sets included in SHCal20 are shown graphically in Figure 3. The new YD measurements refine the shape of the SH atmospheric curve from ~13,100–11,400 cal BP.

Figure 3 Southern Hemisphere YD-age tree-ring ¹⁴C data sets. Tree-ring series include: 8-1 (UCI, Wk and OxA, Agathis australis, Towai, NZ; 8-2 (UCI and Wk, Agathis australis, Dargaville NZ); 5-1 (OZ, Lagarostrobos franklinii, Stanley River, Tasmania, Australia).

Hogg et al. (Reference Hogg, Southon, Turney, Palmer, Ramsey, Fenwick, Boswijk, Büntgen, Friedrich, Helle and Hughen2016a) showed that the ¹⁴C calibration curve IntCal13, and as a consequence SHCal13, was clearly too young from ~12,200–11,900 cal BP. This time interval had previously been occupied by an incorrectly linked larch tree-ring series from Ollon, Switzerland (Ollon505, also termed VOD505, Hogg et al. Reference Hogg, Turney, Palmer, Southon, Kromer, Bronk Ramsey, Boswijk, Fenwick, Noronha, Staff and Friedrich2013c). They also showed that the YD-B chronology was incorrectly linked by dendrochronology to NH Holocene data and this caused the ~10,400 BP ¹⁴C plateau in IntCal13 and SHCal13 to be ~5 decades too short, thus shifting the two ¹⁴C peaks seen at around 12,500 cal BP towards older calendar ages by about 50 cal yrs. This finding has recently been confirmed with new dendrochronological and ¹⁴C measurements on NH trees (Sookdeo et al. Reference Sookdeo, Kromer, Buentgen, Friedrich, Friedrich, Helle, Pauly, Nievergelt, Reinig, Treydte, Synal and Wacker2019, Reference Sookdeo, Kromer, Adolphi, Beer, Brehm, Büntgen, Christl, Eglinton, Friedrich, Guidobaldi, Helle, Muscheler, Nievergelt, Pauly, Reinig, Tegel, Treydte, Turney, Synal and Wacker2020).

3.1 Creating a curve where we have SH observations

We commence by creating a preliminary SH calibration curve based only on the available SH data using the same Bayesian spline with errors-in-variables methodology of IntCal20; see Heaton et al. (Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue) for further details. For the three distinct time periods covered by the observed SH trees dated by dendrochronology (approximately 2140–0, 3520–3453, 3608–3590 cal BP) we select the same knots for our spline, both in terms of number and location, as used in IntCal20. This includes the specific, additional knots placed around 774/5 CE and 993/4 CE used to enable representation of the rapid Miyake “cosmic” events (Miyake et al. Reference Miyake, Nagaya, Masuda and Nakamura2012, Reference Miyake, Masuda and Nakamura2013). This matched knot selection with the NH IntCal20 curve enables equitable comparisons on the N-S hemispheric offset to be made and, since the density of the SH data is lower than the density of the data from the NH, is sufficient for us to capture the detail within the available SH ¹⁴C observations. For the floating SH tree rings (approximately 13,140–11,375 cal BP) we use the same number of knots as in the corresponding section of the IntCal20 curve but space them evenly over the interval. Such a placement allows the Bayesian method to explore the exact calendar location of these floating trees equitably—an uneven knot spacing introduces bias as the method will tend to align the trees so that wiggles in the underlying levels of ¹⁴C fit with higher knot density (Heaton et al. Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue). Outside these time periods, where no SH data are present, we reduce the number of knots in our spline.

Repeat observations on same-sample material are combined as for IntCal20. We also place the same prior on the additive over-dispersion, obtained from the Sixth International Radiocarbon Intercomparison (SIRI, Scott et al. Reference Scott, Cook and Naysmith2017), as used for IntCal20 (see Heaton et al. Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue, for details on how this prior was constructed). This over-dispersion attempts to identify and quantify any potential additional variability in ¹⁴C determinations beyond the laboratory-quoted uncertainty (e.g. growing season, intra-hemispheric regional effects, species) to ensure we do not obtain an over-precise calibration curve. Our SIRI-based prior, together with the posterior we obtain from the SHCal20 data, are shown in Figure 5.

Figure 5 Histogram of the posterior for $\tau $, the level of over-dispersion in the SHCal20 data. We model ${F_i}$, the observed value of F¹⁴C in (annual) datum $i$ with calendar age ${\theta _i}$, as ${F_i} = f\left( {{\theta _i}} \right) + ϵ{_i} + {\eta _i}$. Here $f\left( {{\theta _i}} \right)$ is the underlying atmospheric F¹⁴C level at ${\theta _i}$ cal BP in the SH; $ϵ{_i} \sim N\left( {0,\sigma _i^2} \right)$ the laboratory reported uncertainty; and ${\eta _i} \sim N\left( {0,{\tau ^2}f\left( {{\theta _i}} \right)} \right)$ the over-dispersion used to model any potential further additional variability seen in the SHCal20 data. Note that, as for IntCal20, our model for this over-dispersion/additional variability scales with $\surd f\left( {{\theta _i}} \right)$, the square root of F¹⁴C. We also show the SIRI-based prior in red. As expected, this SIRI-based prior on the level of over-dispersion is larger than our posterior due to the screening criteria for SHCal20 data. The posterior estimate is dominated by the SHCal20 data rather than the prior.

The MCMC Bayesian spline curve fitting process was run for 50,000 iterations with the first 25,000 discarded as “burn-in” (for an introduction to Bayesian analysis and MCMC, see Gelman et al. Reference Gelman, Carlin, Stern and Rubin2013). Due to the predominantly Gibbs update steps, and after a visual assessment of various individual parameters, this length of run was felt to be sufficient to achieve convergence. The output was then thinned, saving only every 10th iteration, to reduce auto-correlation in the fitted location of the floating trees.

The described approach provides curve estimates in these four time intervals based only upon the SH data in such a way that still allows us to share information on ¹⁴C variability over time (as represented by the spline’s smoothing parameter) between the disjoint intervals and jointly estimate the level of over-dispersion. After thinning, we obtain from this MCMC 2500 posterior part-realizations (so called since they only cover these four disjoint time intervals) of the SH atmosphere (Figure 4B) together with corresponding posterior estimates for the over-dispersion in ¹⁴C within the SH. We call these sections of the curve SH-data-based.

Outside of these four time intervals however, the above approach does not work since we do not have SH data on which to base an estimate. To create a complete SH curve from 0–55 cal kBP we need to rely upon importing information from the NH curve based upon a model for the N-S offset. However, we wish to do this in such a way that ensures that, in these four specific intervals, our final SHCal20 curve still corresponds to these SH-data-based estimates.

3.2 Extending the SH curve by modeling

We model the evolution of the N-S offset over time as a first-order autoregressive (AR(1)) process. This time-varying process models ${X_t}$, the N-S offset at calendar age $t$, as linearly dependent upon the size of the offset in the previous year, with the addition of a stochastic/random component. Specifically,

${X_t} = μ + \phi \left( {{X_{t - 1}} - μ } \right) + {ϵ_t},$

where the $ϵ{_t}$ are considered independent and identically distributed from a $N\left( {0,\sigma _e^2} \right)$ distribution. Here $μ $ denotes the mean N-S offset and the parameter $\phi $ encapsulates the N-S offset’s dependence over time. We aim to learn about this offset based upon both our SH-data-based curve and our corresponding IntCal20 curve. We can then extrapolate the NH curve accordingly in the time periods where SH data are lacking.

Each of the 2500 SH-data-based part-realizations is paired with a randomly drawn NH curve realization from the IntCal20 posterior (Figure 4C). Unlike the SH part-realizations, each of these NH realizations extend from 0–55 cal kBP. In the SH-data sections, the offset between the NH and SH realization is calculated and an AR(1) process is fitted to these observed values using maximum likelihood. We then form, conditional on the observed offset in the SH-data sections, predictions for the N-S offset for the calendar ages outside the range of the SH-data covering the complete 0–55 cal kBP period (Figure 4D). To create a full SH-realization this offset is then added to the NH realization (Figure 4E).

This approach provides 2500 complete SH-realizations (Figure 4F) which, in the SH-data-based regions, are identical to their corresponding part SH-realization. Outside they transition smoothly to an offset version of their paired NH realization, with the size of that offset realization dependent and determined according to our best, maximum likelihood, estimate.

3.3 Adding in over-dispersion and creating predictive intervals

Finally we add back in our SH-data-based measure of over-dispersion to these 2500 complete SH-realizations to create predictive intervals for the complete SHCal20 curve. This over-dispersion aims to assess any potential additional ¹⁴C variability, seen within observed SH tree-ring determinations with identical calendar ages, beyond that quoted by the laboratories. Such additional ¹⁴C variability is to be expected since laboratories are only able to quantify those elements within their control and could occur as a consequence of intra-hemispheric regional differences in ¹⁴C levels; between-tree, species or growing season differences; or inter-laboratory differences. By incorporating this additional variability into predictive intervals on the calibration curve we aim to ensure we do not provide overly precise calendar age estimates. For more details see Heaton et al. (Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue).

We note that the level of over-dispersion estimated from the SHCal20 data (shown in Figure 5) is smaller than that within the NH (Heaton et al. Reference Heaton, Blaauw, Blackwell, Bronk Ramsey, Reimer and Scott2020 in this issue: Figure 5). This could be due to the smaller range of potential sources of additional variability within the collated SHCal20 tree rings when compared to the NH IntCal20 data. The SH ¹⁴C tree-ring measurements generally have fewer overlapping trees in any time period, come from a smaller number of laboratories, cover a smaller number of tree species, and arise from a smaller set of locations. However, alternatively it might be that there is less additional variation in ¹⁴C levels within the SH, perhaps due to reduced regional effects. As a consequence of this reduced over-dispersion with the SH tree rings, as we extend past ~40 cal kBP the published SHCal20 curve has slightly narrower predictive intervals than IntCal20.

3.4 Model robustness: An alternative approach with sampling importance resampling (SIR)

The above approach of randomly pairing NH and SH curve realizations does not fully incorporate all the information one could extract from the offset model, as it does not take into account possible dependence between the uncertainties in the two curves. Specifically, we may end up with pairings whereby the calculated offset does not agree well with an AR(1) process in which case we have the potential to overestimate the uncertainty in the predicted offset and consequently the final SHCal20 curve.

To investigate if this was a genuine concern, we therefore also considered an alternative, more statistically rigorous sampling importance resampling (SIR) approach, whereby each SH part-realization was compared to all of the 2500 NH realizations stored for IntCal20. The maximised likelihood of the observed offset between the SH part-realization and each of these NH realizations was calculated under the AR(1) model to provide an importance weight for each match of the SH part-realization to all the NH realizations. The specific NH realization with which to pair the SH part-realization was then randomly chosen according to these importance weights. The offset was then extrapolated as for the random pairing but now with a pair that best approximates an AR(1) process. As for the random pairing, this SIR approach ensures that the complete SHCal20 realizations agree with their corresponding part-realization.

This SIR approach was however extremely slow to run since, unlike the random pairing, it required calculation of the maximised likelihood for all possible combinations of SH and NH realizations. Further, and more significantly, this SIR approach gave a very small effective sample size with the same few NH realizations being chosen repeatedly. In the older pre-Holocene time periods where the underlying NH data have considerable calendar age uncertainty and a large posterior sample is required to accurately represent the curve, the SIR approach consequently failed. The value of the NH calibration curve in these older time periods is almost independent of its value in the calendar periods (e.g. 2000–0 cal BP) where the weighting was performed, and the offset observed. As a result of this independence, during these older periods where we have no direct SH data, the SH curve should simply be a constant offset from IntCal20 since for these older times we are reliant upon the NH curve and our offset model to provide the SH estimate. This was not the case for the SIR based estimate. The random pairing approach did however create an SH curve estimate with a constant offset from IntCal20 as we extend back beyond our observed SH data and towards 55 cal kBP indicating that random pairing worked appropriately for these older time periods.

Further, in regions where the underlying NH curve is based on trees dated by dendrochronology and a small posterior sample is sufficient to reliably summarize the curve, the SIR approach was indistinguishable from the random pairing method. This suggests the random pairing approach is robust and does not introduce additional uncertainty, providing confidence in the reconstruction.

Finally, we note that it would be possible to construct both the NH and SH curves simultaneously, combining data from both hemispheres and a model for the offset into a joint curve creation. This is a consideration for future curve updates although it would raise the issue of potential asymmetry (which hemisphere should be considered offset?) and that, even in regions where SH data were dense, the SH curve would be influenced by the NH data. Such an integrated MCMC approach may also re-introduce some of the mixing difficulties seen in the random walk approach of IntCal13 making curve estimation much more computationally intensive.

3.5 Value of the N-S offset

Each of the 2500 pairings of NH and SH curve realizations provides a slightly different estimate of the N-S offset. This will result in us estimating 2500 different AR(1) processes each with its own mean and standard deviations. In Figure 6, we present the estimated mean offsets and standard deviation over time for all these 2500 realized offsets. Based upon these 2500 paired curves we estimate the mean N-S offset to be 36 ¹⁴C yrs with a standard deviation of 27 ¹⁴C yrs (i.e. 36 ± 27 ¹⁴C yrs).

Figure 6 Histograms of the estimated means and standard deviations (sd) of the AR(1) processes used to model the N-S hemispheric offset.

For time intervals based upon IntCal20 data adjusted for the N-S hemispheric offset, we recommend that readers also consult the complementary IntCal20 paper (Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell, Ramsey, Butzin, Cheng, Edwards and Friedrich2020 in this issue) which details the differences between IntCal20 and IntCal13.

A comparison of the SHCal13 and SHCal20 curves without the underlying data for clarity is given in Supplementary Material (Figure S1).

4.1 The Inter-Tropical Convergence Zone (ITCZ) and inter-hemispheric air mass mixing

The youngest portions of IntCal20 and SHCal20 curves are based on tree rings from temperate latitudes of each hemisphere. An obvious question is which calibration curve should be used in the tropics and subtropics given no long tree-ring ¹⁴C records for these regions are currently available.

Several studies suggest that atmospheric ¹⁴C over tropical and subtropical regions could be influenced by north-south air-mass mixing as a result of monsoon circulation. Hua et al. (Reference Hua, Barbetti, Zoppi, Fink, Watanasak and Jacobsen2004a, Reference Hua, Barbetti and Zoppi2004b) and Hua and Barbetti (Reference Hua and Barbetti2007) discussed influences of SH air masses on atmospheric ¹⁴C over northern tropical and subtropical regions during the pre-bomb and post-bomb periods. Similarly, Hua et al. (Reference Hua, Barbetti, Levchenko, D’Arrigo, Buckley and Smith2012) reported monsoonal influences of NH air masses on atmospheric ¹⁴C over a southern tropical island (latitude 5°S) in Indonesia during the post-bomb period. The use of either IntCal20 or SHCal20 curve for regions which are influenced by monsoon circulation, might be therefore inappropriate, and a mixed curve accounting for north-south air-mass mixing can be used for these regions. This approach is similar to that for tropical South America proposed by Marsh et al. (Reference Marsh, Bruno, Fritz, Baker, Capriles and Hastorf2018).

Monsoon circulation can result in air-mass mixing and high precipitation (McGregor and Nieuwolt Reference McGregor and Nieuwolt1998; Zhang and Wang Reference Zhang and Wang2008), and in several monsoon regions the bands of highest summer precipitation approximately follow the position of the ITCZ (Figure 7). However, climate model simulations and reanalysis data indicate that not all areas experiencing high precipitation are strongly influenced by air masses from the opposing hemisphere. In particular, reanalysis data for the past 3 decades (Figure 7) shows that the intrusion of NH air into the southeast trending South Pacific Convergence Zone (SPCZ) is sharply truncated (but see Section 4.2). We therefore consider it more prudent to use wind data and the seasonal ITCZ positions rather than precipitation to estimate the zonal boundaries of the areas in which IntCal20, SHCal20 and a mixed curve are used.

Figure 7 World map showing means of winds (1981–2010) and precipitation (1979–2018) at 925 hPa for two different periods: December–February (top panel) and June–August (bottom panel). These data are derived from the NCEP reanalysis wind (https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.pressure.html) and the CMAP precipitation (Xie and Arkin Reference Xie and Arkin1997), respectively. Dotted and dashed lines represent the ITCZ during December–February and June–August, respectively.

We therefore recommend the use of (i) IntCal20 for areas north of the ITCZ in June–August (dashed lines in Figure 7) which receive NH air masses all year round, (ii) SHCal20 for areas south of ITCZ in December–February (dotted lines in Figure 7) which see SH air masses all year round, and (iii) a mixed curve for areas between the two seasonal ITCZ positions shown in Figure 7, which receive northern air masses in December–February and southern air masses in June–August.

The degree of north-south air-mass mixing for the mixed curve is another issue. For vegetation, which grows almost all year round, the degree of air-mass mixing should be 50%:50% (northern:southern air masses). For seasonal vegetation, the degree of air-mass mixing might be different than the above value. However, Hua and Barbetti (Reference Hua and Barbetti2007) estimated a 52 ± 13% contribution of southern air masses to northern tropical and subtropical regions influenced by monsoons during the boreal summer months, which was based on atmospheric and tree-ring ¹⁴C data during the bomb-peak period. We thus recommend a simple approach with 50%:50% of north-south air-mass mixing for the mixed curve, and age calibration using that curve can be easily carried out using the “Mix_curve” command in OxCal program (Bronk Ramsey, Reference Bronk Ramsey2009) or selected as a curve option in CALIB 8.2 (Stuiver et al. Reference Stuiver, Reimer and Reimer2020).

4.2 Past ITCZ variations

As the mean seasonal positions of the ITCZ and the strength of monsoon circulation change at different timescales (inter-annual, decadal, centennial, etc.; e.g., Haug et al. Reference Haug, Hughen, Sigman, Peterson and Röhl2001), the boundaries of the mixed-curved areas in the past can be different from those shown in Figure 7, which are based on the modern meteorological data. At least for the remote past, the use of IntCal, SHCal or a mixed curve for a particular region in the tropics and subtropics becomes less critical for many studies when the maximum ¹⁴C offsets of these curves, which are the N-S offsets of 36 ± 27 ¹⁴C yrs (based upon our estimates in Section 3.5) are negligible compared to the ages of the events to be dated.

However, this may not be the case for researchers studying human expansion over short time scales. For example, Goodwin et al. (Reference Goodwin, Browning and Anderson2014) have reconstructed Pacific climate and in particular wind field patterns associated with changing El Niño and La Niña events, over the last ~1000 yrs, using a general circulation climate model extensively calibrated with paleoclimate data. Their model showed large variations on decadal and longer timescales in the mixing of NH air into the South Pacific via the SPCZ, dependent on the state of the Pacific Decadal Oscillation and the decadal-scale predominance of El Niño vs La Niña events, with minimal mixing in some bi-decades and a mixed region extending as far as 30°S across most of the Pacific in others.

It seems highly likely that similarly detailed wind field reconstructions in other areas will also show significant variations in the past boundaries of the zone covered by the movement of the ITCZ. Until new tree-ring ¹⁴C data sets are created from within the tropics and subtropics, only approximations for the boundaries of the mixed region and the regional balance of northern and southern air masses within it are possible, and researchers wishing to calibrate radiocarbon data from these regions should be mindful of these limitations. Future research efforts should focus on ¹⁴C measurement on tree rings from these areas in order to construct proper calibration curves for these regions, and to gain a better understanding of regional ¹⁴C offsets and their temporal variations.