Mineral Element Stocks in the Yedoma Domain: A Novel Method Applied to Ice-Rich Permafrost Regions

With permafrost thaw, significant amounts of organic carbon (OC) previously stored in frozen deposits are unlocked and become potentially available for microbial mineralization. This is particularly the case in ice-rich regions such as the Yedoma domain. Excess ground ice degradation exposes deep sediments and their OC stocks, but also mineral elements, to biogeochemical processes. Interactions of mineral elements and OC play a crucial role for OC stabilization and the fate of OC upon thaw, and thus regulate carbon dioxide and methane emissions. In addition, some mineral elements are limiting nutrients for plant growth or microbial metabolic activity. A large ongoing effort is to quantify OC stocks and their lability in permafrost regions, but the influence of mineral elements on the fate of OC or on biogeochemical nutrient cycles has received less attention and there is an overall lack of mineral element content analyses for permafrost sediments. Here, we combine portable X-ray fluorescence (pXRF) with a bootstrapping technique to provide i) the first large-scale Yedoma domain Mineral Concentrations Assessment (YMCA) dataset, and ii) estimates of mineral element stocks in never thawed (since deposition) ice-rich Yedoma permafrost and previously thawed and partly refrozen Alas deposits. The pXRF method for mineral element quantification is non-destructive and offers a complement to the classical dissolution and measurement by optical emission spectrometry (ICP-OES) in solution. Using this method, mineral element concentrations (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr and Zr) were assessed on 1,292 sediment samples from the Yedoma domain with lower analytical effort and lower costs relative to the ICP-OES method. The pXRF measured concentrations were calibrated using alkaline fusion and ICP-OES measurements on a subset of 144 samples (R 2 from 0.725 to 0.996). The results highlight that i) the mineral element stock in sediments of the Yedoma domain (1,387,000 km2) is higher for Si, followed by Al, Fe, K, Ca, Ti, Mn, Zr, Sr, and Zn, and that ii) the stock in Al and Fe (598 ± 213 and 288 ± 104 Gt) is in the same order of magnitude as the OC stock (327–466 Gt).


INTRODUCTION
The ice-rich deposits of the Yedoma domain hold more than 25% (213 Gt) of the frozen organic carbon (OC) of the northern circumpolar permafrost region, while covering only about 8% of its total soil area (c. 1.4 of 17.8 million km 2 ; Schirrmeister et al., 2013;Hugelius et al., 2014;Strauss et al., 2017). This carbon-and ice-rich region is particularly vulnerable to abrupt thawing processes (Nitzbon et al., 2020;Turetsky et al., 2020). This is the reason why Yedoma deposits are considered as a potential "tipping element" for future climate warming (Lenton, 2012). The key features of Yedoma deposits relative to other permafrost deposits are i) their ground ice properties, with 50-90% volume percent ice  and ii) their large spatial extent and thickness, resulting in a large total volume . Thawing of ice-rich sediments has severe geomorphological consequences for the landscape (Kokelj and Jorgenson, 2013). Striking examples are collapsing river and coastal bluffs (Günther et al., 2015;Kanevskiy et al., 2016;Fuchs et al., 2020). But a spatially more important process is thermokarst lake formation, caused by surface subsidence due to excess ground ice melting in Yedoma deposits and leading to the formation of vast thermokarst depressions. The process not only affects the superficial soil horizons but also deep mineral horizons (Strauss et al., 2013;Walter Anthony et al., 2018). Thermokarst processes were highly active during the Deglacial and early Holocene period (Walter et al., 2007;Brosius et al., 2021) and shaped the Yedoma domain region with numerous lakes. Following drainage of lakes, a widespread process across lakerich lowlands , renewed permafrost aggradation and new soil development could occur in the drained thermokarst lake basins also called Alas. The Yedoma domain therefore includes Yedoma deposits never affected by thaw and frozen deposits that accumulated after Yedoma degradation in Alas landforms (Olefeldt et al., 2016;Strauss et al., 2017).
In permafrost soils, between ∼30% (Mueller et al., 2015) and ∼80% (Dutta et al., 2006) of the total soil OC is associated to minerals through various mechanisms, as detailed below. It is well known that interactions between minerals and OC can have a stabilizing effect on OC (Kaiser and Guggenberger, 2003;Lutzow et al., 2006;Kögel-Knabner et al., 2008;Gentsch et al., 2018;Wang et al., 2019). Mineral-protected OC is less bioavailable than mineral free OC (Schmidt et al., 2011;Kleber et al., 2015), thereby contributing to the long term carbon storage (Hemingway et al., 2019). Mineral protection of OC includes aggregation, adsorption and/or complexation or co-precipitation processes (Kaiser and Guggenberger, 2003). These protection mechanisms involve i) clay minerals, Fe-, Al-, Mn-oxy-(hydr) oxides, or carbonates (aggregation); ii) clay minerals and Fe-, Al-, Mn-oxy-(hydr)oxides using polyvalent cation bridges such as Fe 3+ , Al 3+ , Ca 2+ or Sr 2+ (adsorption); or iii) Fe 3+ , Fe 2+ and Al 3+ ions (complexation). In addition, mineral elements (e.g., Si, Al, Fe, Ca, K, Mn, Zn, Sr) drive nutrient supply to plants or microorganisms, including algae such as diatoms. Nutrient supply is essential for plants growth and/or microbes metabolic activity and therefore indirectly influences C storage or degradation. It remains however unclear how mineral-OC interactions and nutrient supply will evolve upon permafrost thaw (Opfergelt, 2020), especially in ice-rich permafrost regions. This is mainly due to a lack of knowledge about the mineral element content in permafrost regions, relative to the well-known OC stocks in permafrost soils (Harden et al., 2012;Hugelius et al., 2014;Strauss et al., 2017;Hugelius et al., 2020;Kuhry et al., 2020) and the increasing knowledge on N stocks Fuchs et al., 2019;Fouché et al., 2020;Hugelius et al., 2020).
To improve our understanding on the potential effects mineral elements can have on OC from permafrost regions, it is essential to have better knowledge on mineral element stocks in these regions, and on the changes in mineral element stocks upon thawing with changing conditions for mineral weathering (Zolkos and Tank, 2020). Ice-rich permafrost regions are particularly well suited areas to study the impact of thawing processes on the mineral elemental stocks in deposits. Indeed, "never" (since deposition) thawed deposits can be compared with previously thawed deposits resulting from thermokarst processes during warmer and wetter periods in the Deglacial and early Holocene periods (13-9 ka BP; Morgenstern et al., 2013;Walter et al., 2007). Assessing the evolution of the mineral element stocks between never thawed (since deposition) and previously thawed, refrozen as well as newly formed deposits will contribute to a better understanding of the impact of past thawing processes on the evolution of the pool of mineral elements available for mineral-OC interactions. It also provides insights into how ongoing permafrost thaw and thermokarst processes may impact mineral elements and what the potential implications are for the fate of OC in ice-rich deposits . This is particularly relevant given that thermokarst processes are projected to spread across the Arctic and will potentially unlock additional OC stocks (Lacelle et al., 2010;Abbott and Jones, 2015;Schneider von Deimling et al., 2015;Nitzbon et al., 2020;Turetsky et al., 2020).
The objective of this study is to provide a method to assess the mineral element content on a large sample set in order to cover the different types of deposits from the Yedoma domain. Using this method, we provide a Yedoma domain Mineral Concentrations Assessment (YMCA) dataset for ten mineral elements (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr, and Zr) and the stocks for these mineral elements in deposits from ice-rich permafrost regions. This assessment comprises "never" thawed since syngenetic freezing Yedoma Ice Complex deposits (in the following referred as Yedoma) and at least once previously thawed and then refrozen drained thermokarst lake basin deposits (in the following referred as Alas).

