Isotopic distribution of bioavailable Sr, Nd, and Pb in Chungcheongbuk-do Province, Korea

Background Mapping the distribution of bioavailable isotope ratios across landscapes serves as an efficient geochemical tool for delineating the origins and migration trajectories of humans and animals. Chungcheongbuk‑do Province in central Korea, known for its geological diversity and inland location isolated from coastal influences, provides an ideal area to study the contributions of geological and environmental factors to the isotope landscape (isoscape). This study analyzed the distribution of bioavailable Sr, Nd, and Pb isotopes in the province using plant and soil data obtained in this study and from previous works. Findings Chungcheongbuk‑do features diverse geological elements, including Precambrian basement, Paleozoic metasedimentary rocks, and Mesozoic granitoids and volcano‑sedimentary sequences. The 87 Sr/ 86 Sr ratios of bulk soil samples from 44 sites primarily range from 0.781 to 0.706, with two ratios exceeding 0.9 originating from Precambrian basement and Cretaceous granitoid areas. Fractions of soils treated with 1 M ammonium nitrate and acetic acid exhibit indistinguishable 87 Sr/ 86 Sr ratios (R 2 = 0.99, except for one point), spanning from 0.804 to 0.707. Plant 87 Sr/ 86 Sr ratios demonstrate a robust positive correlation with leachate ratios (for ammonium nitrate data, ( 87 Sr/ 86 Sr ) plant = 0.938 × ( 87 Sr/ 86 Sr) leachate + 0.045, R 2 = 0.98). The ε Nd values of bulk soils from Precambrian basement areas (–18 to –30) plot against Sm/Nd ratios around the reference line corresponding to 3 Ga, while other bulk soil samples (ε Nd = –1 to –21) align with a younger (~ 2 Ga) reference line. Plant ε Nd values, ranging from –4 to –24, exhibit a prominent positive correlation with ammonium nitrate leachates ( plant ε Nd = 0.77 × leachate ε Nd – 3.83, R 2 = 0.89). Plant samples do not show consistent variation between 87 Sr/ 86 Sr and ε Nd . The 206 Pb/ 204 Pb and 207 Pb/ 204 Pb ratios of bulk soils show a ~ 2 Ga trend, typical for Korean basement rocks. The Pb isotopic ratios of ammonium nitrate and acetic acid leachates match perfectly with each other (R 2 = 0.99). The 206 Pb/ 204 Pb, 207 Pb/ 204 Pb, and 208 Pb/ 204 Pb ratios of plant samples vary narrowly (19.2–17.9; 15.8–15.5; 39.0–38.0) and are distinctly different from those of bulk soils (24.4–17.9; 16.5–15.6; 42.5–37.9) and their leachates (23.1–17.7; 16.2–15.6; 41.0–38.1). Plant and soil data from this study and previous works were used to construct Sr and Nd isoscapes, employing interpolation models based on inverse distance weighting, simple kriging, empirical Bayesian kriging, and geology and topography‑considered empirical Bayesian kriging regression prediction. These maps await validation through analysis of additional archives. Conclusions The isotope data obtained in this study highlight a strong geological control over bioavailable Sr and Nd, in contrast to a dominant environmental influence on bioavailable Pb. The Sr and Nd isoscapes presented here are potentially valuable for addressing archaeological or forensic inquiries in their current state. Nevertheless, the maps would benefit from additional refinement with increased sample density and enhanced interpolation models.


