Paleo±Dust: quantifying uncertainty in paleo-dust deposition across archive types

Mineral dust aerosol concentrations in the atmosphere varied greatly on glacial–interglacial timescales. The greatest changes in global dust activity occurred in response to changes in orbital parameters (which affect dust emission intensity through glacial activity) and the lifetime of dust in the atmosphere (caused by changes in the global hydrological cycle). Long-term changes in the surface dust deposition rate are registered in geological archives such as loess, peats, lakes, marine sediments, and ice. Data provided by these archives are crucial for guiding simulations of dust and for better understanding the natural global dust cycle. However, the methods employed to derive paleo-dust deposition rates differ markedly between archives and are subject to different sources of uncertainty. Here, we present Paleo±Dust, an updated compilation of bulk and <10 µm paleo-dust deposition rates with quantitative 1σ uncertainties that are inter-comparable among archive types. Paleo±Dust incorporates a total of 285 pre-industrial Holocene (pi-HOL) and 209 Last Glacial Maximum (LGM) dust flux constraints from studies published until December 2022, including, for the first time, peat records. We also recalculate previously published dust fluxes to exclude data from the last deglaciation and thus obtain more representative constraints for the last pre-industrial interglacial and glacial end-member climate states. Based on Paleo±Dust, the global LGM:pi-HOL ratio of <10 µm dust deposition rates is 3.1 ± 0.7 (1σ). We expect Paleo±Dust to be of use for future paleoclimate dust studies and simulations using Earth system models of high to intermediate complexity. Paleo±Dust is publicly accessible at https://doi.org/10.1594/PANGAEA.962969 (Cosentino et al., 2024).

1 Introduction

Mineral dust aerosols interact with climate through multiple mechanisms at different timescales, constituting a long-recognized, relevant component of the Earth system. However, although some of the dust–climate interactions are well understood in principle, many remain relatively unconstrained quantitatively, to the point where it is still unknown whether the net radiative effect of dust implies net cooling or warming in the present day (Kok et al., 2023). Thus, dust remains an important contributor to past and future climate change uncertainty (e.g., Andreae et al., 2005; Sherwood et al., 2020).

Both natural and anthropogenic processes are responsible for present-day dust emissions (Ginoux et al., 2012; Stanelle et al., 2014; Chen et al., 2018). To constrain naturally induced changes in the dust cycle (a prerequisite to disentangle the present-day, human-induced component), one can turn to pre-industrial records of dust activity. Based on such records spanning the last several hundreds of thousands of years, it has been shown that the most dramatic changes in dust activity have taken place at glacial–interglacial timescales (Lambert et al., 2008; Simonsen et al., 2019; Struve et al., 2022), although significant variability is present at all timescales (Lovejoy and Lambert, 2019). Thus, global observational datasets of past dust activity for the pre-industrial Holocene and last glacial periods are key to constraining the natural processes involved in dust–climate interactions through multiproxy comparisons, their use in guiding Earth system simulations of past climates, and as a tool for model–data comparisons.

Paleo-dust archives include continental (i.e., loess-paleosol sequences, lake, peat, and ice) and marine sediments. The dust dynamics parameter that may be calculated from information preserved in these deposits is the dust deposition rate, from which other aspects of the dust cycle must be inferred, such as atmospheric dust concentrations and the dust emission rate, typically through Earth system simulations. Each archive type has its strengths and weaknesses in terms of spatial representativity and how well they preserve dust (Albani et al., 2015), and it is only through a combination of constraints from all these archives that a spatially extensive coverage of paleo-dust activity can be achieved. The degree of coverage is important due to the presence of short-scale variability in dust emissions associated with the uneven distribution of dust sources as well as spatially abrupt changes in dust deposition fluxes associated with precipitation-controlled wet deposition and with dry deposition controlled by the distances to dust sources.

Several compilations of paleo-dust deposition fluxes exist that combine constraints from different archive types (Mahowald et al., 1999, 2006; Kohfeld and Harrison, 2001; Tegen et al., 2002; Maher et al., 2010; Albani et al., 2014, 2015; Lambert et al., 2015) or specialize in marine sediment cores (Kohfeld et al., 2013; Kienast et al., 2016). These compilations have allowed data–model comparisons (Hopcroft et al., 2015; Kienast et al., 2016; Ohgaito et al., 2018; Otto-Bliesner et al., 2020; Braconnot et al., 2021; Krätschmer et al., 2022), the tuning of modeled dust emission rates against observed dust deposition rates (Mahowald et al., 2006; Albani et al., 2014, 2015, 2016; Albani and Mahowald, 2019; Braconnot et al., 2021), and simulations with intermediate-complexity biogeochemistry models (Lambert et al., 2015, 2021; Heinemann et al., 2019). One key aspect missing from these compilations, except the Holocene compilation by Albani et al. (2015), is a quantification of the uncertainty in paleo-dust deposition rates. This is essential for compilations derived from multiple archive types, as each of these is subject to different sources of uncertainty. Explicit site-specific uncertainties in paleo-dust deposition rates may be used for selecting subsets of observations against which to compare dust simulations, weighting observations for dust emission tuning purposes, or deriving distributions of global interpolations, for example through bootstrapping, Monte Carlo experiments, or Bayesian approaches.

During the 7 years since the last paleo-dust compilation was published (that of Kienast et al., 2016), there has been significant progress in the global coverage of paleo-dust proxies. Here, we update the previous compilations of global paleo-dust deposition flux and grain size observations with the newest data. Moreover, three improvements are introduced. First, given the expected high variability of dust deposition fluxes during the transition between the last glacial period and the current interglacial (compared to that within any of those periods), we exclude any data corresponding to the last deglaciation to more precisely quantify dust deposition fluxes during the last glacial and current interglacial end-member climate states. Second, we derive site-specific, quantitative uncertainties of bulk and <10 µm paleo-dust deposition fluxes for both Holocene and last glacial observations. Third, we include, for the first time, dust deposition rate observations from peat bogs.

The manuscript is organized as follows: Sect. 2 details our methodology in constructing Paleo±Dust, Sect. 3 presents the main results, Sect. 4 describes the structure of the datasets, and Sect. 5 presents some concluding remarks.

2 Methods

2.1 Dust mass accumulation rate

The mean dust mass accumulation rate (DMAR) for a given time window at a location on the surface of the Earth may be quantified as follows by dating two horizons along a vertical sedimentary profile (t_top,t_bottom), measuring the profile thickness between these two horizons (h_thick), the mean dry bulk density (DBD) between these two horizons, and the mass fraction represented by atmospherically derived mineral dust not associated with direct volcanic ash fall (i.e., dust), henceforth defined as the eolian content (EC):

$$\mathrm{DMAR}={\displaystyle \frac{h_{\mathrm{thick}}\times \mathrm{DBD}\times \mathrm{EC}}{\mathrm{\Delta }t}},$$