Environmental Settings
The Yedoma domain is part of the permafrost region characterized by organic-and ice-rich deposits as well as thermokarst features. Today, the Yedoma domain predominantly covers areas in Siberia and Alaska ( Figure 1) Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 703304 which were not covered with ice sheets during last glacial period (110 ka -10 ka BP; Schirrmeister et al., 2013).
Recent estimates indicate that the current Yedoma domain has an extent of ∼1.4 million km 2 and contains between 327 and 466 Gt OC . Frozen deposits from the Yedoma domain alone store probably at least as much carbon as the tropical forest biomass (Lai, 2004). These deposits were formed by long-term continuous sedimentation and syngenetic freezing. Yedoma deposits were first described as homogeneous silty fine, ice-, and organic-rich sediments derived from aeolian processes (i.e., loess or loess-related deposits) (Pewe and Journaux, 1983;Tomirdiaro and Chernen'kiy, 1987;Murton et al., 2015). Current evidence indicates that depositional processes are polygenetic including alluvial and aeolian deposition and re-deposition, as well as in situ weathering during the late Pleistocene cold stages (Konishchev and Rogov, 1993;Schirrmeister et al., 2013). Despite evidence of homogeneity of Yedoma deposit aggradation, grain-size analyses show the large diversity of Yedoma deposits, pointing at multiple origins and transport pathways of sediments as well as a diverse interplay of (post)depositional sedimentary processes (Strauss et al., 2012;Schirrmeister et al., 2020;Sizov et al., 2020). For tens of millennia, continuous sedimentation led to the accumulation of several tens of meters thick permafrost deposits with characteristic syngenetic ice-wedge formation. Harsh, cold late Pleistocene winter climate triggered the formation of vertical frost cracks within surface deposits in which snow melt water percolated during spring and refroze into vertical ice veins in the permafrost. Over time and with annual repetition of this process, many meters wide and up to 40 m high ice-wedges were formed. Those large volumes of icy sediments together with the rise in temperature during the Holocene Thermal Maximum (HTM between 11 and 5 ka BP depending on the region; HTM reached at 7.6-6.6 ka BP in North Alaska and variable with time in Siberian regions) (Velichko et al., 2002;Kaufman et al., 2004;Porter and Opel, 2020) lead to a vast reshaping of the landscape with formation of thermokarst lakes and drained Alas basins . Alas deposits in drained thermokarst lake basins are composed of reworked Yedoma deposits as well as Holocene accumulation during sub-aquatic and sub-aerial phases Windirsch et al., 2020). The Yedoma domain area includes 30% of unthawed Yedoma deposits composed of homogeneous silty deposits with polygenetic origins (eolian, alluvial or colluvial deposition) between large ice-wedges. Alas deposits have experienced thawing processes and drainage during early Holocene until more recently and account for 56% of the Yedoma domain area. Deltaic deposits (4%) as well as lakes and rivers (10%) complete the Yedoma domain deposits area distribution ( Figure 2).
In this study, mineral element stocks are evaluated in Yedoma and Alas deposits (for a mean deposit thickness of 19.6 and 5.5 m, respectively; Figure 3A), but the evaluation does not include sediments underlying extant thermokarst lakes, active layer sediments or fluvial sediments ( Figure 3B).

Sample Collection
This large-scale mineral element concentration assessment from Yedoma and Alas deposits is based on samples collected in two major Arctic regions, Alaska and Siberia. More specifically, the dataset includes 22 locations and compiles 75 different profiles from West, North and Interior Alaska, the Kolyma region, the Indigirka region, the New Siberian Archipelago, the Laptev Sea coastal regions, and Central Yakutia (Figure 1). About 64% of the profiles are from Yedoma deposits, 33% from Alas deposits and 3% from deltaic deposits, and the sites cover a range of geomorphological settings (Supplementary Table S1). For each location, Yedoma and Alas deposits profiles were sampled if both types of deposit were present. In total, 1,292 different samples were analyzed for their mineral element concentration. Sampling strategies were as follows: during the cold season, frozen samples were collected by drilling from the surface down with drilling rigs, whereas during the mid-summer season, drilling was performed in frozen ground below the existing thawed active layer. Cliffs or headwall exposures along coasts, rivers, or lakes were cleaned and frozen in situ sediments sampled with hammer and ax. For headwall or cliff sampling, sub-profiles from different vertical exposures were included when needed to reconstruct a complete stratigraphical composite profile. Additional information on specific site sampling techniques can be found in the reference papers cited in Table 1. Recovered samples were air-or freeze-dried before being stored and archived. The number of samples collected in each location and the number of samples analyzed with the different methods presented in Methods are provided in Supplementary Table S1.

Total Elemental Analysis by ICP-OES Measurement After Alkaline Fusion
Inductively coupled plasma optical-emission spectrometry (ICP-OES) is a classical method to assess mineral element concentrations in solutions from the environment accurately. Solid phases such as soil samples should be first digested prior to ICP-OES analysis. Here, soils were digested by alkaline fusion. Briefly, air-or freeze-dried soil samples were carefully milled for homogenization (agate mill). We mixed a portion of the milled sample (80 mg) with Lithium metaborate and Lithium tetraborate and heated it up to 1,000°C for 10 min. Then we dissolved the fusion bead in HNO 3 2.2N at 80°C and stirred until complete dissolution (Chao and Sanzolone, 1992). We measured the mineral element concentrations in that solution by ICP-OES (iCAP 6500 Thermo Fisher Scientific). We assessed the loss on ignition at 1,000°C, the sum of oxides was between 98 and 102%, and the total element content in soils is expressed in reference to the soil dry weight at 105°C. To assess the accuracy (i.e., trueness and precision) of the method, the analytical measurement was validated by repeated measurements on three certified materials; i) USGS basalt reference material BHVO-2 (Wilson, 1997), ii) GBW07401 podzolitic soil, and iii) GBW07404 limy-yellow soil (National Research Center for CRM, 1986; Supplementary Table  S2). To assess the precision of the method on a matrix similar to samples from the Yedoma domain, we conducted three  repetitions on three individual samples from Yedoma and Alas deposits. For each mineral element, standard deviations on the repetitions, expressed in mg kg −1 , are available in Table 2.
Out of the 1,292 deposit samples retrieved from the Yedoma domain (Sample Collection), a subset of 144 samples was analyzed by ICP-OES after alkaline fusion to determine their mineral element concentrations (sample list is provided in Supplementary Table S3). We used this subset of analysis to calibrate element concentrations measured by portable X-ray fluorescence (Total Elemental Analysis by ex situ Direct Measurement by Portable X-Ray Fluorescenc) for accurate determination of concentration values. For this first assessment, we measured the following elements on the subset of samples (except for Zn, which was measured on 119 out of the 144 samples): Si, Al, Fe, Ca, Mg, Cr, Ba, K, Ti, P, Cu, Mn, Ni, Zn, Sr, and Zr.
Total Elemental Analysis by ex situ Direct Measurement by Portable X-Ray Fluorescence X-ray fluorescence (XRF) spectrometry is an elemental analysis technique with broad environmental and geologic applications, from pollution assessments to mining industries (Weindorf et al., 2014a;Rouillon and Taylor, 2016;Young et al., 2016;Ravansari et al., 2020). In addition, there is a growing use of portable XRF 1 | Studied locations from the Yedoma domain with associated labels, the number of profiles assessed for each type of deposits and the associated reference papers for site description. The site numbers (Site Nb) 1-17 are from Siberia, and 18-22 are from Alaska (located on the map in Figure 1).   a Sobo Sise Yedoma profiles were collected from the top of Yedoma hills but radiocarbon dating indicates that those are actually cover deposits of Holocene age.
TABLE 2 | Precision of the element concentrations expressed as pooled standard deviations (i.e., two pooled standard deviations, expressed in mg kg −1 ) for the 10 elements considered. The values are based on three repetitions of three identical samples for both ICP-OES and raw concentrations from pXRF method. The coefficient of variation (CV; expressed in %), defined as the ratio of the standard deviation to the mean is also provided. for soil science (McLaren et al., 2012;Ravansari and Lemke, 2018). Portable XRF is used primarily for solid elemental analysis (soils, sediments, rocks or even plastics) but can also deal with oil chemical characterization (Weindorf et al., 2014a). XRF is based on the principle that individual atoms emit photons of a characteristic energy or wavelength upon excitation by an external X-ray energy source. By counting the number of photons of each energy emitted from a sample, the elements present may be identified and quantified (Anzelmo and Lindsay, 1987; Supplementary Figure S1). XRF-scanning results represent elements intensities in "counts per second" (cps) which are proportional to chemical concentrations in the sample but depend also on sample properties (Röhl and Abrams, 2000), ice and water content (Tjallingii et al., 2007;Weindorf et al., 2014b) and interactions between elements called "matrix effect" (Weltje and Tjallingii, 2008;Fritz et al., 2018).