Research background
The radiogenic isotope system is based on the spontaneous and constant transformation of parent isotopes of one element into daughter isotopes of another element through the emission of particles and radiant energy.As well explained in textbooks (Dickin 2005;Faure and Mensing 2005), radionuclides enable the dating and source tracking of terrestrial and extra-terrestrial rocks and minerals, shedding light on Earth's geological history and the origins of the solar system.These include long-lived isotopes such as 87 Rb (T 1/2 = 48.8Ga), 147 Sm (= 106 Ga), 238 U (= 4.47 Ga), 235 U (= 0.704 Ga), and 232 Th (= 14.01 Ga), which decay to 87 Sr, 143 Nd, 206 Pb, 207 Pb, and 208 Pb, respectively.The isotopic composition of radiogenic daughter elements, particularly Sr, has also been used to provenance biological materials since its application in ecology and archaeology about four decades ago (Graustein and Armstrong 1983;Ericson 1985).This is principally based on the fact that the isotopic signatures of these elements are directly conveyed from surrounding sources into biological materials through the food chain, with insignificant isotopic fractionation resulting from the small relative mass differences of the isotopes (Price et al. 2002;Bentley 2006).
The proportion of the four naturally occurring isotopes of Sr ( 88 Sr, 87 Sr, 86 Sr, and 84 Sr) in rocks and minerals depends on their age and geochemistry (i.e., Rb/Sr ratio).Consequently, Sr isotopic composition is variable across lithologies of different ages and formation histories.Older ages and high Rb/Sr ratios observed in granitic rocks lead to higher 87 Sr/ 86 Sr ratios.Conversely, younger ages, and lower Rb/Sr ratios typically found in mafic rocks, result in relatively lower 87 Sr/ 86 Sr ratios.Besides geological sources, biological strontium is influenced by environmental sources such as aeolian dust, rainfall, sea spray, and modern fertilizers.The 87 Sr/ 86 Sr ratio transferred into biological materials does not change measurably over the time scales relevant to disciplines dealing with the comparatively recent past due to the long halflife of 87 Rb.This stability has long been used for provenancing studies in fields such as archaeology (Bentley 2006;Slovak and Paytan 2011), forensics (Bartelink and Chesson 2019), food science (Voerkelius et al. 2010;Di Paola-Naranjo et al. 2011), and natality information studies (Brennan et al. 2015a, b).
Neodymium, a member of the rare earth elements (REEs), exhibits seven naturally occurring isotopes: 142 Nd, 143 Nd, 144 Nd, 145 Nd, 146 Nd, 148 Nd, and 150 Nd.Given the slow decay rate of 147 Sm and the generally low (~ 0.15) Sm/Nd ratios in rocks and minerals, variations in the isotopic composition of Nd are typically minimal.The evolution of Nd isotopes is generally discussed in relation to a chondritic uniform reservoir (CHUR) representing the bulk earth, using epsilon notation (ε Nd = [{( 143 Nd/ 144 Nd) sample -( 143 Nd/ 144 Nd) CHUR }/ ( 143 Nd/ 144 Nd) CHUR ] × 10 4 , ( 143 Nd/ 144 Nd) CHUR = 0.512638 when normalized to 146 Nd/ 144 Nd = 0.7219).Igneous rocks derived from lithophile element-depleted (i.e., high Sm/Nd) mantle sources have positive ε values, whereas rocks of the continental crust with lower Sm/Nd have negative ε values (DePaolo 1988).Like Sr isotopes, Nd isotopes have long been used to trace surface processes and the circulation of oceanic currents (Goldstein and Jacobsen 1988;Frank 2002).Neodymium isotopes, sometimes combined with Sr and Pb isotopes, are also applicable for provenancing archaeological materials such as ceramics (Li et al. 2006;Di Bonis et al. 2018;Lü et al. 2023) and glass fragments (Ganio et al. 2013).REEs are incorporated in ppb concentrations in human and animal skeletons, and Nd isotope data from bones and teeth can be used for archaeological fingerprinting (Tütken et al. 2011).However, Nd isotope data for plants are rarely reported, even for standard reference materials (SRMs).
The decay of three major U ( 238 U and 235 U) and Th ( 232 Th) isotopes give rise to long decay chains ending in stable 206 Pb, 207 Pb, and 208 Pb.The U-Th-Pb isotope system provides a reliable and widely used dating tool (Dickin 2005;Faure and Mensing 2005).Lead isotope data have long been used to discuss environmental human health (Gulson 2008), pollution (Hansmann and Köppel 2000;Komárek et al. 2008), and archaeological issues related to tooth enamel (Budd et al. 2004;Munkittrick et al. 2023) and metals or ore minerals (Killick et al. 2020).

Significance of study
A basemap of isotope ratios is useful for effectively provenancing biological materials.The term "isoscapes, " coined from "isotope landscapes, " emerged in the mid-2000s to describe maps of isotopic variation generated by predictive models estimating the local isotopic composition of environmental materials (Bowen 2010;West et al. 2010).Despite the increasing availability of regional and global isoscapes for radiogenic heavy isotopes, especially Sr, as well as light stable isotopes of H, O, C, and N, primarily for water and plant samples, challenges persist regarding sampling strategies, mapping and modeling, and interpretation (Bowen 2010;West et al. 2010;Holt et al. 2021, and references therein).
The study of bioavailable isotopes in Korea has only recently gained attention.Strontium isotope data from various sources, including bottled water (Kim et al. 2013), agricultural products (Song et al. 2014), groundwater (Shin et al. 2021), andplant andsoil leachates (Choi et al. 2023), have been used to construct regional isoscapes for southern Korea.Additionally, Park et al. (2015) employed Sr isotope data from human tooth enamel at an archaeological site in southwestern Korea to trace past human migration.Jung et al. (2020) reported on the mineralogical and Sr isotopic composition of soil and plant samples from Jeju Island, situated off the southern coastline, and interpreted the sources of bioavailable Sr.
This study examines the isotopic distribution of bioavailable Sr, Nd, and Pb in Chungcheongbuk-do Province, central Korea, to understand their sources and pathways.This province was selected due to its geological diversity and inland location, which is isolated from coastal influences, making it an ideal region to study the impact of geological and environmental factors on the isoscape.Alongside isotope data, geological features and elevation were considered as external factors in the contour mapping process to construct the isoscapes.The Sr and Nd isoscapes provided are the first local-scale isotopic maps in Korea, holding practical potential for provenance studies as well as archaeological and forensic applications.