where $\mathrm{\Delta }t=t_{\mathrm{bottom}}-t_{\mathrm{top}}$. Given that the bulk DMAR may be very sensitive to small, local dust sources, we also define a finer-grained DMAR, which is more representative of deposition of far-traveled dust derived from the main dust sources globally. A 10 µm dust diameter is used to define this finer-grained dust fraction, as this is either the uppermost- or intermediate-bin threshold diameter that is found in the greatest number of CMIP6 Earth system models with dust representation (Zhao et al., 2022), thus facilitating model–data comparisons. It is also the recommended definition for the threshold between the coarse and super-coarse size ranges of dust (Adebiyi et al., 2023). The fraction of DMAR that corresponds to dust particles less than 10 µm in diameter (DMAR₁₀) is simply

$${\mathrm{DMAR}}_{\mathrm{10}}=\mathrm{DMAR}\times f_{\mathrm{10}},$$

where f₁₀ is the <10 µm dust mass fraction. This is the approach to calculating paleo-dust deposition fluxes used for all continental archives (i.e., loess-paleosol sequences, lakes, peat, and ice) and for a few sediment cores in the marine realm. Most DMAR constraints from marine sediment cores use the ²³⁰Th normalization technique, which is, for the most part, independent of the profile chronology, as discussed further in Sect. 2.8.

The absolute error of DMAR (σDMAR) and DMAR₁₀ (σDMAR₁₀) can be expressed as

$$\mathit{\sigma }\mathrm{DMAR}=\mathrm{DMAR}\sqrt{\begin{array}{c}{\left({\displaystyle \frac{\mathit{\sigma }{h}_{\mathrm{thick}}}{h_{\mathrm{thick}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\sigma }\mathrm{DBD}}{\mathrm{DBD}}}\right)}^{\mathrm{2}}\\ +{\left({\displaystyle \frac{\mathit{\sigma }\mathrm{EC}}{\mathrm{EC}}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathit{\sigma }{t}_{\mathrm{bottom}}^{\mathrm{2}}+\mathit{\sigma }{t}_{\mathrm{top}}^{\mathrm{2}}}{{\left(t_{\mathrm{bottom}}-t_{\mathrm{top}}\right)}^{\mathrm{2}}}}\end{array}}$$

and

$$\mathit{\sigma }{\mathrm{DMAR}}_{\mathrm{10}}={\mathrm{DMAR}}_{\mathrm{10}}\sqrt{{\left({\displaystyle \frac{\mathit{\sigma }\mathrm{DMAR}}{\mathrm{DMAR}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\sigma }{f}_{\mathrm{10}}}{f_{\mathrm{10}}}}\right)}^{\mathrm{2}}},$$

where σh_thick, σDBD, σEC, σt_bottom, σt_top, and σf₁₀ are the absolute errors of h_thick, DBD, EC, t_bottom, t_top, and f₁₀, respectively. All errors reported in Paleo±Dust are Gaussian 1σ uncertainties.