Si
In situ pXRF measurement often involves variability from uncontrolled environmental factors, such as water content, organic matter content or sample heterogeneity (Shand and Wendler, 2014;Weindorf et al., 2014b;Ravansari et al., 2020). To avoid such variability in water content, measurements were performed on dried samples in laboratory (ex situ) conditions with the handheld device in the laboratory. The particle size distribution of these Yedoma domain deposits (described in reference papers from Table 1) is below 2 mm: therefore no sieving was necessary. For the pXRF measurement, the dried sample is placed on a circular plastic cap (2.5 cm diameter) provided at its base with a transparent thin film (prolene 4 µm). To avoid the underestimation of the detected intensities, sample thickness in the cap needs to be higher than 5 mm to 2 cm, depending on the element of interest (Ravansari et al., 2020). For a precise measurement, the sample thickness in the cap is set to >2 cm. Above 2 cm, the width is considered as "infinitely thick" for all elements. Measurements on the 1,292 samples from the Yedoma domain were performed using the pXRF device Niton xl3t Goldd+ (Thermo Fisher Scientific), which has two specific internal calibration modes called mining and soil. Each internal calibration is dealing with different energy range and filters to scan the complete energetic spectrum from low to high-energetic fluorescence values. Both modes were used on each sample and depending on the element, the calibration with the best correlation with the ICP-OES method (Total Elemental Analysis by ICP-OES Measurement After Alkaline Fusion) was kept for further calculations (Mineralogy by X-Ray Diffraction).
To standardize the analysis, total time of measurement was set to 90 s. We conducted the analysis in laboratory conditions, using a lead stand to protect the operator from X-rays.
In theory, the pXRF device used to generate this dataset can measure simultaneously elements of atomic mass from Mg to U. Because ambient air annihilates fluorescence photons that do not have enough energy, low atomic mass elements from Na and lighter cannot be quantified by pXRF. Note that Na quantification would be possible in controlled void conditions during analysis (Weindorf et al., 2014a). In this study, we focussed on the concentrations in 16 elements (Si, Al, Fe, Ca, Mg, Cr, Ba, K, Ti, P, Cu, Mn, Ni, Zn, Sr, Zr) by pXRF and by the ICP-OES method (Total Elemental Analysis by ICP-OES Measurement After Alkaline Fusion) to allow for quality check, calibration and correction. Some elements are at the limit of detection (LOD) for pXRF device (e.g., Cu, Ni). The LOD was reached when the sample was lacking a specific mineral element. LOD concentrations were set to 0.7 times the minimal concentration measured for this element, which is an arbitrary number but conventionally used for data statistical treatment (Reimann et al., 2008). Depending on the considered element, pXRF measured concentrations were highly precise but not always accurate (far from the true value) (Figure 4). Trueness was achieved after correction using a regression with concentration values measured by the ICP-OES method to avoid systematic bias (Figure 4; Linear Regression for Accurate Mineral Element Concentration Measurements). Using a well-defined regression to correct pXRF measurements for trueness allowed using the pXRF method to measure mineral element concentrations on a large number of Yedoma and Alas sample (n 1,292), i.e., providing a valuable method to assess the mineral element content on a large sample set ( Figure 4). Raw pXRF concentrations cannot be used for absolute quantification if not corrected with a reliable and accurate method. Here, only pXRF concentration values corrected using a well-defined regression were used (Linear Regression for Accurate Mineral Element Concentration Measurements). Because of poor correlation (R 2 < 0.5) between pXRF and accurate ICP-OES method for Mg, Cr, Ba, P, Cu and Ni, as discussed in Linear Regression for Accurate Mineral Element Concentration Measurements, these elements were not considered further in this study.
To assess the precision (i.e., random errors) of the pXRF method, three repetitions were conducted on three individual samples from different locations. Between each repetition, instrument/sample repositioning is used to mitigate the "nugget effect" due to heterogeneity in the sample (Ravansari and Lemke, 2018). Pooled standard deviations (SD pooled , weighted average of standard deviations), expressed in mg kg −1 , of the repetitions for the ten elements used for stock calculation (Linear Regression for Accurate Mineral Element Concentration Measurements) are available in Table 2. Given the influence of sample matrix on pXRF measurements, the precision of the pXRF method was also evaluated based on three to five repetitions on 20 individual samples and on average 2.6 times larger than based on three repetitions on three samples (Supplementary Table S4). The coefficient of variation was 20% smaller for Al but 6.5 times larger for Ca based on 20 samples. To ensure a cautious evaluation of the dataset, we decided to report the precision on the data in Supplementary Presentation S1 based on the precisions from Table 2 for ICP-OES measurements, and from Supplementary  Table S4 for pXRF measurements to use the largest set of sample with precision data available.