Geological and geochemical background
Chungcheongbuk-do Province is situated in the central part of the Korean Peninsula, approximately between 36° 00′ 35" to 37° 15′ 20" latitude and 127° 16′ 40" to 128° 38′ 15" longitude (Fig. 1).It covers an area of approximately 7407 km 2 and has a population of around 1.6 million people.Chungcheongbuk-do does not border the sea, with its closest point to the western sea being approximately 47 km away.Two major rivers, the South Han and Geumgang Rivers, flow through the region.Mountain ranges are located in the east and north, and low-mountainous and hilly areas are found in the southwest (Fig. 1a), shaped by the erosion from the Geumgang River.The province also features numerous basins primarily formed by the South Han and Geumgang Rivers.Karst landforms are prevalent in the northeastern part of the province.
Geotectonically, Chungcheongbuk-do primarily belongs to the Precambrian Gyeonggi Massif and the Neoproterozoic to Paleozoic Okcheon Fold Belt, with some southern parts falling within another Precambrian terrane, the Yeongnam Massif.It encompasses Precambrian metamorphic rocks, Paleozoic sedimentary sequences, and Mesozoic igneous rocks.In this study, the geology of Chungcheongbuk-do is classified into six units based on their ages and lithologies (Fig. 1b), which affect Sr, Nd, and Pb isotopic compositions: Precambrian schist and gneiss complex (unit 1) mainly along the boundary of the province, Paleozoic carbonates (unit 2) in the northeastern part, Paleozoic sediments and metavolcanic rocks (unit 3) extending to the northeast, Triassic to Jurassic granitoid rocks (unit 4) in the southern and northwestern parts, Cretaceous granitoid rocks (unit 5) along the middle eastern boundary of the province, and scattered Cretaceous volcano-sedimentary rocks and dykes (unit 6).