Our approach to assigning uncertainties to the components of DMAR combines objective and subjective considerations. On the one hand, when measurement uncertainties are reported in the original studies, these uncertainties are used. One example is the use of ²³²Th as a dust proxy to calculate EC in marine sediments. This method requires the normalization of measured ²³²Th in marine sediments to the mean global concentration of ²³²Th in dust: the variability of this normalizing value provides a means to calculate the uncertainty in EC. Instead, when uncertainties are not reported, these are defined based on the distribution of reported relative uncertainties for sites of the same archive type (conservatively choosing the 75th percentile of this distribution). On the other hand, when the uncertainty of any component of DMAR is not reported, it is assigned. Assigning uncertainty is more problematic, as the true uncertainty is usually not known. This is the case for example for uncertainties in chronologies of loess-paleosol sequences derived through pedostratigraphy. In such cases, we assign higher subjective uncertainties to methods deemed more uncertain (based on the literature and the authors' experience). This combination of approaches implies that reported DMAR uncertainties in Paleo±Dust should be considered to reflect the relative uncertainties among sites (more so than the exact uncertainties for each site). Sections 2.3–2.8 expand on these considerations for each archive type and DMAR component.

2.2 Criteria for the inclusion of sites

Previous compilations of paleo-DMAR served as the starting point for the construction of Paleo±Dust. A literature search was performed to include new sites for studies published until (and including) December 2022. Recent advances in the study of peat bogs as a dust archive (e.g., De Vleeschouwer et al., 2014; Kylander et al., 2016, 2018; Marx et al., 2018; Sjöström et al., 2020) have allowed the inclusion of this archive type in a dust compilation for the first time.

2.2.1 Time span

One of the main objectives of Paleo±Dust is to compile global constraints on DMAR for the ongoing interglacial and last glacial periods. In particular, Paleo±Dust compiles observations for the pre-industrial Holocene (pi-HOL; between the year 1850 CE and 11.7 ka) and the Last Glacial Maximum (LGM; between 19.0 and 26.5 ka) (Clark et al., 2009). We thus exclude observations of DMAR for the present day and recent past as well as for the last deglaciation (11.7–19.0 ka) to represent DMAR due to natural processes for the two main quasi-equilibrium climate states at glacial–interglacial timescales. Whether this may be achieved or not for any particular site will depend on the sampling resolution along the vertical profile for dating. In order to include a site in Paleo±Dust, it needs to be possible to define a sub-section of the full sampled vertical section of a sedimentological archive whose time interval falls mostly (i.e., ≥75 %) within either the pi-HOL or LGM as previously defined. Another criterion that needs to be fulfilled is that the site's time window is resolved within uncertainty: if $t_{\mathrm{bottom}}\le t_{\mathrm{top}}+\mathit{\sigma }{t}_{\mathrm{top}}$ and $t_{\mathrm{top}}\ge t_{\mathrm{bottom}}-\mathit{\sigma }{t}_{\mathrm{bottom}}$, then that DMAR constraint is discarded.

Each entry in Paleo±Dust is associated with a specific time window within pi-HOL or LGM defined by t_bottom and t_top. This allows the selection of subsets of data that include or exclude a specific time range of interest. As most of the previous compilations only reported general time periods (e.g., Holocene, Late Holocene, LGM, last glacial), the original references of each entry in these compilations were revisited to extract exact values for t_bottom and t_top. Given the imposed age criteria described above, most of the reported DMAR values in these compilations are different from those reported in Paleo±Dust for the same site, as a different set of samples from each site may have been used in this study to calculate DMAR.

2.2.2 Geomorphological setting

The most important consideration when selecting sites for Paleo±Dust is that the lithogenic fraction of the sediment should be dominated by atmospheric deposition or that it is possible to estimate the atmospherically derived lithogenic fraction quantitatively. The site's geomorphological setting largely controls this. Sites of loess deposits with a massive structure reflective of the deposition of dust from the atmosphere are included in the database, while loessoid sites consisting of loess reworked by either fluvial or slope processes are excluded (e.g., Pye, 1995). Peat usually consists of an upper ombrotrophic section where lithogenic materials are exclusively supplied by the atmosphere and a lower minerotrophic section where lithogenic materials are also supplied laterally by groundwater (e.g., Shotyk, 1996). Only the ombrotrophic sections of peat sections are considered in this study. In the case of lake sediments, non-eolian lithogenic material may be advected into a studied site through fluvial inlets or via slope processes within the lake basin, a process known as lake sediment focusing (e.g., Blais and Kalff, 1995). We only included lake sites from endorheic basins, and we discarded sites that lack a quantification of the fraction of lithogenic material advected through lake sediment focusing. With respect to marine sediments, we excluded sites located in continental margins, as they may be influenced by riverine inputs, as well as high-latitude sites potentially influenced by ice-rafted debris (e.g., Kienast et al., 2016). In both cases, exceptions were made for sites where these non-eolian lithogenic inputs are quantified. Other marine sites may be affected by sediment focusing; that is, the lateral submarine transport of non-eolian sediment. The ²³⁰Th normalization method can be used to isolate the eolian lithogenic component in these cases (e.g., Francois et al., 2004). Finally, polar ice cores are retrieved from upland landscape positions, and, for those sites that are located far from local lithogenic sources (those included in Paleo±Dust), their lithogenic materials are assumed to be wind blown (i.e., Albani et al., 2015).

For loess, the local geomorphology and topography can enhance the dust deposition rate by focusing wind-blown material onto the surfaces of windward escarpments (e.g., Comola et al., 2019), which may translate into a higher local loess DMAR than the mean regional value (e.g., Xiong et al., 2015). This effect cannot be easily identified in the profile, particularly when a regional loess stratigraphy against which to compare the local stratigraphy is not available in the literature. In general, upland loess sites are considered to be less prone to topography-induced enhancements in mass accumulation rates (Kohfeld and Harrison, 2003). We therefore performed a case-by-case evaluation of the geomorphological setting of each loess site and discarded sites only when this potential problem was suggested explicitly in the original studies.

Conversely, the preservation potential of loess sites located in topographic lows within high-relief environments may be limited due to slope or fluvial erosional processes that disrupt the continuity of the loess stratigraphy, generating hiatuses. These erosional hiatuses may be evident in the field or may only be discerned through high-resolution dating (e.g., Stevens et al., 2018; Volvakh et al., 2022). For these reasons, we restricted the analyzed vertical extents of loess-paleosol profiles to sections with no intervening erosional hiatuses. We also completely excluded sites that have a loess stratigraphy that considerably deviates from a well-defined regional loess stratigraphy in regions where such a reference regional loess stratigraphy is available (based on the original authors' assessment).

2.3 Uncertainty in top and bottom ages

Direct absolute dating of the sedimentary units of interest is the main tool for defining a chronology for dust records. The two most popular techniques are ¹⁴C and optically stimulated luminescence (OSL) dating. When t_bottom and t_top were obtained by direct absolute dating, the reported 1σ uncertainties are used (Fig. 1). If the dating uncertainty was not reported, relative uncertainties of 3.4 %, 9.1 %, and 13.1 % are assigned to ¹⁴C-, OSL-, and thermoluminescence-based age determinations. These values represent the 75th percentile of the distribution of the relative uncertainty in the age measured via ¹⁴C (N=83), OSL (N=129), and thermoluminescence (N=20), respectively, for the sites included in Paleo±Dust. Absolute ages are usually reported as years before present, where the reference age (0 yr BP) may be the year of sampling or a fixed year (e.g., 1950 CE). We do not homogenize reported ages across studies so that they are referenced to the same year. Given that we focus on mean LGM and pi-HOL dust deposition rates, this does not introduce a significant uncertainty.

Figure 1 Criteria for assigning uncertainty to the top and bottom ages. OSL: optically stimulated luminescence, TL: thermoluminescence, AICC2012: Antarctic Ice Core Chronology 2012 (Veres et al., 2013), GICC05: Greenland Ice Core Chronology 2005 (Svensson et al., 2008).

If a study presents a continuous age model and this model includes, for example, the full LGM age interval (26.5–19.0 ka), then t_bottom is defined as 26.5 ka and t_top is defined as 19 ka (Fig. 1). If the study reports the modeled errors in these interpolated ages, we use these reported values. If it does not, we calculate the errors in t_bottom and t_top as the L2 norm of the errors of the two closest bracketing measured ages. The advantage of this error is that it is higher than the error of either measured age. This extra uncertainty can be thought of as being due to the interpolation.

When no absolute ages are available and the chronology is defined based on the correlation of a sediment parameter with any given reference record – for example, the correlation of the magnetic susceptibility signal in marine sediment cores or loess-paleosol profiles to the oxygen isotope record of marine benthic foraminifera (e.g., Lisiecki and Raymo, 2005) – then a considerably higher relative error is assigned (Kohfeld and Harrison, 2003; Fig. 1). Wiers et al. (2019) quantified the uncertainty in the chronology of a Late Pleistocene Arctic marine sediment core obtained by correlating magnetic properties in the core with global patterns of δ¹⁸O in benthic foraminifera. They found a high absolute error during the ongoing interglacial and last glacial periods of approximately ± 6 kyr, which remained relatively constant during this time span. We apply this same absolute error to t_top and/or t_bottom obtained through such correlations (Fig. 1). When these chronologies are also supported by absolute dating, we instead assume a considerably lower absolute error of ± 3 kyr. In the case of loess, chronologies may also be defined based on pedostratigraphy, by assuming that loess and paleosol units can be correlated to glacial and interglacial periods, respectively, as defined by marine isotope stratigraphy (Kohfeld and Harrison, 2003). Only pi-HOL dust deposition rates were obtained through this method (N=11, all from China; Kohfeld and Harrison, 2003), as the LGM period cannot be discriminated from the encompassing last glacial period in terms of any pedostratigraphic unit. An evaluation of the uncertainty in the chronologies of these sites would require validation against absolute dates, which are not available for these sites. Moreover, to the best of our knowledge, there are no similar validation studies for other sites that we could use as a model. We thus compensate for this lack of validation by assigning a higher absolute uncertainty compared to magnetic correlations (±9 kyr).

For the specific case where the surface of a sedimentary unit was one of the bounding surfaces for a DMAR constraint and no continuous age model was available, t_top is assumed to be 0 ka, and the uncertainty in this assumption is related to the difficulty in defining the surface level and quantified in the following way (Fig. 1): if, for example, t_bottom is dated at 1.0 ± 10 % ka at 15 cm depth, then, in 1 cm centered at the top surface, the time span covered is 0.067 kyr, assuming a linear relationship between age and depth. Then the age of the surface is fully contained within 0 ka ± (0.067/2) ka. Translating this absolute time interval to an equivalent 1σ uncertainty interval, we multiply by 68.27 % and get 0 ka ± (0.067/2) × 0.6827 ka. Finally, we multiply this 1σ interval by 1.1, which corresponds to the relative age uncertainty for t_bottom in our example (=10 %). The final surface age is 0.000 ka ± 0.025 ka (1σ). This procedure includes the uncertainty due to the difficulty in defining the surface level and the uncertainty due to the dating method.

Finally, all polar ice core DMAR constraints in Paleo±Dust were obtained based on correlation to a common chronology for each of the polar regions: the 2012 Antarctic Ice Core Chronology (AICC2012, Veres et al., 2013) for Antarctic ice cores and the 2005 Greenland Ice Core Chronology (GICC05, Svensson et al., 2008) for Greenland ice cores. In both cases, uncertainties in ages are as reported for each chronologic framework.

2.4 Uncertainty in sediment profile thickness

As $h_{\mathrm{thick}}=h_{\mathrm{bottom}}-h_{\mathrm{top}}$, where h_bottom and h_top are the measured depths of the bottom and top layers, respectively, the error of the sedimentary thickness between dated layers is

$$\mathit{\sigma }{h}_{\mathrm{thick}}=\sqrt{\mathit{\sigma }{h_{\mathrm{top}}}^{\mathrm{2}}+\mathit{\sigma }{h_{\mathrm{bottom}}}^{\mathrm{2}}}=\mathit{\sigma }h\sqrt{\mathrm{2}},$$

where $\mathit{\sigma }h=\mathit{\sigma }{h}_{\mathrm{top}}=\mathit{\sigma }{h}_{\mathrm{bottom}}$ is the error of the measured depth to a dated layer beneath the surface level.

The uncertainty in measuring layer depth is associated with the fact that the sample obtained for dating has a finite vertical height, defined by the amount of sample mass required to perform dating. If, for example, loess sampling for dating is carried out by inserting horizontal corers 2 cm in diameter into the vertical face of the sedimentary profile, then the full sample within ±1 cm of the center of the corer is included. We may translate this total depth range into an equivalent range with 1σ uncertainty by multiplying by 0.6827, resulting in an absolute 1σ uncertainty in depth (σh) of ±0.6827 cm. If not reported, we assumed a horizontal corer diameter of 7.9 cm, which is the mean value reported for loess sites in Paleo±Dust. For lake, marine, peat, and polar ice sampling where coring is performed vertically from the surface, the relevant quantity from which to calculate σh is the height of each sample in the vertical direction along the corer, or, in other words, the thickness obtained after field operations in the laboratory. If not reported, we assumed a sample height of 2 cm for marine sediment and peat bog records and 7.9 cm for lake sediment samples (the mean reported value for sample height in Paleo±Dust in all cases). All ice core sites in Paleo±Dust have a reported sample height.

When the site chronology is based on magnetic susceptibility, σh is calculated as $\mathit{\sigma }h=\mathrm{0.6827}\times h_{\mathrm{MS}}/\mathrm{2}$, where h_MS is the depth interval between susceptibility measurements.

2.5 Uncertainty in dry bulk density

In a study of forest soils in California, it was found that between three and 17 samples were needed to estimate the mean soil DBD to ±10 % at a 95 % confidence level (Han et al., 2016). Here, a significantly higher DBD relative uncertainty of 15 % (for 1σ) is assumed for sites with no reported DBD uncertainty, given that several sites have less than three DBD determinations, and a greater variety of sediment and soil types are considered in this compilation compared to the study by Han et al. (2016). If, instead, the measurement error is reported in the original study, then this value is used. Finally, if DBD was measured but could only be estimated visually from a figure in the original study, a greater uncertainty of 20 % is used. The same uncertainty is assigned when the DBD was assumed equal to that at a site near (<100 km away) to where it was measured.

If DBD was not measured and no nearby sites with measured densities are available, a DBD value of 1.45 g cm⁻³ and an uncertainty of 30 % are used for loess and lake sites, except for loess sites in the Chinese Loess Plateau and in East Central Europe (west of 21.5^∘ E), where mean values of 1.48 g cm⁻³ (Kohfeld and Harrison, 2003) and of 1.497 g cm⁻³ (Újvári et al., 2010; Peric et al., 2020) are preferred, respectively. For peat bog sites, when DBD is not measured (N=5), the mean value of the DBD among all sites with DBD measurements globally is used (=0.12 g cm⁻³, N=18). No assumption about the DBD of marine sediment cores is required, as DBD was always measured at sites that require the DBD for the derivation of DMAR (those not based on ²³⁰Th normalization).

For polar ice cores, the density is taken to be that of ice: 916,750 g m⁻³. This value is assumed to have a small variation of ±50 g m⁻³, which translates into a 1σ variability of ±50 g m${}^{-\mathrm{3}}\times \mathrm{0.6827}=\mathrm{34.14}$ g m⁻³. This approach is justified by the fact that the density of ice remains almost constant with depth below approximately 100 m depth from the surface (e.g., Gerland et al., 1999).

2.6 Uncertainty in the mass fraction of dust

The calculation of EC and its uncertainty depends strongly on the type of dust archive considered (Fig. 2).

Figure 2 Criteria for assigning uncertainty to the non-volcanic dust fraction. OM: organic matter, bSiO₂: biogenic silica, PCA: principal component analysis.

2.6.1 Loess-paleosol sequences

In loess studies with a focus on dust dynamics, EC is usually assumed to be 1; that is, loess is assumed to be fully composed of eolian dust. However, while organic particles present in dust sources may be transported by wind and later deposited in the same manner as lithic particles (Muhs et al., 2014), here we assume that organic matter in loess is post-depositional in origin (e.g., Hatté et al., 2001), and we subtract the total organic carbon (TOC) to calculate EC. Instead, in some regions, carbonates constitute a significant fraction of airborne dust (Scheuvens and Kandler, 2014), which supports the interpretation that carbonates in loess are dominantly primary; that is, they are derived from dust sources (e.g., Li et al., 2013; Meng et al., 2019). However, carbonates in loess have also been shown to be authigenic (e.g., Da et al., 2023). Unfortunately, for most published studies it is not possible to calculate the contributions of primary and authigenic carbonates to the total carbonate content (TCC) of loess. We thus assume that 50 % of the carbonate present in loess is primary while 50 % is authigenic. Furthermore, previous research has shown that wet climates during interglacials tend to deplete 10 % of the carbonates in soils (by weight), while no loss occurs during glacial periods (Meng et al., 2015, 2018; Zhang et al., 2023). We account for this post-depositional loss of carbonates in Holocene loess-paleosol sections in our EC calculations.

The relative 1σ uncertainty in EC is assumed to be 10 % (20 %) if both (either) TOC and (or) TCC are (is) reported (Fig. 2). If both TOC and TCC are available and there is also a quantification of volcanic inputs, then the relative uncertainty is reduced to 1 %. Instead, when neither TOC nor TCC data are available, loess-paleosol units are classified as either organic carbon rich or organic carbon poor based on the physical description of the unit of interest, and the sum of TOC and TCC is assigned a value of 6 wt % or 2 wt %, respectively. These correspond, respectively, to the first and third quartiles of the sum of TOC and TCC for sites in Paleo±Dust where both TOC and TCC were determined (N=28). In addition, when TOC is not determined, values of 0.2 wt % and 1 wt % are assumed for LGM and pi-HOL sites, based on previous studies (Yang et al., 2015). In these cases where both TOC and TCC are assumed, the relative uncertainty in EC is highest at 30 %.

2.6.2 Lake sediment cores

The procedure for calculating EC in lake sediment cores is similar to that for loess deposits, except that biogenic silica is an extra sediment component that can be very relevant in terms of mass and needs to be corrected for (Fig. 2). Another effect to consider in lake sediment cores is sediment focusing, by which local sediments from the lake catchment can contribute to the siliciclastic mass. Some lake sediment studies isolate long-range dust from local catchment sediment, for example, through grain size end-member modeling, trace element geochemistry (Petherick et al., 2009), or corrections that consider lake bathymetry (Arcusa et al., 2020). Given the greater number of corrections required to calculate EC, the maximum potential relative uncertainty in EC is greater than that for loess records (i.e., 50 % versus 30 %).

2.6.3 Marine sediment cores

A few studies of dust using marine sediment records calculate EC based on the same principles as used with loess and lake sediment cores: by subtracting the carbonate, organic matter, and biogenic silica mass fractions and, when available, correcting for volcanic inputs (Fig. 2). However, most studies that look at dust archived in marine sediment cores use isotope ²³²Th measurements to calculate EC, assuming a global mean concentration of ²³²Th in dust of 14.0 ± 4.6 ppm (Kienast et al., 2016; Ouyang et al., 2022). The relative 1σ uncertainty in this mean value (i.e., 33 %) is used as the uncertainty in EC when calculated based on ²³²Th. Studies that used a different global mean value of ²³²Th concentration in dust are recalculated with this value.

When cores are extracted close to the continents, a significant fraction of the siliciclastic mass may be due to riverine inputs. Following Singh et al. (2011) and Kienast et al. (2016), we excluded sites located less than 300 km from the coast except for the equatorial Atlantic Ocean off the coast of Brazil, for which this distance is 600 km (Holocene) and 800 km (LGM). This rule does not hold for studies that isolate dust from riverine contributions to terrigenous sediment (e.g., McGee et al., 2013).

Another potential terrigenous component that may obscure dust signals in marine sediment cores is ice-rafted debris. Here, we follow Kienast et al. (2016) and exclude marine sediment cores that are potentially affected by ice-rafted debris with high probability; that is, sites located poleward of 55^∘ N in the North Pacific, poleward of 50^∘ N in the North Atlantic, and poleward of 50^∘ S in the southern oceans. These same latitudinal thresholds are applied to both pi-HOL and LGM paleo-dust sites. While the polar fronts may have changed latitudinally over time (and thus the extent of the influence of ice-rafted debris may have too), there is no evidence of this from the relationship between latitude and paleo-lithogenic fluxes in marine sediment cores (Kienast et al., 2016). Sites located poleward of these latitudinal thresholds are kept when lithogenic fluxes are confirmed to be mostly due to eolian dust, for example, through the measurement of n-alkanes, a proxy for continent-derived materials (e.g., Lamy et al., 2014).

2.6.4 Polar ice cores

The value of EC was calculated for the two Antarctic ice cores included in Paleo±Dust (for both pi-HOL and LGM) based on Coulter counter insoluble particle volume concentration data, assuming an insoluble particle density of 2.5 g cm⁻³ (Delmonte et al., 2005). The relative 1σ uncertainty in Antarctic EC is 15.3 %, which represents the L2 norm of a 11.4 % error component due to uncertainty in the insoluble particle density (2.0–2.8 g cm⁻³, Kohfeld and Harrison, 2001), a 10 % error component due to the lack of a volcanic input correction, and a 2 % component derived from replicate Coulter counter measurements (Fig. 2).

In the case of Greenland ice cores, EC is based on δ¹⁸O and Ca²⁺ concentration data. Sample depths with δ¹⁸O <−40 ‰ (cold periods) are assigned a Ca-to-dust ratio of 0.095, while sample depths with δ¹⁸O >−37 ‰ are assigned a Ca-to-dust ratio of 0.26 (Steffensen, 1997; Ruth et al., 2002), with a linear interpolation for in-between values of δ¹⁸O (Albani et al., 2015). The uncertainty in EC is assumed to be 22.4 %, which is the L2 norm of a 10 % uncertainty component due to possible volcanic inputs (just as for the Antarctic cores) and 20 % due to a combination of analytical and Ca²⁺ proxy uncertainties (Fig. 2; Albani et al., 2015).

2.6.5 Peat cores

Multiple ways to calculate EC in peat bogs are reported in the literature. One way is to measure the fraction of non-combustible mineral ash in total dry peat mass, which we can approximate as the siliciclastic mass fraction (e.g., Martínez Cortizas et al., 2020). A preferred approach is to use the concentration of one (e.g., Sharifi et al., 2018) or multiple (e.g., Pratte et al., 2020) conservative lithogenic elements that are not subject to post-depositional mobilization and have no anthropogenic source, such as Sc, Zr, Y, and the rare earth elements (REEs), typically normalized to the mean concentration of these elements in the upper continental crust. Irrespective of the method employed, the lowest relative uncertainty in the calculation of EC is assigned to cases where the choice of the dust proxy or combination of dust proxies is based on a principal component analysis (PCA, Fig. 2). This is because EC is sensitive to the choice of the geochemical proxy for dust (Shotyk et al., 2002; Eléonore Resongles, personal communication, 2023), and a PCA analysis usually provides the best way to identify the element or set of elements that best represents atmospherically derived lithogenic fluxes. However, high discrepancies in DMAR calculated with different elements that were all found to be associated with dust based on PCA analyses have also been reported (Kylander et al., 2016). This stresses the need to carefully select the elements to be used for DMAR calculations based on the site-specific element behavior. Specifically, it is recommended that when multiple elements are identified as the optimal dust proxy, the DMAR should be calculated based on the combination of those elements (e.g., based on the sum of the REEs, Sc, Y, and Zr if all these elements are identified as being equally suitable to represent dust).

To be consistent with the choice of uncertainties in EC for loess-paleosol sequences, lakes, and marine sediments, a relative uncertainty of 10 % is assigned to PCA-supported EC calculations for peat bogs (Fig. 2). If no PCA analysis is available, DMAR calculated based on multiple elemental proxies is preferred over a single-proxy approach, including the case where the single proxy is the non-combustible ash fraction. Moreover, when independent information is available that allows the identification of samples with high inputs of direct volcanic ash fall, such as that based on Nd isotopes (e.g., Vanneste et al., 2015), a correction is applied in which the high DMAR values for those samples are replaced by background DMAR values in the same profile.

2.7 Uncertainty in the <10 µm dust mass fraction

We assign the lowest uncertainty to f₁₀ when it is calculated from a full volumetric grain size distribution (Fig. 3). This uncertainty is associated with the reproducibility of measurements (∼ 0.7 % after the propagation of each bin's uncertainty) based on laser-diffraction determinations of 135 South American loess samples from three sites dated to between 8–53 ka (Coppo et al., 2022a), as well as with the use of different laser-diffraction devices and assumptions about their optical settings (19.9 %) based on measurements of loess samples from East Central Europe (Varga et al., 2019). We thus assign a 20 % relative uncertainty to f₁₀ when it is calculated from full grain size distributions using laser diffraction devices. Coulter counters measure particle volume more accurately than laser diffraction devices (Simonsen et al., 2018), so volumetric grain size distributions and f₁₀ values that are measured by this technique are arguably more accurate as well. We thus assign a lower relative uncertainty of 5 % to f₁₀ when it is calculated from Coulter counter measurements. This specific value of 5 % is, however, arbitrary, as the authors are not aware of studies that quantified sources of uncertainty of Coulter counter-derived grain sizes. In many cases where such measurements were carried out, the full-size distribution data were not published and f₁₀ was not reported in the original study. Albani et al. (2015) compiled these distributions for a number of Holocene sites published previously to 2015. Also, Albani et al. (2014) compiled f₁₀ values for Holocene and LGM dust archive sites published prior to 2014. For all these sites, f₁₀ is retrieved from these studies. For the rest of the sites in Paleo±Dust, if data are not available but the grain size distribution is plotted in the original study, we estimated f₁₀ visually, adding an extra 10 % of uncertainty (30 % in total). To perform this visual estimation, we used a vector graphics editor (Adobe Illustrator) to draw two polygons: one that encompassed the area under the curve of the volumetric grain size distribution for all measured particle sizes and one for sizes <10 µm. We calculated f₁₀ as the ratio of the latter to the former area. Greater uncertainty (30 %) is associated with f₁₀ when it is calculated from reported grain size bin volumetric abundances (e.g., clay, silt) because either grain size was determined using the sieve and pipette method or the full grain size distribution is not reported. Also, a higher uncertainty was assigned when f₁₀ is calculated from a reported mean or median value (40 %). To derive f₁₀ in these cases (from the mean, median, or bins), we used the average grain size distribution from the Coppo et al. (2022a) dataset to obtain linear least-square regression equations for f₁₀ vs. mean (R²=0.71), f₁₀ vs. median (R²=0.81), and f₁₀ vs. f₂₀ (R²=0.88), among others (see site-specific notes for more details), except when f₁₀ calculated from the mean grain size was retrieved directly from Albani et al. (2014). When no grain size measurements are available for a given site, only for nearby sites (<100 km away) that are comparable in terms of their geomorphological setting, then the same value for f₁₀ is used in both sites, with an extra 10 % uncertainty added for the site with no data. The same is true for sites that include grain size data for a different time window than the one considered.

Figure 3 Criteria for assigning uncertainty to the <10 µm grain size fraction. GSD: grain size distribution.

For marine sediment cores, it is rarely the case that grain size information is available (only 6 % of sites). For sites with no grain size information, if they are within 500 km from a site with grain size measurements, and provided the two sites have similar bathymetries and are not separated by significant bathymetric features, then the same f₁₀ value is used for both sites, with an extra 10 % uncertainty added for the site with no measurements. For all other marine sediment sites, f₁₀ values are assumed based on their downwind distance from known dust sources. If sites are located >2000 (1000–2000, <1000) km downwind from known dust sources south of 15^∘ S, we assume that f₁₀=1.00 (0.75, 0.50) ± 60 %. For locations in the ocean north of 15^∘ S where dust sources are more intense, these threshold distances are instead 3000 and 1500 km. The downwind directions from dust sources are qualitatively defined based on Holocene and LGM dust deposition maps from Mahowald et al. (2006), although dust sources not considered in this study, such as Alaska, are also considered.

2.8 The use of thorium-230 normalization for total sediment flux calculations in marine sediment cores

It has long been recognized that the traditional method of obtaining lithogenic mass accumulation rates between two dated horizons in marine sediment cores (e.g., Mortlock et al., 1991) cannot be directly attributed to pelagic sedimentation from the water column due to the process of sediment focusing, which is the lateral transport of sediment acting on the seafloor and along different depths in the ocean. Instead, the majority of lithogenic flux estimations from studies of marine sediment cores that are interested in the vertical fluxes of particles use the ²³⁰Th normalization method, which accounts for sediment focusing and provides estimations of pelagic mass accumulation rates. This method was first proposed by Bacon (1984) and is fully described in Francois et al. (2004). Briefly, the main assumption of the method is that the flux of scavenged ²³⁰Th to the seafloor is equal to the decay production rate of ²³⁰Th from ²³⁴U dissolved in the overlying water column. This is a reasonable assumption given the short residence time of ²³⁰Th in ocean water and the fact that its removal from the dissolved pool occurs mostly through adsorption to sediments derived vertically in the water column, a process known as proximal scavenging.

For the case of marine sediment cores, DMAR is calculated using the ²³⁰Th normalization technique as the product of the sediment bulk mass accumulation rate (SBMAR) and EC:

$$\mathrm{DMAR}=\mathrm{SBMAR}\times \mathrm{EC},$$

where

$$\begin{array}{}\text{(1)}& \mathrm{SBMAR}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{230}}\times z}{\mathrm{Th}\mathrm{\_}{\mathrm{230}}_{\mathrm{xs}}^{{}^{\circ }}}}\text{(e.g., Francois et al., 2004)}.\end{array}$$