Mineralogy by X-Ray Diffraction
The X-ray diffraction (XRD) method allows the characterization of the presence of crystalline mineral phases. We assessed the mineralogy of 39 finely ground bulk samples and six clay fraction samples (samples selected in Siberia and Alaska detailed in September 2021 | Volume 9 | Article 703304 6 Supplementary Table S1). The mineralogy of bulk samples was determined on powder (Cu Kα, Bruker Advance D8). Clay fraction mineralogy was assessed after K + and Mg 2+ saturation, ethylene glycol (eg) solvation and thermal treatments at 300 and 550°C (Robert and Tessier, 1974). The clay size fraction (<2 µm) was recovered after sonication, sieving at 50 µm to remove the sand fraction, and dispersion with Na +resins to separate silt (2-50 µm) and clay fractions (Rouiller et al., 1972).

Mean-Bootstrapping Technique for Mineral Element Stocks Calculations
Calculations of mineral element stocks are based on a meanbootstrapping technique. From the corrected mineral element concentrations obtained via linear regression of pXRF measurements (YMCA dataset; Yedoma Domain Mineral Concentrations Assessment (YMCA) Dataset), the aim is to calculate the mineral element stock in Yedoma and Alas deposits individually, in order to estimate total mineral element budget at the Yedoma domain scale. This technique has been used to estimate OC budget from Yedoma and Alas deposits (Strauss et al., 2013) and was further improved by Jongejans and Strauss (2020). Because Yedoma domain deposits contain large ice volumes, stock calculations need to take into account not only the bulk density (BD) of the samples but also total thickness and total wedge-ice volume (WIV) of the deposits. In order to avoid overestimation of the mineral stocks, WIV was subtracted from the total Yedoma domain deposit volume, since the proportion of mineral elements locked inside ice-wedges is negligible. Indeed, Opel et al. (2018) indicate the minor contribution of mineral particles during ice-wedge formation. Fritz et al. (2011) further indicate that ion concentrations in ice is generally low, dominated by HCO 3 − (55%) and Cl − (37%) for anions and Na + (58%) and Ca 2+ (30%) for cations. Note that a flush of highly labile mineral elements (e.g., Na + , Ca 2+ , Cl − ) locked inside ice-wedges may increase nutrient supply for a short time scale upon Yedoma deposits degradation compared to long-term solid-liquid interactions upon thaw. Assumptions and calculations for BD determination, WIV estimations, deposits thickness, and coverage are fully explained in Strauss et al. (2013) and are summarized here. The bootstrapping statistical method use resampled (10,000 times) observed values (mineral element concentrations, BD, WIV and deposits thickness (Supplementary Figure S2)) and derive the mean afterward. The BD determination for each sample was obtained by using an inverse relationship with porosity (ϕ) Eq. 1 (Strauss et al., 2013): whereas ρ s is the solid fraction density. We assume that pore volume in ice-saturated samples is directly measured with segregated ice volume. For samples where no BD determination is available, BD is inferred from the total If neither BD or TOC is available, BD is fixed to 0.88 103 and 0.93 103 kg m −3 for Yedoma and Alas deposits, respectively, i.e., the mean BD measured in such deposits (Strauss et al., 2013). These values are comparable with other studies on BD of Yedoma deposits (0.98 103 kg m −3 ; Dutta et al., 2006). The WIV estimations were performed using ice-wedge width from field measurements (Strauss et al., 2013), ice-wedge polygon size determinations from high-resolution satellite, and additional geometrical tools (assuming ice-wedges have an isosceles triangular or rectangular shape depending on the type of icewedge; see Strauss et al. (2013), Ulrich et al. (2014)). Thickness used for mineral element stock estimations in Yedoma domain deposits are based on mean profile depths of the sampled Yedoma (n 19) and Alas (n 10) deposits (Table 3; Strauss et al., 2013). Since the 10,000 step bootstrapping technique randomly picks one thickness at a time with replacement, we evaluated the stock estimation for a mean thickness of 19.6 m deep in Yedoma deposits and 5.5 m deep in Alas deposits.
The core Yedoma domain extent was estimated to ∼1,387,000 km 2 , based on digital Siberian Yedoma region map (Romanovskii, 1993) and the distribution of Alaskan ice-rich silt deposits equivalent to Yedoma (Jorgenson et al., 2008). The Yedoma deposit extent is estimated to 410,000 km 2 , i.e., 30% of the Yedoma domain, based on the fact that 70% of the Yedoma domain area is affected by degradation (Strauss et al., 2013). Considering 10% of the area of the Yedoma domain covered with lakes and rivers, 4% covered with other deposits including deltaic and fluvial unfrozen sediments, this leaves 56% (780,000 km 2 ) of the Yedoma domain covered by frozen thermokarst deposits in drained thermokarst lakes ( Figure 2). Using these parameters, the overall mineral element stock is determined with Eq. 3 included into the bootstrapping calculation:

100
[ − ] × mineral concentration mgMin element kgdry soil 1, 000, 000, 000 Because Yedoma and Alas deposits have different properties for BD, ice content, and deposit thickness, mineral element stocks were estimated for each deposit type individually. For each bootstrapping step, the sample mineral element concentration and specific BD were paired (as they are not independent) and a specific stock was estimated based on one value for WIV and thickness. The stock was then multiplied by total coverage for total mineral stock estimation to the Yedoma domain scale. From these input parameters, multiple mineral element stocks were computed, one for each bootstrapping step. Eventually, from those multiple steps (n 10,000), a mineral element stock distribution was estimated from which the mean represents the best stock estimation of the considered mineral element. The error estimates in this study represent 2 standard deviations (mean ±2σ), with the assumption of a normal distribution. Computations were performed using R software (boot package; R Core Team, 2018). Supplementary information on input parameters (deposit thickness, WIV) are available in Supplementary Figure S2. Tables S1, S3). Linear regressions can be computed from this subset of data and parameters of the regression are estimated and based on a robust linear regression (alpha 0.95; Table 4). Among the 16 mineral elements measured, 10 elements (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr and Zr) display high correlations between pXRF and ICP-OES method with R 2 ranging from 0.725 (Al) to 0.996 (Ca ;  Table 4). Linear regression plots for these elements are presented in Supplementary Presentation S1. The other six elements (Mg, Ba, Cr, Cu, Ni and P) present weak (R 2 < 0.5) or no correlations between the two methods. According to these correlations, the mineral element stock quantification was performed on the 10 elements reliably measured by pXRF (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr and Zr), and the other six elements are not discussed further in this study.

Mineral element concentrations were measured with both the ICP-OES method (Total Elemental Analysis by ICP-OES Measurement After Alkaline Fusion) and the pXRF method (Total Elemental Analysis by ICP-OES Measurement After Alkaline Fusion) on a subset of samples (n 144) from Yedoma domain deposits (Supplementary
The uncertainty based on three identical samples on the ICP-OES measurement after alkaline fusion was lower than pXRF measurements, except for Ca, K, and Zr (Table 2). Nonetheless, the advantages of a non-destructive, rapid and cheaper method predominate for a large-scale mineral concentration assessment given that the coefficient of variation (2.1-8.7%) on the pXRF measurement was satisfying to differentiate between high and low concentrations values measured in these deposits (i.e., comparing the horizontal error bar with the total range of values represented on the X-axis in Supplementary Presentation S1).
In this dataset, the risk of overestimating the element concentration by pXRF in organic-rich samples (Shand and Wendler, 2014) is limited. Among the 1,292 Yedoma and Alas samples used in this study to build the YMCA dataset (Yedoma Domain Mineral Concentrations Assessment Dataset), less than 0.5% of the samples were characterized by TOC content similar to or higher than 40 wt% ( Table 5). Overestimation of the concentration was only observed for a single sample for Fe (Supplementary Presentation S1C) and for three samples for Ca (Supplementary Presentation S1D) out of the 144 samples from the subset (for samples with TOC content >40 wt%), and not observed for Si, Al, K, Ti, Mn, Zn, Sr and Zr.
Since the points corresponding to these organic-rich samples were excluded from the robust linear regressions (Supplementary Presentation S1), we assume that the matrix effect is not a source of significant bias to correct pXRF measurements of element concentrations using these regressions. Thus, the pXRF method can be applied for this first large-scale mineral element assessment of Yedoma domain deposits.

Yedoma Domain Mineral Concentrations Assessment Dataset
Robust linear regression equations with reliable R 2 ( Table 4) were applied for each pXRF-measured concentration to correct concentrations values for trueness (Figure 4 and Supplementary Presentation S1). The resulting YMCA dataset, freely available under doi:10.1594/PANGAEA.922724 (Monhonval et al., 2020), compiles 1,292 samples for the following elements: Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr and Zr.
In addition to the corrected mineral element concentrations, the YMCA dataset combines a wide range of relevant existing information on these 1,292 samples. Site and samples properties are integrated from AWI (Alfred Wegener Institute) Potsdam and Stockholm University data. We associated the lithology of the underlying bedrock, and the type of unconsolidated sediments at the surface characterizing each site from a spatial join using Arc Map 10.4. Specifically, coordinates from each site are joined based on their spatial location to Global Lithology Map (GLiM; Hartmann and Moosdorf, 2012) and Global Unconsolidated Sediments Map (GUM; Börker et al., 2018). Typically, we included the lithology of the underlying bedrock to evaluate its potential control on the mineral element concentrations in the deposits.
The attribute columns of this dataset are organized as follows: i) all site and sample properties (sample ID, description (identifying active layer (AL)), type of deposit, site location, number of profile, GPS coordinates, country, GLiM-lithology, GUM-unconsolidated sediment type, GUM-age, sample depth below surface level (b.s.l) or height above sea/river level (a.s.l), sediment characteristics, BD, gravimetric and absolute ice content, TOC values); and ii) corrected pXRF concentrations based on linear regressions from Linear Regression for Accurate Mineral Element Concentration Measurements. Table 5 summarizes corrected mineral element concentration values as well as BD, TOC and ice content from the whole Yedoma domain region, with no differentiation between Yedoma and Alas deposits samples (n 1,292).