Materials and methods
As depicted in Fig. 1, soil and plant samples were collected from 44 individual sites exhibiting no signs of anthropogenic disturbance.Whenever possible, sampling locations were chosen to be close to those documented in the literature regarding bedrock research.A total of 5, 7, 11, 9, 3, and 9 samples were collected from geological units 1, 2, 3, 4, 5, and 6, respectively.The detailed coordinates of the sampling sites are listed in Table S1.
Soil was extracted from depths exceeding 30 cm below the ground surface using asrm soil auger.Plant samples were individually bagged in Ziploc bags at each site, with multiple leaves collected from locations near the soil sampling points.The majority of leaf samples were obtained from densely wooded areas, including deciduous trees such as Quercus acutissima, Quercus mongolica, Quercus aliena, Quercus serrata, and Quercus dentata, as well as evergreen trees like Pinus densiflora.These samples were homogenized in the sites using the quartering method.Rainwater samples were collected from the rooftop of the Korea Basic Science Institute (KBSI) in Ochang campus on a monthly basis, spanning three periods: November 2021, March 2022, and April 2022, using rainfall collectors (Rain sampler RS-1C, PALMEX).
The accuracy of isotopic measurements in this study was validated through repeated analyses of botanical SRMs, including NIST (National Institute of Standards and Technology) SRM 1515 (Apple Leaves), 1547 (Peach Leaves), 1570a (Spinach Leaves), 1573a (Tomato Leaves), and 1575a (Pine Needles).
All chemical procedures, including sample preparation and instrumental analyses, were carried out at the KBSI Ochang campus.Acid treatment and chemical experiments for the samples were conducted in a clean laboratory area (Class 1000).HCl, HF, and HNO 3 used for sample digestion, dilution and column chemistry were double-distilled by sub-boiling distillation systems.Commercially available high-purity reagents were utilized for HClO 4 (Ultrapure analytical reagent, TAMA CHEMI-CALS), HBr (Suprapur, Merck), acetic acid (CH 3 COOH) (optima grade, Fisher chemical), and ammonium nitrate (NH 4 NO 3 ) (99.999% trace metals basis, Aldrich).Milli-Q water (18 ΜΩ, ELGA Labwater) was used for diluting acids and vial cleaning processes.
Each bulk soil sample was dried overnight at room temperature and then sieved through a 2 mm sieve.A 0.1 g aliquot of each sample was placed into a 60 mL PFA vessel, digested in 3 mL of mixed acid (HF:HNO 3 :HClO 4 = 10:3:1) at 160 °C, and subsequently dried and redissolved in 5 mL of 2.5 N HCl.Additionally, ~ 1 g subsamples were placed into a 15 mL conical tube and leached by adding 2.5 mL of 1 M ammonium nitrate and 2.5 mL of 1 M acetic acid, then shaken for 8 h.The leachates were centrifuged at 10,000 rpm for 15 min, the supernatant (~ 2 mL) was extracted and evaporated to dryness, then redissolved in 5 mL of 2.5 N HCl.Plant samples were washed with Milli-Q water to remove exogenous dust and dried overnight at room temperature.An aliquot of 0.3 g of the sample was digested in 5 mL of a mixed acid (HNO 3 :HF:HClO 4 = 3:1:0.1)at 160 °C, then dried and redissolved in 5 mL of 2.5 N HCl.About 150 mL of rainwater samples were evaporated to dryness and then dissolved in 2.5 N HCl.An aliquot of the liquified subsamples was analyzed using inductively coupled plasma mass spectrometry (ICPMS, Agilent 5800) to determine trace element concentrations.
The Sr, Nd, and Pb fractions for isotope ratio measurement were separated through column chromatography equipped with laminar flow benches.The 2.5 N HCl solutions were evaporated again and redissolved in 0.5 mL of 2.5 N HCl.After centrifugation of the sample solution for 10 min at 10,000 rpm, it was loaded onto a column containing cation exchange resin (H + form, 200-400 mesh, AG50W-X8, Bio-Rad) in a quartz column for the separation of Pb (group), Sr, and REE fractions.For further separation of Nd and Sm in bulk soil samples, plants, and ammonium nitrate leachates, a 0.25 N HCl solution and Ln resin (100-150 μm, Eichrom Technologies Inc.) were used.Pb fractions were further extracted using an anion exchange resin (100-200 mesh, AG1-X8, Bio-Rad) with a 1 N HBr and 6 N HCl solution.
Isotopic measurements were conducted using multicollector inductively coupled plasma mass spectrometry (MC-ICPMS) with a Nu Plasma II instrument.To enhance analyte sensitivity, a desolvating nebulizer system (Aridus III, Teledyne CETAC) was utilized to convert the liquid sample into aerosol form for the measurement of Nd and Pb isotopes.Instrumental operating parameters are detailed in Table S2.Sr, Nd, and Pb isotopes were collected using static mode in one block of 30 cycles, each with an integration time of 4 s per cycle.Instrumental mass fractionation of Sr isotopes was corrected using a factor of 86 Sr/ 88 Sr = 0.1194.Isobaric interference from Kr and Rb was found to be negligible.Replicate MC-ICPMS analyses of NBS 987 resulted in average 87 Sr/ 86 Sr values of 0.710319 ± 0.000012 (n = 94, 2 standard error (SE)).The 87 Sr/ 86 Sr data obtained in this study and from the literature were corrected for interlaboratory bias by adjusting to the recommended NBS 987 value of 0.710248 (Thirlwall 1991).Our laboratory standard solution ([Sr] = ~ 100 ng/g, single component ICP standard solution ICP-55N-1, AccuStandard) gave an average 87 Sr/ 86 Sr of 0.708645 ± 0.000016 (n = 24, 2 SE), which was consistent with the thermal ionization mass spectrometry (TIMS) value of 0.708695 ± 0.000007 (2 SE, internal error).
The Pb isotopic analysis involved correction for instrumental mass fractionation and consideration of mass bias effects.A thallium solution was added to the Pb fraction at a 4:1 ratio of Pb and Tl.The analytical results indicated average values of 206 Pb/ 204 Pb = 16.935 ± 0.006, 207 Pb/ 204 Pb = 15.489 ± 0.005, and 208 Pb/ 204 Pb = 36.701± 0.005 (n = 94, 2 SD), closely aligning with the reported values of NBS 981 (Todt et al. 1996).Total procedural blank levels were below 300 pg for Sr and Nd, and around 1 ng for Pb, which were negligible considering the sample amount.
Spatial analysis and geostatistical interpolation were conducted using ArcGIS Pro software version 3.2.0(Esri, USA).To create a Sr isoscape for Chungcheongbuk-do, previously published data (Song et al. 2014;Shin et al. 2021;Choi et al. 2023) were combined with data obtained in this study.For the Nd isoscape, only data from this study were used.A Pb isoscape is not provided in this study, considering the different origins of bioavailable Pb (as discussed below).The shapefile utilized in the ArcGIS analysis was a revised 1:250,000-scale geological map provided by the Korea Institute of Geoscience and Mineral Resources through the Geo Big Data Open Platform (www.data.kigam.re.kr).River network data were obtained from the Spatial Information Industry Promotion Agency (www.vworld.kr).

Results
The Sr, Nd, and Pb isotope data are presented in Table S3, with Rb, Sr, Sm, Nd, U, Th, and Pb concentrations.It is noted that the concentrations in the soil leachates treated with 1 M ammonium nitrate (AN) and acetic acid (AC) are presented as recovery percentages.The results discussed below pertain only to the data analyzed in this study, unless otherwise specified.