Here, β₂₃₀ is the decay constant of ²³⁴U (and the production constant of ²³⁰Th) throughout the water column, with a value of 2.562 × 10⁻⁵ dpm cm⁻³ kyr⁻¹ (Costa et al., 2020); z is the water depth to a given site; and Th_230${}_{\mathrm{xs}}^{{}^{\circ }}$ is the decay-corrected excess ²³⁰Th activity (in dpm g⁻¹).

Another advantage of the ²³⁰Th normalization method is that it provides DMAR estimates that are, for the most part, independent of the chronology of sediment accumulation. This is particularly relevant for high-resolution DMAR studies, as large uncertainties exist for DMAR when using the traditional method if estimates are attempted for two horizons located close to each other in age.

The relative uncertainty in deriving SBMAR based on the ²³⁰Th normalization method is 30 % (Fig. 4), based on calibration studies (Henderson et al., 1999; Scholten et al., 2001; Yu et al., 2001a, b).

Figure 4 Criteria for assigning uncertainty to the sediment bulk mass accumulation rate when using the ²³⁰Th normalization technique for deriving dust fluxes from marine sediment cores.

2.8.1 Extra uncertainty in DMAR due to downslope sediment flow

If particles settling into the studied site directly from the overlying water column have the same ²³⁰Th activity as particles that have previously been advected laterally for some distance, then Eq. (1) can be used as discussed before to calculate the pelagic SBMAR. This is typically the case when laterally advected sediment is resuspended by bottom currents from a position in the seafloor at a similar water depth compared to the studied site (Francois et al., 2004). However, when the ²³⁰Th activities of the sediment components differ, inaccuracies may appear in the calculation of pelagic SBMAR. In principle, this may be the case when laterally advected sediment is originated at positions on the seafloor that are shallower than the studied site (note the dependence of SBMAR on the water depth in Eq. 1) and transported by either downslope bottom currents or by intermediate nepheloid currents (Francois et al., 2004). This is favored in regions of the seafloor with high regional bathymetric gradients, such as along the continental slope or in the open ocean close to bathymetric features (e.g., seamounts, aseismic ridges).