Yedoma Domain Deposits Mineralogy
The main mineral phases (i.e., primary and secondary minerals) identified in selected Yedoma, Alas and fluvial deposits are presented in Table 6. The following minerals were identified in all bulk samples: quartz, feldspar plagioclase, micas and kaolinite. Chlorite was identified in Siberian deposits from Sobo Sise, Buor Khaya and Kytalyk. Calcite and dolomite were only detected in Alaskan Yedoma deposits from Colville and Itkillik. The mineralogy is generally similar along the profile depth for each location (Supplementary Figure S3). In these silty-fine sediments, the clay content ranges between 6 and 18% (silt content 67-79%; sand content 3-55%), i.e., well in line with Strauss et al. (2012). The diffractograms on clay fractions from Buor Khaya highlighted the presence of kaolinite, illite, and smectite (Supplementary Figure S3D-I). These data highlight the co-existence in these deposits of highly stable minerals such as quartz with more weatherable minerals such as dolomite, calcite, or feldspar plagioclase (Goldich, 1938;Wilson, 2004;Cornelis et al., 2011), and the presence of clay minerals characterized by higher specific surface area relative to the other mineral constituents (i.e., smectite, kaolinite, illite) (Kahle et al., 2004;Saidy et al., 2012).

Mineral Element Stocks in the Yedoma Domain
The quantification of mineral element stocks is a major step to assess the distribution of mineral elements in the Yedoma domain. We performed mineral elements stock estimations with the bootstrapping technique (Mineralogy by X-Ray Diffraction). For example, the mean Fe stock from Yedoma deposits (n 814) is 147 ± 43 Gt (±2σ) using the theoretical normal distribution (mean-bootstrapping technique; Figure 5).
Based on the YMCA dataset, the bootstrapping technique was applied to calculate the stocks in Yedoma deposits for the 10 elements (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr, and Zr), basing on the same input parameters for BD, thickness, WIV and coverage. The estimate mean and standard deviation derived the best estimation of the Yedoma deposit stocks for each specific mineral element. With a similar approach, mineral element stocks (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr, and Zr) were estimated in Alas deposits. The element stocks in the Yedoma domain were obtained by adding element stocks in Yedoma and Alas deposits ( Table 7). The results highlight that the mineral element with the largest stock in the Yedoma domain is Si (2,739 ± 986 Gt) followed by Al, Fe, K, Ca, Ti, Mn, Zr, Sr, and Zn ( Figure 6). In comparison, the estimated OC stock in the Yedoma domain reaches around 324-466 Gt ; Figure 6). Overall, Yedoma and Alas deposits ( Table 7) represent 86% of the total Yedoma domain area (Figure 2): this means that mineral element stocks are based on 86% of the deposits. The remaining surfaces of the Yedoma domain correspond to deltaic deposits (4%) or lakes and rivers (10%). Element stocks in deltaic deposits have not been estimated in this study due to the scarce number of available mineral concentration data in deltaic deposits (0.6% of the YMCA dataset) and the lack of empirical estimations of their thickness, or their ice volume (even if ice volume is probably negligible in these deposits). Waterbodies, lakes and rivers are not considered in this assessment focusing on the mineral element content in surface solid materials vulnerable to thaw or already thawed in the past. Sediments underlying lakes and rivers are therefore outside the scope of this evaluation and would only become relevant if drainage occurs. Absolute stock estimates (in Gt) allow direct comparison between Yedoma and Alas deposits, despite their different thickness (19.6 and 5.5 m, respectively), WIV and coverage. Mineral element density estimations (kg m −3 ) are available (Supplementary Table S5) in which thickness and WIV have been neglected: this estimation is the product of mineral element concentrations and BD using bootstrapping technique.

Portable X-Ray Fluorescence
The main advantage of the pXRF method is to ensure a nondestructive, more rapid (less steps involved; Figure 4) and cheaper method for mineral element concentration measurement compared to the ICP-OES measurements after alkaline fusion. This is beneficial when dealing with large-scale assessments involving thousands of archive samples from the Yedoma domain. The concentrations of 10 elements could be reliably measured in these deposits (i.e., Si, Al, Fe, K, Ca, Ti, Mn, Zr, Sr, and Zn). The main limitation of the pXRF method is that raw concentration data cannot be used before applying a linear regression to correct for systematic error. This drawback related to inaccurate raw pXRF measurements can be rectified by correcting the raw pXRF values with accurate values obtained with the ICP-OES method after alkaline fusion (calibration based on 144 samples in this study). This correction can only be applied when the correlation between pXRF and ICP-OES concentrations is satisfactory (here, R 2 > 0.7). In these deposits, the concentrations of some elements (Mg, P, Cu, Ni, Ba) could not be assessed due to poor pXRF-ICP-OES correlation (R 2 < 0.5). Linear regressions used to correct raw pXRF concentrations depend on internal geometry of the pXRF device used. This means that each pXRF device needs its own linear regression and that a single linear regression equation cannot be used with different pXRF devices. Moreover, pXRF measurements are also matrix-dependent. The matrix (i.e., organic content, bulk density) of a sample can affect TABLE 5 | Statistical summary of YMCA dataset (n 1,292) for bulk density (BD), total organic carbon content (TOC), absolute ice content, and corrected mineral element concentrations (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr and Zr) in Yedoma domain deposit samples (grouping Yedoma and Alas deposits together). Columns stand for N (number of data); Missing (missing data); MIN (minimum value reported); Q (quantiles); Q1 (first quantile); Q3 (third quantile); MEDIAN; MEAN; MAX (maximum value reported); SD (standard deviation); MAD (median absolute deviation standardised to conform with the normal distribution); IQR (interquartile range standardised to conform with the normal distribution); CV (coefficient of variation (%)); CVR (robust coefficient of variation (%)). pXRF-measured concentrations and therefore linear regressions must be calibrated with samples from the same matrix (here, we used a subset of 11% of samples from the Yedoma domain). For some other elements, pXRF measurements are not possible due to the low atomic mass of these elements (N, Na).