Rb-Sr data
The Rb-Sr data are displayed in Fig. 2 according to the geological subdivision considered in this study.The average concentrations of Rb and Sr across the 44 sites were 74 ± 62 μg/g and 100 ± 82 μg/g in bulk soil, 4.4 ± 4.8 μg/g and 51 ± 46 μg/g in plants, 1.8 ± 3.8% and 3.8 ± 4.4% in AN leachates, and 0.2 ± 0.2% and 2.7 ± 3.3% in AC leachates.As shown in Fig. 2a, bulk soil samples exhibit a broader range of Rb and Sr concentrations compared to plant samples, clustering towards higher Rb concentrations.AN leaching exhibited higher extraction levels of Rb compared to AC leaching.The six geological units overlap in Rb and Sr concentrations and do not exhibit a clear correlation.Rainwater samples have Sr concentrations from 5.37 to 1.87 ng/g.
Most samples were plotted on the 2 Ga to 100 Ma reference lines on the 87 Sr/ 86 Sr versus Rb/Sr plot (Fig. 2b).AN and AC leaching yielded almost identical Sr isotopic compositions, and leachate samples showed a close correlation with the bulk soil samples (R 2 = ~ 0.9) (Fig. S1).The average 87 Sr/ 86 Sr values for AN leachate in each geological subdivision were as follows: 0.732 ± 0.012 (unit 1, n = 4), 0.714 ± 0.004 (unit 2, n = 7), 0.718 ± 0.006 (unit 3, n = 11), 0.714 ± 0.002 (unit 4, n = 9), 0.714 ± 0.001 (unit 5, n = 2), and 0.715 ± 0.002 (unit 6, n = 9).Leachate samples from sites 24 and 25 in unit 5 had distinctly lower 87 Sr/ 86 Sr ratios compared to their corresponding bulk soil samples.As with the bulk soil samples, AN and AC leachates from samples S-32 and S-29 were excluded from the average calculation.The 87 Sr/ 86 Sr ratios of unit 1 showed distinctly higher and more scattered values compared to the other units, mirroring the trend observed in bulk soil samples.
The bulk soil, AN and AC leachate, and plant data are summarized in Fig. 3 using box and whisker diagrams.These diagrams include literature data from 115 bedrocks, 32 plants, and 44 soil leachates combined with data obtained in this study.Groundwater samples from unit 1 exhibit distinctly higher 87 Sr/ 86 Sr ratios (~ 0.75) compared to ground and stream water samples from the other units.Rainwater samples collected in this study approach the lower end of 87 Sr/ 86 Sr ratios.

Sm-Nd data
The concentrations of Sm and Nd were, on average, 5.3 ± 2.7 μg/g and 29 ± 16 μg/g in bulk soil (n = 44), 0.3 ± 0.3 μg/g and 1.7 ± 1.9 μg/g in plants (n = 35), and 3.0 ± 4.6% and 4.2 ± 7.3% in AN leachates (n = 24).As the Fig. 2 Rb-Sr plots.a The correlation of Sr and Rb concentrations for bulk soil and plant samples on a log scale.b 87 Sr/ 86 Sr versus Rb/Sr plot for bulk soil samples.Reference lines corresponding to 2 Ga and 100 Ma are also shown concentration of Sm increases, a corresponding increase in Nd concentration was clearly observed for both bulk soil and plant samples (Fig. 4a).In the AN leachates, the ratios of Sm and Nd concentrations also exhibited a positive correlation, indicating that Sm and Nd were leached in a similar manner.
The ε Nd ranged from -1.8 to -29.3 in bulk soil, -4.4 to -23.8 in plants, and -3.2 to -25.6 in AN leachates.The average ε Nd values for bulk soil in each geological subdivision were: -24.1 ± 4.4 (unit 1, n = 5), -12.7 ± 2.1 (unit 2, n = 7), -15.3 ± 2.2 (unit 3, n = 9), -16.5 ± 1.5 (unit 4, n = 9), -15.3 ± 4.4 (unit 5, n = 3), and -15.1 ± 1.9 (unit 6, n = 9).Unit 1 samples had ε Nd values ranging from -18.2 to -29.3, while the remaining samples mostly had values within -9.9 to -20.4.Samples S-6 and S-10 from unit 3, with ε Nd ratios of -1.85 and -4.06 respectively, were excluded from the average calculation.The plot between 143 Nd/ 144 Nd and the Sm/Nd ratio of bulk soil from unit 1 shows clustering near the reference line corresponding to 3 Ga, while other soil Fig. 3 Box and whisker diagrams for Sr isotopic compositions of bedrock, bulk soil, leachate, plant, and water samples.Different colors identify the six geological units classified in this study.The dashed line indicates the average 87 Sr/ 86 Sr ratio of rainwater samples collected in this study.Note that a different scale was used for values above 0.9 Fig. 4 Sm-Nd plots.a The correlation of Nd and Sm concentrations for bulk soil and plant samples on a log scale.b 143 Nd/ 144 Nd versus Sm/Nd plot for bulk soil samples.Reference lines corresponding to 1, 2, and 3 Ga are also shown samples align with a younger (approximately 2 Ga) reference line (Fig. 4b).Samples from Triassic to Jurassic and Cretaceous granite areas are close to the reference line between 3 and 2 Ga.The samples with the highest Nd isotope values, S-6 and S-10, are positioned near the 1 Ga reference line.The ε Nd values of the AN leachates are generally concordant with those of their corresponding bulk soil samples (R 2 = 0.90) (Fig. S3).