Only a few studies have carried out bathymetric analyses to derive explicit corrections to their SBMAR estimates to account for this potential effect of downslope sediment flow. One such example is from a site located at the foot of the Sierra Leone Rise in the equatorial Atlantic Ocean (EN066-29GGC, Francois et al., 1990). Based on a detailed analysis of the surrounding bathymetry, those authors defined a worst-case scenario in which all laterally transported sediment originated from the topmost part of the rise (with the greatest possible difference in water depth from the studied site and thus the greatest deviation in ²³⁰Th activity). In this scenario, the true pelagic SBMAR was overestimated by 36.5 % compared to the SBMAR calculated with Eq. (1) (Francois et al., 1990). Nonetheless, the true overestimation was probably lower, as a more realistic scenario is one in which laterally transported sediment was not fully derived from the top of the rise but partially from different steps along the rise at different water depths (Francois et al., 1990).

Instead, a study that looked at six cores <50 km from each other in the Juan de Fuca Ridge found no systematic differences in ²³⁰Th-normalized sediment fluxes, despite the hundred-meter-scale relief between sites (Costa and McManus, 2017). Whether this is an indication of the lack of sensitivity of ²³⁰Th-normalized sediment fluxes to downslope sediment flow at a global level remains to be determined. Because of this uncertainty, we did not attempt to apply corrections to ²³⁰Th-normalized DMAR estimates from marine sediment cores to account for this effect here. Another reason for not doing so is that it would require an in-depth, site-specific bathymetric analysis, which is out of the scope of this study. Instead, we raised the uncertainty in SBMAR for those sites that have probably experienced sediment focusing through downslope flow. We did so for sites that (i) are located at positions on the seafloor with high regional bathymetric gradients, (ii) are located at relatively deep positions compared to their surroundings, and (iii) have experienced sediment focusing. We evaluated the first two criteria by using GEBCO_2022, a global bathymetric grid at 15 arcsec horizontal resolution (GEBCO Compilation Group, 2022). For each marine sediment core site in Paleo±Dust, we took a 5^∘ latitude by 5^∘ longitude area centered at the site and calculated the difference between the 5th and 95th percentiles in bathymetry for all GEBCO_2022 grid cells located in that area. We also calculated this value for a reference marine sediment core site: EN066-29GGC (Francois et al., 1990). Any given marine sediment core site in Paleo±Dust where this 5th–95th percentile difference was equal to or higher than half of the 5th–95th difference for the reference site was considered to satisfy condition (i). Condition (ii) was satisfied for those sites located at a water depth of at least the median of that of the area defined for each site. Finally, condition (iii) was satisfied for those sites with a sediment focusing factor greater than 1, where the focusing factor is the ratio of laterally advected to vertical sediment flux, as defined by Francois et al. (2004) and retrieved for each site either from the original study or from the compilation by Costa et al. (2020).