Mean-Bootstrapping Technique
The main advantage of the mean-bootstrapping technique is to ensure a reliable calculation of non-parametric uncertainty estimates. The bootstrapping technique allows to include variability in mineral element concentrations and BD, WIV and thickness of the deposits. However, standard deviations from stock estimations ( Table 7) are markedly large because of the conservative (observation-based) approach using all measurement for bootstrapping and deriving the mean afterwards. The only fixed parameter included in mineral element stock calculation (Eq. 3) is coverage. The coverage is set to 410,000 and 780,000 km 2 for Yedoma and Alas deposits, respectively. These estimations are subject to uncertainty but are the best current estimates, identical to the ones used in Strauss et al. (2013) to allow direct comparison with OC stocks. However, with the ongoing Arctic warming, these coverage estimations are dedicated to change with time. Because mineral element stocks are based on different WIV, thickness and coverage, the mineral element density expressed in kg m −3 allows to compare Yedoma and Alas deposits for an identical one cubic meter of sediments (without ice-wedges volume) (Supplementary Table S5).

IMPLICATIONS UPON PERMAFROST THAW Implications for the Evolution of Mineral-Organic Carbon Interactions
The YMCA dataset allows investigating the evolution of selected element to OC ratio upon thermokarst processes resulting from permafrost thaw, i.e., between Yedoma and Alas deposits. Organo-mineral associations with Fe-, Al-, Mn-oxides or clay  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 703304 minerals using polyvalent cations bridging between OC and mineral surfaces (e.g., Ca 2+ , Sr 2+ in neutral and alkaline soils or Fe 3+ and Al 3+ in acid soils), or organo-metallic complexes involving Fe 3+ , Fe 2+ and Al 3+ ions are well known to contribute to OC stabilization (Lutzow et al., 2006). Some mineral elements such as Fe or Mn are known for an inhibition effect on methane production (Lovley and Phillips, 1987;Beal et al., 2009;Herndon et al., 2015;Sowers et al., 2018), despite microbial function being considered as a larger part of the limitation for methane production (Monteux et al., 2020). The present dataset highlights a decrease of the Fe/TOC ratio and the Al/TOC ratio from Yedoma to Alas in Alaska, and virtually no change in these ratios in the Siberian Yedoma domain ( Figures 7A,B). Providing evidence for changes in Fe/TOC or Al/TOC ratios upon thermokarst processes in a portion of the Yedoma domain highlights the need for further investigations regarding the evolution of interactions between OC and Fe or Al upon thawing.
In addition, the clay minerals present in these Yedoma domain deposits (confirmed by XRD; Table 6) can be associated to the formation of aggregates, thereby contributing to spatial inaccessibility of OC for decomposer organisms due to occlusion in aggregates (Lutzow et al., 2006). Iron oxides are also known to contribute to the formation of aggregates (Eusterhues et al., 2005;Kleber et al., 2015). Changing Fe/ TOC ratio between Yedoma and Alas in Alaska ( Figures  7A,B) reinforces the need to further investigate the potential role of the evolution of aggregates for OC availability upon thawing (Monhonval et al., 2021). Changing conditions for mineral protection of OC upon thawing is likely to influence OC microbial degradation (Gentsch et al., 2015;Herndon et al., 2017;Kögel-Knabner et al., 2010;Opfergelt, 2020), thereby contributing to modulate the permafrost carbon feedback (Schuur et al., 2015).