Source of bioavailable Sr, Nd, and Pb
Various sample types have been proposed as bioavailable element archives for building isoscapes, including ground and surface water (Evans et al. 2010;Voerkelius et al. 2010;Ryan et al. 2018;Willmes et al. 2018;Adams et al. 2019), soil and soil leachates (Ryan et al. 2018;Willmes et al. 2018;Adams et al. 2019;Moffat et al. 2020;James et al. 2022;Lugli et al. 2022), plants (Evans et al. 2010;Hartman and Richards 2014;Ryan et al. 2018;Willmes et al. 2018;Adams et al. 2019;James et al. 2022), and insects and animal remains (Evans et al. 2010;Hartman and Richards 2014;Coutu et al. 2016;Adams et al. 2019).Each archive type has its strengths and weaknesses.This study employed plants as the archives of bioavailable elements because they are easily collected and form the base of the food chain, allowing elements to enter the diets of animals and humans.The use of homogenized sampling of deciduous tree leaves, as adopted in this study, can mitigate issues related to point-bias and variations in root systems (Holt et al. 2021).Although rainwater samples exhibit distinctly lower 87 Sr/ 86 Sr ratios compared to bedrock and soil archives (Fig. 3), their contribution to bioavailable Sr is not clearly established.Below, the Sr, Nd, and Pb isotope results from plants are compared with those from bedrock, bulk soil, and soil leachate across geological subdivisions to examine their correlations.
The Sr, Nd, and Pb isotopic compositions of bedrocks are governed by their formation ages and parent/daughter isotope ratios.The relatively high 87 Sr/ 86 Sr and low ε Nd values of unit 1 rocks, despite the limited number of samples, reflect their geologically old (~ 2 Ga) formation ages in the crust.Conversely, the comparatively high 87 Sr/ 86 Sr ratios of unit 5 samples, which are geologically young (~ 90 Ma; Cheong and Chang 1997) and highly fractionated (i.e., high in SiO 2 ) granites, result from their anomalously high Rb/Sr ratios (188-28; Cheong and Chang 1997;Lee et al. 2010).The ε Nd values of unit 5 bedrocks are similar to those of the other unit bedrocks.The reported 87 Sr/ 86 Sr ratios of unit 2 bedrocks (0.716-0.707) are distinctly different from those expected from the seawater value (McArthur et al. 2001), indicating alteration or a metamorphic imprint on the samples.The bedrock isotopic compositions are not directly reflected in biomaterials because only the labile fractions of the bedrocks are transferred to the ecosystem.Despite this, as displayed in Figs.8a and b, the bedrock Sr and Nd isotopic compositions from Chungcheongbuk-do are moderately compatible with bioavailable values represented by plant data (R 2 = ~ 0.8), suggesting a broad geological control on biomaterial compositions.In contrast, Pb isotope data between plant and bedrock samples show no correlation (Fig. 8c).
The Sr and Nd isotopic compositions of bedrock, bulk soil, and plant samples are not correlated with each other (Fig. S5).This lack of correlation primarily arises because different geological units belong to separate isotope systems with varying initial daughter isotope compositions and ages.On the other hand, the modeled Sr, Nd, and Pb ages of bulk soils (~ 3-2 Ga; Figs.2b, 4b, 7a) suggest a geological control from the bedrocks.
Sr, Nd, and Pb in the minerals composing the bedrocks and weathered soil profile are not equally transferred to the ecosystem.For example, Sr-rich carbonates and plagioclase, which weather rapidly, would dominate the exchangeable Sr pool, particularly in young soils (Vitousek et al. 1999;Chadwick et al. 2009).This selectivity could explain the difference in 87 Sr/ 86 Sr between the bulk soil and leachate/plant data for sites 24 and 25 in unit 5. Compared to bedrock data, bulk soil data are more closely correlated with plant data (Figs.8d, e).The R 2 values of 87 Sr/ 86 Sr and ε Nd data are calculated to be 0.83 and 0.87, respectively.Like the bedrock data, Pb isotope data for bulk soil are not correlated with the plant data (Fig. 8f ).The ammonium complex solution primarily facilitates the desorption of cations, effectively suppressing the formation of colloids that are not absorbed by plants, thereby enabling the extraction of bioavailable trace elements that are transferable to plants.Acetic acid is primarily used to dissolve carbonates, allowing the tracing of the origin of carbonates within the soil (DePaolo et al. 1983;Capo et al. 1998;Gryschko et al. 2005).As mentioned earlier, AN and AC leachates in this study yielded nearly identical Sr and Pb isotope data (Figs.S1, S4).The correlation of Sr (R 2 = 0.98) and Nd isotope data (R 2 = 0.89) is prominent between the AN leachate and plant data (Figs. 8g,h), confirming an important contribution of labile soil fractions to bioavailable Sr and Nd.However, no Pb isotopic correlation is observed in this comparison (Fig. 8i).The distinctly different Pb isotopic compositions of plants, compared to those of bulk soil samples and their leachates, suggest that the source of bioavailable Pb includes other environmental archives, likely including historical gasoline usage and/or mining-related pollution.Identifying the precise Pb sources requires further study.