Those sites that passed the above criteria were assigned an extra uncertainty component for SBMAR, which was calculated for each site based on a realistic scenario for reference site EN066-29GGC (half of the worst-case scenario, Francois et al., 1990), and proportionally to the site's focusing factor. The L2 norm was used to combine this extra uncertainty component to the base relative uncertainty component of 30 % common to all SBMAR estimates for marine sediment cores arising from ²³⁰Th normalization (Fig. 4).

3 Results

Paleo±Dust consists of a total of 285 pi-HOL and 209 LGM sites, of which approximately a third are sites published since 2016 and not included in previous paleo-dust deposition flux compilations (Table 1). Of all the sites, 52 % correspond to loess, 39 % to marine sediment cores, 4.7 % to peat cores, 2.6 % to lake cores, and 1.6 % to polar ice cores. All peat sites are new to this compilation, with 65 % of the sites published since 2016. Loess is the archive type with the highest number of new sites published since 2016 (89), and, except for peat sites, it is also the archive type with the highest percentage of new sites (35 %). All DMAR determinations from peat sites are for the pi-HOL, which can be explained by the fact that the vast majority of peat bogs globally only started to form during the last deglaciation (Yu et al., 2010).

Table 1 Number of sites in Paleo±Dust for the pre-industrial Holocene (pi-HOL, year 1850 CE–11.7 ka) and Last Glacial Maximum (LGM, 19.0–26.5 ka) with a dust deposition rate determination. For the mean LGM/pi-HOL dust flux ratio (<10 µm fraction) calculations, only sites with both pi-HOL and LGM determinations were considered.