Implications for Mineral Weathering and Nutrient Supply to Ecosystems
The YMCA dataset allows investigating the change in concentration of soluble elements such as Ca upon thermokarst processes resulting from permafrost thaw, i.e., between Yedoma and Alas deposits. As shown on a density plot ( Figure 7C), sediments from Alas deposits are characterized by lower Ca concentrations compared to Yedoma deposits. This was observed between Yedoma and Alas deposits in Alaska, and in Siberia ( Figure 7C).
Higher Ca concentrations in Yedoma and Alas deposits from Alaska relative to the deposits from Siberia can be explained by the local lithology, i.e., carbonate rocks from the northern Brooks Range contributing to the deposits in Alaska (Till et al., 2008). The lithology of the bedrock, beneath Quaternary deposits, inferred from the GLiM map, is similar between Yedoma and Alas deposits (Supplementary Figure  S4). This can likely be explained by the fact that the Yedoma and Alas deposits from the Yedoma domain mainly originate from the same source material on local to regional scales Strauss et al., 2017). Since most of the studied regions are covered by thick Quaternary deposits, mineral element concentrations are likely more influenced by the mixing of the unconsolidated sediments contributing to the Quaternary deposits rather than by the lithology of the underlying bedrock. The lithological similarity between FIGURE 6 | Mineral element stocks (in Gt) in the Yedoma domain deposits (including Yedoma and Alas deposits) for Si, Al, Fe, K, Ca, Ti, Mn, Zr, Sr, and Zn. The surface area of the circles is proportional to the stock of the element considered. The organic carbon stock in the Yedoma domain  is provided for comparison (the dotted line represents the lower range of the estimate, i.e., 327 Gt). Yedoma and Alas deposits can also be explained by the fact that Alas deposits are dominated by reworked sediments from former Yedoma deposits . Therefore, Ca depletion in Alas deposits relative to Yedoma deposits from one region is not lithology dependent but potentially results from mineral weathering and leaching processes of soluble elements such as Ca during former thawing periods. Indeed, Alas formation history includes lake formation and drainage, and the dynamism of such formation over the past thousands of years may lead to leaching processes of soluble elements, a commonly observed process in non-cryogenic soils (Stumm and Morgan, 1995) and cryogenic soils (Ping et al., 2005). The co-existence of more weatherable (dolomite, calcite, or feldspar plagioclase) and less weatherable (quartz) mineral constituents in Yedoma deposits ( Table 6) is likely to lead to an increasing contribution from mineral weathering to the inorganic carbon budget (Zolkos and Tank, 2020) upon thermokarst processes. The balance between carbonate or silicate mineral weathering is known to influence the shortterm and long-term atmospheric CO 2 consumption (Berner et al., 1983). Locally, the oxidation of sulfides such as pyrite present in low amount, not detected here but reported in beach deposits (Schirrmeister et al., 2010), can lead to the formation of sulfuric acid and the subsequent weathering of carbonate mineral without involving atmospheric CO 2 consumption (Zolkos and Tank, 2020).
The YMCA dataset provides total mineral element concentrations, including i) a mineral element reservoir known to be rapidly available to vegetation and microorganisms in permafrost soils (Ping et al., 1998(Ping et al., , 2005, and ii) a mineral element reservoir providing a longer term supply of nutrients through mineral weathering (Zolkos and Tank, 2020). Indeed, the weathering of minerals such as dolomite, calcite, or feldspar plagioclase is expected to release soluble cations such as Ca available to supply nutrients to terrestrial and aquatic Arctic ecosystems, or to interact with OC (Implications for the Evolution of Mineral-Organic Carbon Interactions). Among the 10 elements included in the YMCA dataset, macro-(e.g., K, Ca) and micro-nutrients (e.g., Fe, Mn, Zn) are required for plant nutrition and also regulate other vital processes for plants and microorganisms Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 703304 growth and metabolic activity (DalCorso et al., 2014). Briefly, K regulates vital processes, such as photosynthesis, water and nutrient transportation, or protein synthesis (Marschner, 2012). Ca is a major second messenger in plant signal transduction, mediating stress-and developmental processes (Liese and Romeis, 2013). Micro-elements are required as cofactor for some essential proteins (Fe; Morgan and Connolly, 2013), play an essential role as a photosynthetic function or in metallo-protein conformation (Mn; Yang et al., 2008) or can have enzymatic functions (Zn; Lindsay, 1972). Some non-essential elements, such as Si and Al, can stimulate plant growth by playing with abiotic-biotic stress resistance and symbiosis (Richmond and Sussman, 2003;DalCorso et al., 2014). In aquatic ecosystems, silicon is a limiting nutrient for diatoms and other siliceous organisms, thereby controlling diatom abundance and community structure in the ocean, and as a result, food web and CO 2 uptake by photosynthesis (Smetacek, 1999;Yool and Tyrrell, 2003). The leaching of the more mobile cations such as Ca or K can be estimated by comparing ratios between mobile and immobile elements: this is why assessing the concentration of elements considered as immobile such as Ti or Zr can be useful to evaluate the advance of weathering in a soil relative to its parent material, and thereby the soil mineral reserve (Kurtz et al., 2000;Hodson, 2002;Jiang et al., 2018). However, the reserve in important nutrients such as P, N, S, and Mg in the Yedoma domain could not be estimated with the method used in this study (Sect. 5). The YMCA dataset is a first step needed in order to evaluate the impact of widespread rapid permafrost thaw through thermokarst processes (Turetsky et al., 2020) on the mineral element concentrations in the deposits and the potential implications for OC and mineral nutrient supply.

CONCLUSION
This study presents a method combining portable X-ray fluorescence with a bootstrapping technique to generate the first mineral element inventory of permafrost deposits from the ice-rich Yedoma region, i.e., never thawed Yedoma deposits and previously thawed Alas deposits for a mean thickness of 19.6 and 5.5 m, respectively. In the resulting Yedoma domain Mineral Concentrations Assessment (YMCA) dataset, the total concentrations of 10 mineral elements (Si, Al, Fe, Ca, K, Ti, Mn, Zn, Sr, Zr) in Yedoma domain deposits have been quantified in 75 different profiles. The associated mineral elements stocks are shown to be in the same order of magnitude for Al and Fe than for OC, and to decrease from Si, Al, Fe, K, Ca, Ti, Mn, Zr, Sr, to Zn. This dataset allows tracking dynamic processes controlling mineral element concentrations in thawing environments, as illustrated by lower Ca concentration in Alas deposits relative to Yedoma deposits highlighting potential Ca leaching upon thawing. Providing the YMCA dataset is contributing to improve the knowledge on the mineral constituents in ice-rich permafrost, i.e., a necessary step to better understand the evolution of mineral-OC interactions, mineral weathering and mineral nutrient supply to ecosystems upon permafrost thaw. The YMCA dataset is particularly relevant given the increasing occurrence of abrupt thaw in ice-rich permafrost regions.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.pangaea. de/10.1594/PANGAEA.922724.

AUTHORS CONTRIBUTION
AM, SO, and JS conceived the project. JS, GG, LS, MF, and PK retrieved samples from Alaska and Siberia from many different field expeditions. AM did the pXRF measurements with the help of EM. AM, SO, EM, BP, and AV contributed to set up the YMCA database. AM analyzed the data and calculated to the stocks based on the code developed by JS for carbon stock estimation using mean-bootstrapping. AM wrote the manuscript with contributions from all co-authors.

ACKNOWLEDGMENTS
The authors acknowledge Hélène Dailly, Anne Iserentant, Elodie Devos and Laurence Monin from the MOCA analytical platform at UCLouvain for mineral elemental analysis. We thank Waldemar Schneider (AWI logistics) and Dmitry Melnichenko (Hydrobase Tiksi) for many years of logistical support, and Catherine Hirst, Maxime Thomas, and Pierre Delmelle for fruitful scientific discussions. We thank the two reviewers for their constructive comments.