Mapping and modeling
The principal methods for constructing isoscapes involve isotopic estimation and geostatistical modeling.Empirical isotope data can be collected from various sample types representing bioavailable element pools within a specific study area.When empirical data are scarce, predictive models of radiogenic isotope composition can be used.These models are based on the geochemical knowledge of bedrocks or waters and employ the basic equation of isotopic evolution (Beard and Johnson 2000;Bataille and Bowen 2012).While this approach is relatively straightforward and does not require extensive empirical data, it is limited to predicting isotope ratio distributions over large regions or globally.Since this study aims to create a practical, finescale isoscape, this mechanistic approach (Bataille et al. 2020) is not considered here.Additionally, machine learning is increasingly being used to integrate isotope data from different sources based on training data (Bataille et al. 2018).However, this approach also provides rough estimations over broad regions and is therefore not considered in this study.
Mapping protocols are classified into domain and contour mapping.Domain mapping relates the analyzed isotope ratios of targeted archives to specific areas.In this study, this means associating plant data with geological units.This approach, as exemplified by the domain maps of Evans et al. (2010) and Kootker et al. (2016), is useful in regions where superficial deposits largely follow the underlying geology and exhibit unique isotopic compositions.However, as illustrated in Figs. 3 and 5, the Sr and Nd isotopic compositions of plants from the subdivided geological units in Chungcheongbuk-do overlap significantly, making this approach less applicable for provenance studies.Contour mapping methods estimate unknown values at specific locations using known values from surrounding sample points.These methods include inverse distance weighting (IDW), simple kriging, empirical Bayesian kriging (EBK), and kriging with external drift (Emery et al. 2018;Wang et al. 2018;Willmes et al. 2018;Adams et al. 2019).In this study, these approaches were used to create Sr and Nd isoscapes for Chungcheongbuk-do using plant and soil data obtained in this study and from the literature.For sites where both plant and soil samples were available, the plant values were used as the site values.For locations in the literature with multiple reported values for either plants or soils, the average value was used.For Sr, the isoscape was constructed using values from 44 sites in this study and 46 sites from previous studies (32 sites with plant data, and 14 sites with soil leachate data), totaling 90 data points.For Nd, the isoscape was based on plant values measured at 35 sites in this study.Ground and surface water data were not used for building the isoscapes because they may represent integrated large areas rather than specific sample locations.A Pb isoscape is not included in this study because it is evident that bioavailable Pb comes from extraneous environmental sources, as explained earlier with Pb isotope data (Figs.7 and 8).
The IDW method operates under the assumption that closer locations share greater similarity, thus assigning higher weights to neighboring points.For the IDW employed in this study, the power parameter was set to 2, and the neighborhood was configured with a maximum of 15 and a minimum of 10 points, divided into eight sectors.The IDW-based Sr and Nd isoscapes are presented in Figs.9a and b, respectively.They exhibit the Bull's eye effect resulting from the prominent influence of nearby points.This effect is shown as concentric circles of varying colors around measurement points, which is commonly observed in IDW maps.As displayed in Fig. 9, the western and northeastern parts of Chungcheongbukdo are marked by high 87 Sr/ 86 Sr ratios and low ε Nd values.The geological classification adopted in this study (Fig. 1b) only weakly explains the isotopic distribution in these maps.
Kriging is a geostatistical interpolation method that provides probabilistic predictions.The predicted values align with the known values on average, minimizing prediction errors.The semivariogram models how sample values change with distance and direction, reflecting spatial autocorrelation.Simple kriging assumes a known constant mean throughout the study area and uses this assumption to predict values at unmeasured locations.For simple kriging applied in this study, the data were non-transformed, and a nugget effect was included in the semivariogram.The neighborhood was configured with a maximum of 15 and a minimum of 5 points, divided into eight sectors.These settings were chosen to minimize the root mean square error (RMSE).The Sr and Nd isoscapes based on simple kriging, as shown in Figs.10a-d, generally follow the features displayed by the IDW maps.The prediction error map from simple kriging reveals pronounced concentric circles around the data points, with higher error values in regions with low sample density, such as the central western area of the Nd map (Fig. 10d).(Krivoruchko 2012).In this study, the EBK function was set with 100 simulations, and the data were transformed none.The neighborhood type for EBK was configured with a maximum of 15 and a minimum of 10 points, to minimize the RMSE.The EBK isoscapes are presented in Fig. 10e-h.Compared to the simple kriging maps, the EBK maps show similar patterns but do not emphasize single data points as strongly.The predicted 87 Sr/ 86 Sr ratio has a comparatively narrow range (0.756-0.712), compared to the IDW and simple kriging models.The prediction error in EBK is more gradually expressed as distance from the data points increases, compared to simple kriging.While most prediction errors around data points are low, in regions where isotope ratios change rapidly over a small spatial area-such as areas with high 87 Sr/ 86 Sr ratios (P-32 from unit 1 and P-29 from unit 5, see Fig. 10e)-the errors (Fig. 10f ).
It should be noted that data smoothing in contour mapping creates geologically unrealistic gradual gradient boundaries between areas with different lithologies.In nature, very different isotopic compositions can exist in close proximity.To address this issue, we used a type of kriging with external drift protocol, the EBK regression prediction (EBKRP).This method involves adding rasterized weighted factor data into the interpolation dataset and utilizing regression analysis.Data analysis employed an exponential semivariogram type with 100 simulations, a maximum of 15 and a minimum of 10 neighborhood points, and no data transformation.
Geological features were used as the explanatory variable in the EBKRP protocol.The EBKRP-based Sr and Nd isoscapes, considering geological classification (Fig. 1b), are presented in Fig. 11a-d.The maps closely reflect the geological features, and the predicted 87 Sr/ 86 Sr ratio exhibiting a notably limited range (0.730-0.714).Geology and topography are closely related: Mesozoic granites, being heavily weathered, typically form low basins, while Precambrian gneisses tend to form high mountains, as observed in the north and east of Chungcheongbuk-do (Fig. 1).The prediction error of the Sr isoscape exhibited a narrower range of error values (0.0075-0.0033) and more gradual variations around the data points.Similar to simple kriging and EBK, regions with fewer data points in the northern and southern areas show higher errors in the Sr and Nd isoscapes.When both geological classification and elevation are used together as weighted factors, the isoscapes are displayed as shown in Fig. 11e-h.In this case, the 87 Sr/ 86 Sr range (0.762-0.710) is much wider than that of the geology-only Sr isoscape.As expected from the geological classification, high 87 Sr/ 86 Sr areas correspond approximately to the boundary of the province.Conversely, Cretaceous volcano-sedimentary rocks and dykes (unit 6) show distinctly lower 87 Sr/ 86 Sr.The Nd isoscape, almost identical to the geology-only map, reflects the geological features more strongly compared to the geology/topography-considered Sr isoscape.The prediction error patterns were similar to those observed in the geology-considered maps.The maps in Fig. 11 await validation through the analysis of ancient or modern faunal archives, particularly small animals that are mobile within a limited range (Price et al. 2002).These maps hold potential for provenance studies and archaeological/forensic applications in Chungcheongbuk-do as baseline maps in their present form.However, further efforts are necessary to accurately predict fine isotopic variations, including increasing sample density, especially in the central and southern parts of the province, and refining interpolation models.

Fig. 1
Fig. 1 Geographic and geological information on Chungcheongbuk-do, central Korea, with sample locations.a Digital elevation model of Chungcheongbuk-do Province derived from ASTER satellite data captured on 30 November 2013, illustrating the topography and indicating the sample site.The black-filled area in the inset figure indicates the geographic position of Chungcheongbuk-do.b Simplified geological subdivision of the province with the sample sites.Stream networks are also shown

Fig. 5
Fig. 5Box and whisker diagrams for Nd isotopic compositions of bedrock, bulk soil, leachate, and plant samples.Different colors identify the six geological units classified in this study

Fig. 6
Fig. 6 Th versus U (a) and Pb versus U (b) concentration plots for bulk soil and plant samples on a log scale

Fig. 8
Fig. 8 Plant Sr, Nd, and Pb isotope data compared with bedrock (a-c), bulk soil (d-f), and ammonium nitrate leachate (g-i) data.Regression lines and their equations are shown for Sr and Nd data.1:1 lines are shown in the Pb isotope plots