While there are several studies that report dust measurements in non-polar ice caps and mountain glaciers for the pi-HOL and LGM, the reported quantity is the particle number concentration in all cases (Fisher, 1979; Thompson et al., 1989, 1995, 1997, 1998; Clifford et al., 2019; Beaudon et al., 2022) except for a site on the Penny ice cap for which the mass concentration is reported (Zdanowicz et al., 2000). Many of these studies also lack grain size data. The derivation of mass deposition rate from particle number concentration alone is not straightforward and requires several assumptions (Kohfeld and Harrison, 2001), including the grain size distribution and the distribution of particle density with grain size. We thus follow Kohfeld and Harrison (2001) and do not derive dust deposition rates in these cases. Thus, no ice cores from non-polar sites are included in Paleo±Dust.

There is a clear Southern versus Northern Hemisphere asymmetry in the number of sites in Paleo±Dust (Fig. 5a–b). This is mostly evident in the number of continental dust archives. There are only a few loess studies from the Pampean region in South America, the main loess belt in the Southern Hemisphere, from which DMAR estimates can be obtained (Kemp et al., 2004; Torre et al., 2019; Coppo et al., 2022b). South American loess is restricted to latitudes <40^∘ S, and thus it only archives dust from low- to mid-latitude sources in South America as well as from northernmost Patagonia. Emissions from southernmost Patagonia's main dust sources during the pi-HOL are captured by peat bogs on the Malvinas Islands (Monteath et al., 2022) and possibly also peat bogs located to the east of the Andes (Vanneste et al., 2015, 2016). These latter studies report on sites that are located upwind of the main present-day dust sources but have captured pi-HOL dust from smaller sources in the region. Instead, the bulk of dust emissions from central and southern Patagonia (40–50^∘ S) since the LGM are largely unconstrained due to a lack of loess deposits downwind. It is thus critical that marine sediment cores are drilled in the open ocean off the eastern Patagonian coast with a focus on constraining dust fluxes. Elsewhere in the Southern Hemisphere, there is a complete dearth of DMAR estimates from loess or loess-like deposits in Australia, and only one pi-HOL DMAR estimate in southern Africa (Brunotte et al., 2009). In the Northern Hemisphere, the high latitudes in Asia (i.e., Siberia) constitute the area with the scarcest coverage of DMAR constraints. In turn, Europe is the region with the highest number of new DMAR estimates since 2016.

Figure 5 Distributions of (a) the sites for the pre-industrial Holocene (pi-HOL) time period and (b) the sites for the Last Glacial Maximum (LGM) time period, the (c) pi-HOL and (d) LGM <10 µm dust deposition rates, and the (e) pi-HOL and (f) LGM <10 µm dust deposition rate uncertainties. The coastlines in all panels correspond to those of the present day.

Based on 126 sites with paired LGM and pi-HOL DMAR determinations (Fig. 5c–d), the global LGM:pi-HOL ratios of DMAR and DMAR₁₀ are 3.3 ± 0.7 and 3.1 ± 0.7, respectively (1σ, Table 1). The latter value is equivalent within error to that determined exclusively based on loess DMAR₁₀ determinations (3.2 ± 0.7, N=47). Based exclusively on intermediate-range dust deposition as archived in marine sediments, this ratio is significantly lower, albeit within the same order of magnitude (2.2 ± 0.5, N=72). Finally, based only on dust archived in polar ice cores, this ratio goes up dramatically to 16.8 ± 4.2, although this is based on only four sites. Three processes are responsible for higher LGM/pi-HOL dust deposition rate ratios for polar ice cores compared to loess: (1) the effect of dust transport in the atmosphere and dry deposition, by which small changes in dust emission intensity translate to small changes in close-to-source dust deposition and to bigger changes in remote dust deposition (Lambert et al., 2008); (2) the effect of mid-latitude precipitation on wet scavenging of dust en route to the poles (Markle et al., 2018), whose variability at glacial–interglacial timescales may imply an amplification of dust deposition rate variability at the poles compared to the low latitudes; and (3), specifically for the Southern Hemisphere, the activation or intensification of dust emissions from Patagonian sources at higher latitudes (compared to the present day) during the LGM implied more efficient transport of dust to Antarctica (Andersen et al., 1998; Petit et al., 1999; Lambert et al., 2008; Albani et al., 2012).

Relative 1σ uncertainties in DMAR and DMAR₁₀ mostly range between 23 %–99 % and 38 %–105 %, respectively (5 %–95 % percentiles, Fig. 5e–f). The least uncertain archive type for DMAR and DMAR₁₀ is polar ice (median: 23 % and 25 %, respectively, N=8), while the most uncertain types are marine sediments (44 % and 74 %, N=191) and loess (45 % and 68 %, N=258). However, the differences in median relative uncertainty among archive types are significantly lower than the differences among individual sites across archive types. There are no significant differences in uncertainty for LGM vs. pi-HOL DMAR and DMAR₁₀. The DMAR₁₀ determinations with the greatest uncertainties (the top 5 % percentile, with values >105 %) all correspond to loess deposits that lack several determinations, and so many assumptions are made (e.g., <10 µm grain size fraction, density, organic content, carbonate content), and/or have t_top and t_bottom values that are distinct within error but very close to each other.

4 Structure of Paleo±Dust

All Paleo±Dust files are included in the Supplement (supplementary_data.zip). This consists of two main tab-delimited text files containing the most important variables for each site (i.e., main_piHOL_tab.txt, main_LGM_tab.txt), two extra supporting tab-delimited text files with an expanded set of intermediate variables used to calculate the main set of variables (i.e., supporting_piHOL_tab.txt, supporting_LGM_tab.txt), and two separate text files (i.e., site_specific_notes_references.txt, site_ references.txt) containing site-specific observations and a list of references from where the necessary data were extracted.

The two main text files each contain 10 variables: locality (“siteName_region”), type (i.e., “ice-core”, “marine”, “loess”, “lake”, or “peat”), lat_N (i.e., latitude in degrees north between −90 and 90 to two decimal places), lon_E (i.e., longitude in degrees east between −180 and 180 to two decimal places), top-age_kaBP and bottom-age_kaBP (i.e., t_top and t_bottom, respectively, in thousands of years before present to two decimal places), DMAR_g/m2/a and DMAR-1sigma_ g/m2/a (i.e., mean DMAR and 1σ σDMAR for the defined age bracket, respectively, in g m⁻² a⁻¹ to three significant figures), and DMAR10_g/m2/a and DMAR10-1sigma_g/m2/a (i.e., mean DMAR₁₀ and 1σ σDMAR₁₀ for the defined age bracket, respectively, in g m⁻² a⁻¹ to three significant figures).

The two supporting text files each contain 24 variables, corresponding to the same 10 variables as the main files plus 14 extra variables: first_appeared_in_dataset (i.e., reference to the compilation where the site first appeared), top-age-1sigma_kaBP and bottom-age-1sigma_ kaBP (i.e., 1σ σt_top and σt_bottom, respectively, in thousands of years before present to two decimal places), thickness_m (i.e., h_thick defined by the age bracket, in m to two decimal places), depth-1sigma_ cm (i.e., 1σ uncertainty of the depth below the surface of the dated layers in cm to two decimal places), DBD_g/m3 and DBD-1sigma_g/m3 (i.e., DBD and 1σ σDBD, respectively, in g m⁻³ to no decimal places), SBMAR_g/m2/a and SBMAR-1sigma_g/m2/a (i.e., SBMAR and its 1σ uncertainty, respectively, only for type “marine”, in g m⁻² a⁻¹ to three decimal places), EC_adim and EC-1sigma_adim (i.e., EC (as a number between 0 and 1) and 1σ σEC, respectively, both to three significant figures), f10_adim and f10-1sigma_adim (i.e., f₁₀ (as a number between 0 and 1) and 1σ σf₁₀, respectively, both to three decimal places), and flag-marine-sed-downslope (i.e., a flag exclusive to type “marine”, where 1 means that the site is prone to contamination by downslope sediment movement and 0 means the site is not prone to such contamination).

5 Data availability

The Paleo±Dust data files are available in the Supplement (supplementary_data.zip). They are also available in the PANGAEA open access repository (https://doi.org/10.1594/PANGAEA.962969, Cosentino et al., 2024).

6 Conclusions

Paleo±Dust is an updated global paleo-dust deposition rate compilation for mean pi-HOL and LGM climate states that includes quantitative estimates of uncertainties. Paleo-dust flux measurements from peat bog cores are included in a paleo-dust compilation for the first time. By excluding deglaciation dust fluxes, Paleo±Dust better isolates mean interglacial and glacial dust fluxes than previous datasets. Site-specific age brackets allow the sub-sampling of mean dust deposition constraints for specific time windows of interest. Grain size information was extracted from the original studies to derive fluxes of dust with <10 µm diameter particles.

The main feature of Paleo±Dust is the inclusion of site-specific dust flux uncertainties that are consistent across paleo-dust archive types, i.e., that can be used to gauge the relative accuracy of dust flux constraints between sites of different geologic nature. These new uncertainty data may also be of use as a criterion for selecting a subset of samples for comparison against dust simulation output, to tune dust emission in Earth system models using dust deposition flux data and a weighted approach based on the uncertainty of proxy data, or to construct distribution-based global interpolation maps of paleo-dust deposition rates through Bayesian approaches or Monte Carlo or bootstrapping experiments for use, for example, as input for biogeochemical models.