Atmospheric dust flux in northeastern Gondwana during the peak of the late Paleozoic ice age

The silicate mineral fraction of shallow marine carbonates archives dust contributions to the Central Persian Terranes along the northeastern margin of Gondwana ( ∼ 30 ° S paleolatitude), enabling reconstruction of atmospheric dust loading and circulation for intervals of the late Paleozoic ice age. The Central Persian Terranes hosted cyclic deposition of warm water carbon-ates from middle Pennsylvanian to earliest Permian time, and our data set includes two ∼ 28 m sections from the Moscovian and As-selian sampled at 20 cm intervals. Bounding surfaces between successive cycles (high-frequency sequences) are recognized by either abrupt basinward shifts in facies or subtle exposure features; these high-frequency sequences range from 1 m to 5 m thick and are interpreted to record glacioeustatic variations. Time series analysis of the dust fraction through the studied interval supports the hypothesis of orbital forcing for the dust signal. The stratigraphic pattern of the dust flux indicates minimal flux during interglacial highstands (0.19–0.27 g/cm 2 / kyr) and peak flux during glacial lowstands (3.77–4.57 g/cm 2 /kyr) after accounting for hiatal time at sequence boundaries. Grain size analysis of the dust for all samples (n = 230) reveals modal sizes (volume-based) of 1–15 µ m through the Moscovian interval and 10–75 µ m through the Asselian interval. Dust deposition increased during glacial times relative to interglacial times by a factor of 16 to 19. Additionally, the Asselian interval exhibits higher dust flux overall relative to the Moscovian interval, which is interpreted to reflect the more extreme ice-house conditions of the Asselian. Variation in the dust content through the studied sections provides an indicator of temporal changes in atmospheric loading that varied at both glacial–interglacial and higher-frequency ( < 10 4 yr) scales. Geochemical data reveal that the Arabian–Nubian Shield and southwestern Pangaea (South America) are the most likely sources of dust deposition in the Central Persian Terranes, with sources shifting during different phases. Increased dust flux during glacials likely reflects multiple factors, including enhanced aridity in the source region, exposure of shelf regions, and potential changes in winds. However, the discrepancy in model reconstructions of the amplitude of glacial–interglacial dust variations indicates that increased production of dust sourced by dynamic glaciation played a large role in enhancing dust flux during glacial phases.


INTRODUCTION
The late Paleozoic Earth is our most recent analog for the icehouse conditions that have existed during the late Cenozoic; both are distinguished by repeated, high-magnitude glacialinterglacial climate shifts. The Cenozoic record contains Earth's best-known loess deposits, and the study of these as well as dust recovered from ice-and ocean cores demonstrates higher dust fluxes during glacials relative to interglacials (2-20 times for the Last Glacial Maximum; e.g., Kukla, 1987;DeAngelis et al., 1987;An et al., 1991;Rea, 1994;Petit et al., 1999;Kohfeld et al., 2005;Marino et al., 2005;Anderson et al., 2006). A growing data set indicates that the Permo-Carboniferous was marked by a particularly dusty atmosphere that was recorded in the form of abundant and widely distributed loess deposits (e.g., Murphy, 1987;Johnson, 1989;Soreghan et al., 2008;Giles et al., 2013;Sweet et al., 2013;Foster et al., 2014;Pfeifer et al., 2020) and in shallow marine units with large eolian contributions (e.g., Carroll et al., 1998;Sur et al., 2010a,b), albeit most of the data supporting this come from low paleolatitudes (western equatorial Pangaea).
Pennsylvanian-Permian strata are cyclic in many regions; within marine and paralic systems of the low-mid paleolatitudes, these cycles are commonly interpreted as glacioeustatic in origin (e.g., Heckel, 1986;Algeo and Heckel, 2008), reflecting repeated oscillations between glacial and interglacial climates. In shallow marine carbonate systems that formed isolated from fluvial influx, carbonates commonly include a fine detrital silicate fraction (e.g., Cecil et al., 2003;Soreghan et al., 2008;Sur et al., 2010b;Carvajal et al., 2018;Oordt et al., 2020) that records significant fluctuations in atmospheric dust input. For example, data from Sur et al. (2010b) on dust from Pennsylvanian carbonates revealed large variations in glacial (lowstand)-interglacial (highstand) dust flux ratios with absolute values that greatly exceed those of the Plio-Pleistocene by up to an order of magnitude. Given that nutrients carried by dust can potentially fertilize primary production (e.g., Boyd et al., 2004;Krishnamurthy et al., 2010;Kanakidou et al., 2012), such large shifts in dust flux might have had similarly large repercussions for carbon cycling during this interval of Earth history Sur et al., 2015;Sardar Abadi et al., 2020).
In this work, we look beyond these lowpaleolatitude records to assess atmospheric dust flux and temporal shifts in this flux in the subtropics of northeastern Gondwana (Fig. 1). To assess whether atmospheric dustiness extended beyond the equatorial region and determine any links between atmospheric dust flux and glacial-interglacial cycles, we examined the Permo-Pennsylvanian strata from the eastern subtropics of Pangaea in present-day Iran. The Central Persian Terranes near the Arabian margin of northeastern Gondwana were situated along the southern margin of the Paleo-Tethyan Ocean at ∼30°S and record deposition of warm water, shallow marine carbonates throughout middle Pennsylvanian-earliest Permian time (Sardar Abadi et al., 2019). The broad distribution of carbonates of this age across the Central Persian Terranes indicates the presence of a widespread carbonate platform isolated from siliciclastic influx except for that delivered by eolian influx. To establish the record of dust flux in this region of northeastern Gondwana, we combine sedimentologic measurements with geochemical proxies to reconstruct atmospheric dust concentrations through the study intervals and assess the magnitude, source(s), and possible forcing(s) of the dust variations. Results augment a growing data set on global dust records for the late Paleozoic ice age and illuminate glacial-interglacial variations in atmospheric dust loading and circulation in a region relatively proximal to Gondwanan ice centers.

GEOLOGICAL SETTING
The Central Persian Terranes encompass continental segments including the Alborz Basin, Sanandaj-Sirjan Zone, and Central Iran positioned along the Arabian margin of Gondwana at ∼30° S during Pennsylvanian-Early Permian time. Sedimentologic and paleontologic proxies recovered from the Central Persian Terranes record the formation of warm, shallow marine carbonates in proximity to Gondwanan glaciations throughout the Late Pennsylvanian-earliest Permian (Sardar Abadi et al., 2019). Sea surface temperatures (SST) modeled for the Paleo-Tethys and Panthalassic Ocean indicate significant contrasts between the two regions, with relatively warm sea surface temperatures (SSTs) affecting the Central Persian Terranes and cooler SSTs along the northwestern margin of Gondwana. Despite the warm SSTs of the Central Persian Terranes (present-day Iran) recorded in both data and modeling results for this interval (Sardar Abadi et al., 2019), coldclimate conditions prevailed immediately south of the study area (∼1000 km), evinced by coeval glacial deposits preserved across the Arabian region throughout the upper Pennsylvanian and Lower Permian record (Melvin and Sprague, 2006;Bussert, 2014).

DEPOSITIONAL SETTING
In the Central Persian Terranes, Pennsylvanian-Permian strata accumulated in mid-to inner-ramp facies belts across an extensive shallow marine carbonate ramp (Sardar Abadi et al., 2019). The Moscovian Absheni Formation (middle Pennsylvanian) in the S-S Zone preserves a variety of phylloid algal deposits and microbial carbonates that formed at or near wave base and broad oolitic shoals that accumulated above wave base. The Asselian Emarat Formation (Lower Permian) in the Alborz Basin preserves a mid-ramp facies association including a brachiopodal-bryozoan packstone facies and hummocky, crosslaminated, siltstone facies occurring in thin layers (<10 cm), which is indicative of deposition above storm wave base. The Emarat Formation includes a shallower-water facies association represented by thick, bioclastic to oncolitic shoals with micritic mudstone, which is indicative of variable high-to lowenergy conditions across an inner-ramp set-ting. Table 1 provides a descriptive summary and interpretation of each facies.
Testing for temporal variations in dust flux requires establishing a sequence-stratigraphic framework, which in turns requires detailed facies analysis. For samples from the Central Persian Terranes (Iran), we relied on the detailed facies analyses by Sardar Abadi et al. (2019), which are summarized in Table 1.

Extracting the Dust Signal
To extract the dust component, we used the protocol developed by Sur et al. (2010a). In brief, all samples were cleaned of exterior  debris, crushed to gravel size, rinsed with distilled water, and subjected to dissolution with 2N hydrochloric acid (HCl) to eliminate carbonate. The insoluble residue was then combusted at 500 °C for 12 h to remove organic matter and oxidize pyrite iron. Subsequently, samples were treated with the citrate-bicarbonate dithionite (CBD) method (Sur et al., 2010a) to remove iron oxides. The resultant silicate mineral fraction (SMF) was freeze-dried and then examined for the presence of chert or other authigenic material, which occurred rarely. Only 1.1% (three of 266 samples) contained authigenic components, and in each case they occurred in association with micritic mudstone (see Table 1). Once authigenic components were removed (by hand), the SMF is considered to represent the detrital (dust) fraction; hereafter we refer to this fraction as "dust."

Measuring Dust Composition and Grain Size
The dust was examined under reflected light to assess the presence of any clearly authi-genic components; smear slides and microprobe analyses were used to assess mineralogy. Grain size distributions were measured with a Malvern Mastersizer 3000 (laser particle size analyzer [LPSA]) using the Hydro SM dispersion unit. Prior to LPSA, each sample was suspended in 20 ml of distilled water with a few drops of sodium hexametaphosphate (dispersant) and sonicated for one minute. Level plots (Cleveland, 1993) to visualize number-and volume percent distributions of grains were produced using the filled.contour function in the Program R ver. 3.4.3 (R Core Team, 2019).

Bulk Geochemistry
Samples (n = 266) were analyzed for bulk geochemistry using a handheld Bruker Trace IV-SD X-ray fluorescence spectrometer (XRF) scanner on cut pieces of hand samples that were at least 2 cm thick. Major elements were acquired using an accelerating voltage of 15 kv, current of 30 μA, vacuum of <10 Torr, and a standing time of 90 seconds with no filter. Trace elements were measured using an accelerating voltage of 40 kv, current of 17.1 μA, vacuum of <10 Torr, and a standing time of 60 s with a titanium aluminum filter.
Variations in the elemental composition were summarized using principal component analysis (PCA). Broken stick models (Legendre and Legendre, 2012) were used to assess the number of informative principal component axes.
Redundancy analysis (RDA) is a direct ordination technique allowing for tests of multiple drivers of a matrix of compositional data (Legendre and Anderson, 1999). We performed an RDA to examine how the chemical composition of samples changed with percent dust, proximity to cycle boundary, dust volume distribution (percent of dust volume <10 μm), and interval (Moscovian or Asselian). To limit collinearity in our models we tested for Variance Inflation Factors (VIF) and only included predictor variables with VIF < 4 (Gross, 2003). Driver variables demonstrated no variance inflation (all VIF < 4). All ordination analyses were conducted using the R package vegan (Oksanen et al., 2015) in program R ver. 3.5.1 (R Core Team, 2019).

Sr Isotope Methods
The dust fractions derived from eight samples representing different facies from both intervals were selected for provenance analysis using Sr isotope geochemistry. Approximately 50 mg of material for each sample was digested for radiogenic isotopic analysis using a combination of double distilled HF-HNO 3 , HClO 4 , and HCl acids in screw-top teflon beakers under clean lab conditions (Gleason et al., 2002). Two overnight hotplate digestions at 120 °C in concentrated HF-HNO 3 (5:1 ratio) were employed to break down silicates, followed by perchloric acid digestion at 160 °C to destroy fluorides and any remaining organics. The residues were taken up in 6N HCl and converted to nitrides prior to matrix separation. Sr was separated from matrix using miniaturized chromatographic columns (Eichrom Sr-Spec resin) following procedures of Gleason et al. (2002) and Stancin et al. (2006). Sr isotope ratios were determined at the University of Michigan on a Triton Plus (Thermo Scientific) thermal ionization mass spectrometry (TIMS) instrument equipped with eight Faraday collectors in static mode configuration. The mean value for the Sr standard NBS 987 was 87 Sr/ 86 Sr = 0.710245 ± 0.000014 (n = 11), requiring no corrections to data by comparison with the accepted 87 Sr/ 86 Sr ratio for NBS987 (0.710248; McArthur et al., 2001). Total process Sr blanks were ∼50 pg (<1/1000 of the sample processed) requiring no corrections to the data.

Time Series Analysis
To determine sedimentation rate and thus temporal length of the study interval, we used time series analysis to assess the presence of Milankovitch forcing in the dust record. Astrochronologic testing was conducted using the astrochron package (Meyers, 2014). Prior to time series analysis, dust values were log-transformed, detrended, and interpolated to a uniform spacing of 0.25 cm. Evolutive harmonic analysis (EHA, Meyers, 2014) was conducted to determine significant frequencies, defined as those exceeding the 90% significance level (Supplementary Figs. S1C and S2C 1 ; see vertical dashed lines). Next, astronomical target frequencies were set using the La2004 astronomical solution (Laskar et al., 2004), extrapolated to Asselian and Moscovian times. For the Moscovian, target cycle periods (inverse frequencies) are "E" (long orbital eccentricity) = 405 k.y., "e" (short orbital eccentricity) = 128 and 97 k.y., "o" (obliquity) = 32 k.y., and "p" (precession) = 21.6 and 17.6 k.y. The significant frequencies from the dust records obtained using EHA were tested with the set of astronomical target frequencies using the Average Spectral Misfit (ASM) method (Meyers and Sageman, 2007). The ASM method was applied to the SMF records (Supplementary Figs. S1F and S2F), testing sedimentation rates ranging from 1 cm/kyr to 10 cm/kyr.
To calculate the sample-specific sedimentation rate, we adjusted sedimentation rate to match an element expected to have a constant deposition rate. We selected zinc (Zn) as Zn exhibited low variance across the sampling intervals, Zn does not correlate (positively or negatively) with dust, and stratigraphic variations in Zn should result mainly from variable influx of other depositional components, notably, carbonate (Rothenstein et al., 2012).
Sample-specific sedimentation rate was calculated as where SR i is sedimentation rate for given sample (i), Zn A is average %Zn across the interval, ASR is average sedimentation rate calculated with ASM, and Zn i is %Zn for given sample (i). Time between samples was calculated as where TBS is time between samples, D is distance between samples (cm), and SR is sedimentation rate, all for a given sample i. ∑TBS i = total temporal length of the interval. Next, the total mass accumulation rate was calculated as where TMAR is total mass accumulation rate, SR is sedimentation rate, and fDBD is faciesspecific dried bulk density. Finally, the mass accumulation rate of SMF was calculated as where SMAR is the mass accumulation rate of SMF, TMAR is total mass accumulation rate, and %SMF is the weight percent of silicate mineral fractions. The calculation and data are provided in Supplementary Information 1 (see footnote 1).

Dust Cycle Modeling
To obtain reference estimates of dust deposition during characteristic glacial and interglacial intervals of the late Paleozoic, we adopted the reference modeling of late Paleozoic dust cycling during interglacial and glacial intervals described by Marshall et al. (2016). Marshall et al. (2016) used simulations with the Community Climate System Model version 3 (CCSM3) of dust emission, transport, and sedimentation with radiatively inactive dust and tuned them to one observation: a reconstructed dust deposition rate of 32 g/m 2 /yr -(3.2 g/cm 2 /kyr) in Middle Permian western Brazil. Marshall et al. (2016) was used to generate a reference estimate for glacial-interglacial dust deposition, because it is to date the only study to tune dust cycle simulations to reconstruct dust deposition for both glacial and interglacial intervals of the late Paleozoic. The dust sizes tracked in the dust cycle simulations are in four size bins covering the interval from 0.1 μm to 10 μm (0.1-1 μm, 1-2.5 μm, 2.5-5 μm, and 5-10 μm). In present day climates, dust this size has a mean atmospheric lifetime ranging from 1.4 days for 5-10 μm size dust to 7.6 days for 0.1-1 μm size dust, and thus it is transported over regional to global scales (Mahowald et al., 2006). Further details about the reference modeling are included in Table 2 and the Supplementary Materials (see footnote 1) for this article.

Depositional Cycles
The carbonate strata of the Moscovian Absheni Formation and the Asselian Emarat Formation occur in repetitive successions (cycles, Fig. 3) relatively analogous to the cyclothems recorded in Pennsylvanian-Permian sections in paleo-equatorial regions (e.g., Heckel, 1986;Algeo and Heckel, 2008). North American cyclothems record facies migrations that predominantly reflect upwardly shallowing patterns (references above). Although carbonate cycles in the Central Persian Terranes consist exclusively of shallow marine facies in an interval exhibiting an overall upwardly shallowing trend, individual high-frequency sequences record facies associations exhibiting an upwardly thickening stacking pattern (ranging from 2 m to 5 m in thickness) that is inferred to reflect upward shallowing. Additional evidence for upward shallowing within each sequence is commonly recorded by an upward increase in abundance of skeletal clasts, recording shallowing into higherenergy (subtidal) environments. This resembles the cyclostratigraphic pattern preserved in approximately coeval Bolivian cool-water carbonates (Carvajal et al., 2018) that developed at equivalent paleolatitudes, but along the northwestern margin of Gondwana.
The high-frequency sequences are separated by bounding surfaces that record either abrupt deepening of facies (i.e., flooding surfaces) or abrupt shallowing of facies (those reflecting subaerial exposure). No peritidal facies occur. Where (mostly micro) evidence for subaerial exposure occurs, it is "abnormal" in that the succession records incomplete facies progradation entirely within subtidal facies, consistent with the influence of high-frequency glacioeustatic oscillations (e.g., Wright, 1992;Soreghan, 1997;Burgess, 2016).

Prominent Surfaces
The pedogenesis (detailed below), dissolution fissures, micro-scale discontinuity, and evidence for over-compaction, but that also preserve evidence for flooding immediately superjacent to surfaces with exposure features. Key phenomena are discussed below.

Incipient Pedogenesis
Incipient pedogenesis is evident in centimeter-scale horizons that exhibit a combination of alveolar microstructures and circumgranular cracking that is commonly associated with ferruginous siltstone. Alveolar microstructures and circumgranular cracks are typically developed within thin (<2 cm) iron-rich micritic laminae commonly overlain by ferruginous siltstone. Alveolar texture is recognized as clusters of very small (<200 μm diameter) voids separated by thin micritic septae (Fig. 4A). These alveolar structures are typically surrounded by brownish micrite and partially filled with calcite spar. Circumgranular cracks (Fig. 4B) occur as thin, commonly curvilinear, calcite spar-filled joints surrounding regions (<2 cm) of micrite that give rise to a glaebular fabric (Goldstein, 1988). Pedogenic surfaces encountered throughout the Absheni Formation and Emarat Formation share some attributes analogous to paleosols that cap fourth-order cycles of low-latitude Pennsylvanian carbonates (e.g., Goldhammer et al., 1991), but are much more weakly developed relative to these examples with recognition mostly in thin section.

Dissolution Fissures
Subtle dissolution fissures occur commonly at the tops of the meter-scale cycles and consist of oolitic grainstone throughout the Moscovian Absheni Formation (Fig. 4C). These dissolution features comprise narrow (millimeter-to centimeter-scale) subvertical cracks penetrating to a depth of several centimeters that truncate grains and are secondarily filled with a mixture of material originating from overlying sediments such as silt-sized quartz and lithoclasts. Relevant

indicated by white arrows). (C) Oolitic grainstone with subvertical cracks that is secondarily filled with a mixture of silt-sized quartz and lithoclasts (indicated by white arrows). (D) Overcompacted skeletal fragments with iron-rich cement filling intergranular interstices.
examples of small-scale solution features produced at cycle tops are associated with highfrequency sea-level fluctuations for the Pennsylvanian-Permian intervals (Foos, 1996).

Micro-Scale Discontinuity
Micro-discontinuity takes the form of centimeter-scale abraded surfaces that cross-cut lithified sediments and are covered with ferruginous siltstone, fine-grained sandstone, or intraclasts associated with host carbonate (Fig. 5C in Sardar Abadi et al., 2019). Such surfaces are interpreted to reflect deflation that occurred upon subaerial exposure during sea-level fall.

Flooding Surfaces
Flooding surfaces are sharp horizons across which abrupt deepening of facies occurs; they are commonly associated with subjacent evidence for incipient pedogenesis. Such surfaces have thin layers (2-3 cm) exhibiting evidence for overcompaction such as grain-to-grain suturing (of skeletal or oolitic grains) and grain cracking with iron-rich cement filling intergranular interstices (Fig. 4D). Lack of early cement rims (particularly fine cloudy drusy cements) suggests a low carbonate production environment that resulted in a low sedimentation rate (Mamet and Boulvain, 1990). Dysoxic bottom waters may have supported the proliferation of iron bacteria in the sediment, thereby causing the distinctive red sediment color. Analogous features such as abrupt contacts in the absence of subaerial exposure associated with high-frequency sea-level fluctuations were described from the upper Pennsylvanian of Arrow Canyon (Nevada; Bishop et al., 2010) and the Swope Limestone (Kansas; Heckel, 1983).

Stratigraphic Distribution and Dust Composition
Amount (weight %) of the dust fraction through the studied sections varies systematically, with dust content increasing in proximity to inferred high-frequency sequence boundaries. Excluding sequence boundaries, intracyclic variation of dust amounts through carbonate facies of the Moscovian interval range from 0.14% to 9.73% and from 0.04% to 10.14% through carbonate facies of the Asselian interval. In contrast, dust contents vary from 10.84% to 80.80% (constituting siltstone or very fine sandstone in some cases) in proximity to sequence boundaries- either coinciding with boundaries or peaking within a few centimeters above or below (further details in Supplementary Fig. S3; see footnote 1). Grain-size distributions (in volume %) of the dust through the Moscovian interval positively track variations in dust content (weight percent), wherein coarser modes correspond to larger amounts. The majority of the dust exhibits modes between ∼1-15 μm through the Moscovian interval. However, near sequence boundaries, particle modes (as measured by volume %) increase to 20-80 μm, with the coarsest material occurring in the most siliciclastic-rich intervals near the top of the Moscovian (Fig. 3). Similarly, in the Asselian, size modes and amounts of the dust fraction positively track one another; however, overall grain sizes are significantly coarser, with the majority of the dust between ∼10-75 μm and exceeding 100 μm near sequence boundaries. Samples are commonly fine-skewed with a significant population of grains in the 1-10 μm range (Fig. 3). The occurrence of this fine fraction (1 μm to 10 μm) indicates considerable numbers of these particles in the atmosphere, as volume percent distributions are biased by low numbers of larger particles.
Petrographic study (binocular, smear-slide analysis, microprobe) of the dust extracted from the carbonate facies of both study sections reveals quartz and aluminosilicates (clays) as predominant components, followed by muscovite, zircon, and rutile, in decreasing abundance. Most of the sand-sized (62-250 μm) quartz grains recovered from the carbonate facies occur as sub-rounded to rounded grains, and some have iron-oxide coatings (Fig. 5). However, finer (silt and smaller) particles tend to be angular.

Geochemistry (Detrital and Paleo Productivity)
Bulk-rock geochemistry (from handheld XRF) was compared using ordination techniques (PCA and RDA) to examine the origin of the dust and provenance variation in dust sources.
To assess correlations between selected chemical elements and dust and to measure statisti-cal relevance of correlations, we used principal component analysis (PCA). The PCA of the XRF elements includes two significant axes (Fig. 6). The first axis explains 40.3% of variation in the measured elements and was highly correlated with dust weight % (Pearson's r = 0.88). Additionally, the first PCA axis correlates positively with many elements commonly used as proxies for detrital material, including Ti (r = 0.71), Zr (r = 0.8), Si (r = 0.92), K (r = 0.78), Al (r = 0.92), and Rb (r = 0.64). The high statistical correlation of these typically terrigenous elements with both dust weight % and the first PCA axis is consistent with a detrital origin for the dust component. Calcium (Ca) is strongly negatively correlated with the first PCA axis (r = -0.93), reflecting the dilution effect of carbonate on the detrital contribution.
The second PCA axis explains 12.7% of the variation in the measured elements and is strongly positively correlated with Sr (r = 0.55), Zn (r = 0.76), Sn (r = 0.50), Sb (r = 50), and negatively with Mg (r = 0.88) and V (r = 0.76) (Fig. 6, Supplementary Table S1; see footnote 1). The elements correlated with the second axis are usually related in some manner to carbonate production and are negatively correlated with dust weight % (Fig. 6).
Zr and Ti are considered insoluble elements, and thus their ratio can be a sensitive index of shifts in provenance. The Zr to Ti ratio shifts between intrasequence and sequence boundary positions, with higher Zr:Ti at the cycle boundaries, suggesting provenance shifts from glacial to interglacial intervals (Fig. 8). Examination of grain size together with Zr:Ti reveals no relationship, strengthening the interpretation of a provenance control on this variation (Fig. 8). Absolute concentrations of Zr and Ti values also increase as much as two orders of magnitude at or near sequence boundaries (prominent surfaces).
To depict possible source regions for extracted dust from the Central Persian Terranes, we conducted a literature search for the Precambrian and Paleozoic sample chemistry from potential proximal source regions (i.e., regions either uplifted or sites of continental clastic deposition at the time). This review revealed 18 discrete source regions including regions of Arabia, east-central Gondwana, northern Gondwana, western Gondwana, India, and Australia (data compilation and related references are listed in Supplementary Information 2; see footnote 1). We then conducted a PCA of the chemical elements that are highly correlated with detrital proxies (r>0.5) using both our samples and those from the 18 potential sources from the literature to assess chemical similarities and thus likely provenance regions (Supplementary Fig. S4; see footnote 1). The resulting PCA has one significant axis and of the 18 potential sources, seven were selected as more likely sources as they have PC1 axis scores with values overlapping those of our samples (Supplementary Fig. S5; see footnote 1).
To further examine the relationships between our samples and the most probable source regions, we conducted an additional PCA using only those seven regions deemed to be the most likely dust sources (Fig. 9). This PCA resulted in four significant axes explaining 30.2%, 19.1%, 15.8%, and 10.1% respectively.

Sr Isotope Geochemistry
87 Sr/ 86 Sr i initial ratios (corrected for Rb/ Sr ratio and age of deposition) show the classic grain-size fractionation trend of Sr isotopic compositions commonly observed in eolian dust deposits, where finer fractions have higher 87 Sr/ 86 Sr than coarse fractions (Aarons et al., 2013). This is explained as a sediment-unmixing phenomenon, wherein more radiogenic (high 87 Sr/ 86 Sr) clay minerals are concentrated in the distal parts of eolian systems versus the less radiogenic or "feldspar dominated" (low 87 Sr/ 86 Sr) coarse fraction (Grousset and Biscaye, 2005). 87 Sr/ 86 Sr i ratios in this study range from 0.709 to 0.718 across middle Pennsylvanian and Early Permian glacial-interglacial stages and also overlap values for modern eolian dust deposits (Aarons et al., 2013). The 87 Sr/ 86 Sr values, grain size, and dust flux values are highly correlated (Table 3, Fig. 10, Supplementary Figs. S6 and S7; see footnote 1), suggesting the presence of a provenance signal in the Sr isotopic compositions (see discussion below).

Origin of the Dust
Paleogeographic evidence from the Central Persian Terranes documents widespread marine  carbonate deposition, far-removed from fluviodeltaic delivery of siliciclastic material, and thus the extracted silicate minerals represent atmospheric dust. The size distribution of these dust particles varies from a minority (<4%) of coarse grains (silt to sand, >10 μm) to a major-ity of fine grains (<10 μm). The coarse grains likely represent short-range transport particles of eolian sand that accumulated in near-coastal regions and were sourced ultimately from fluvial systems emanating from the northeastern Gondwanan margin, including proglacial systems of glaciers occupying the Arabian margin of northeastern Gondwana (e.g., Martin et al., 2008;Stephenson, 2008;Sardar Abadi et al., 2019). Finer grains represent long-range transport dust particles from sources ranging up to thousands of kilometers distant (Prospero, 1999). Sources for this material could include exposed plutonic and magmatic rocks associated with crustal accretion in the Pan-African belts from eastern Africa to southern Arabia or sources in southwestern Gondwana (South America; consistent with the 87 Sr/ 86 Sr data; see below, Fig. 10).
Samples from the two study intervals (Moscovian and Asselian) separate in ordination space and vary with the second RDA axis (Fig. 7). Specifically, elements associated with felsic sources correlate positively with the second RDA axis whereas those associated with mafic sources correlate negatively with the second RDA axis. The Moscovian samples tend to be associated with higher percentages of elements such as Si, K, Zr, and Ti, which are inferred to represent more felsic sources. In contrast, the Asselian samples represent potentially more heavy elements (Co, Mn, Mg, Y, Na, Fe) associated with mafic sources. This suggests a potential provenance shift between the late Moscovian and early Asselian study intervals (Fig. 7).
Furthermore, RDA analysis reveals a significant relationship between high-frequency sequence boundaries (prominent surfaces) and elemental composition of the dust fraction (Supplementary Table S2). Composition of dust in proximity to high-frequency sequence boundaries is more variable than dust composition within sequences (intrasequence, Fig. 7). This indicates a shift in dust provenance between glacial and interglacial phases within both the late Moscovian and early Asselian study intervals. This result is consistent with high Zr/Ti values (independent of grain size; Fig. 8) near the sequence boundaries, indicating a provenance shift rather than in situ weathering. The observed provenance shift could reflect a shift in wind direction and/or strength and/or derivation from expanded source regions as a result of sea level fall and associated exposure of former epeiric sea and shelf regions in siliciclastic-dominant depositional systems.
The Sr isotope signatures of extracted dust are consistent with Arabian-Nubian Shield and southwestern Pangean sources (western Gondwanan, South American; Fig. 10). The provenance of the coarse dust fractions is reinforced by their Sr isotopic compositions ( 87 Sr/ 86 Sr ∼0.710), coinciding with more proximal Arabian crustal sources (Fig. 10). The more radiogenic samples ( 87 Sr/ 86 Sr ∼0.718) contain mostly finer-grained dust and so have two likely, though non-exclusive, explanations. Based on the Sr   isotopic compositions of the source regions, the more radiogenic samples could reflect derivation from more distal southwestern Pangean sources (Fig. 11). Sr isotopic heterogeneity is well known to occur as a function of grain size, which is supported by our grain size data. Sediment "unmixing" during eolian sorting processes, where the finer-grained (clay) fraction concentrates in more distal eolian deposits, results in more radiogenic 87 Sr/ 86 Sr signatures (e.g., as seen in the modern central North Pacific Ocean pelagic clay province; Stancin et al., 2006). Therefore, we argue that all samples containing the finest fraction (<10 μm) measured here have Sr isotopic compositions consistent with several provenance domains (Fig. 10), with most of the isotopic variation being explained either by grain size variations observed across glacial-interglacial cycles, source switching between Arabia and distal southwestern Pangean sources, or a combination of both.

Timing of Dust Influx
Dust flux variation through the sections studied provides an indicator of temporal changes in atmospheric dust loading at both glacial-interglacial and higher-frequency (<10 4 yr) scales that is interpreted to be related to waxing and waning of the southern hemisphere ice sheets controlled by Milankovitch-scale cyclicity. Time series analysis of the dust fraction through the intervals studied supports a hypothesis of orbital forcing with cyclic patterns in dust flux matching frequencies of estimates of Milankovitch cycles for this time period. The stratigraphic pattern of the dust flux indicates minimal flux during highstands (0.19 g/cm 2 /kyr and 0.27 g/cm 2 /kyr in Moscovian and Asselian, respectively) and peak flux during glacial lowstands (3.77 g/cm 2 /kyr and 4.57 g/cm 2 /kyr in Moscovian and Asselian, respectively). The average total mass of dust preserved in intra-cycle strata of the Asselian exceeds that of the Moscovian by 10%; however, at cycle boundaries, the mass of the dust in the Asselian exceeds that of the Moscovian by 30%. Additionally, the majority of dust through the Asselian interval is coarser than that of the Moscovian interval. As the Asselian coincides with the peak of icehouse conditions, the greater volume of the southern ice sheet likely resulted in drier climate in the subtropics . This heightened source aridity combined with greater exposure of shelves, steeper meridional temperature gradients, stronger jet streams in the upper troposphere, and possibly stronger and gustier winds mixing to the surface increased dust loading in the Asselian atmosphere (e.g., McGee et al., 2010;Löfverström et al., 2014).
During the glacial (lowstand) phases, sea-level fall enabled extensive subaerial exposure in the source region (Arabian-Nubian Shield) and subsequently higher atmospheric dust influx to the area studied (Central Persian Terranes). Increasing atmospheric dust loading and grain size coinciding at or near the cycle boundaries affirm a 16-to 19-fold increase in dust deposition during incipient glaciation to full glaciation phases relative to interglacials. In summary, glacial phase dustiness in the Central Persian Terranes could reflect several nonexclusive factors, including greater aridity in emitting regions, decreased vegetation cover and land-area exposed, the possibility of changes in windiness (Supplementary Figs. S8-S10; see footnote 1), and, importantly, increased sediment (especially fines) production driven by glaciation. These large glacial stage increases in dust flux are also associated with increased delivery of highly reactive iron and associated productivity spikes in these records (Sardar Abadi et al., 2020).
Dust cycle modeling broadly confirms and further clarifies this interpretation. During the modeled interglacial interval, a weak dust source is modeled along the coastline to the west of the site (Fig. 11A), but the main substantive sources of sub-10 μm dust lie outside the region in western Gondwana (e.g., Brazil, Bolivia) Marshall et al., 2016;Carvajal et al., 2018). Sub-10 μm dust deposition at the site, as found by interpolation, is 0.30 g/cm 2 /kyr, within a factor of 1.6 of estimated dust deposition at the site during interglacial highstands (0.19-0.27 g/ cm 2 /kyr) (Fig. 11C), when the deposition flux is almost entirely below 10 μm in size. During the glacial interval modeled, drier conditions (mean annual precipitation <400 mm) advanced into eastern sub-tropical Gondwana, making central African terranes a strong dust source and strengthening dust emission along the coast (corresponding to the Arabian-Nubian shield) by two orders of magnitude or more (Fig. 11B), but simulated average winds did not increase ( Supplementary Figs. S8-S10). Thus, the model simulates the activation of proximal coastal dust sources primarily attributable to increased aridity/decreased vegetation that can explain the influx of larger particles at glacial lowstands. Dust deposition at the site is simulated to be 1.94 g/cm 2 /kyr (Fig. 11D), a factor of 2 or more less than dust deposition estimated at the site during glacial lowstands (3.77-4.57 g/cm 2 / kyr). Given that deposition at glacial lowstands remained dominated by sub-10 μm particles, the Figure 11. Simulated dust emission and deposition (g/cm 2 /kyr) near the study site for dust cycle simulations during characteristic interglacial and glacial intervals are shown. The dotted white line indicates the 30°S latitude line, while the solid white line indicates the coastline in the simulation paleogeography. A black circle denotes the approximate location of the study site. Gray Xs indicate land areas that experienced mean annual precipitation less than 400 mm. The approximate location of the Arabian-Nubian Shield is labeled.

A B C D
dust deposition at the site in the 0.1-10 μm size range tracked by the model is approximately equal to the total dust deposition. Thus, the modeling underestimates the dust deposition by a factor of ∼2. The glacial/interglacial ratio in dust deposition at the site is 16-19 but only 6.5 in the model (Fig. 12). The underestimate of glacial-interglacial variability is initially surprising. The underlying baseline climatic simulations span a much wider dynamic range of forcing than is expected for glacial-interglacial variability during the Moscovian or the Asselian. The interglacial simulation has pCO 2 of 2500 ppmv and no land ice, while the glacial simulation has pCO 2 of 250 ppmv and extensive high-latitude ice sheets in Gondwana. While glacial-interglacial pCO 2 variability was much higher during the late Paleozoic than during the late Cenozoic (Montañez et al., 2016), it seems to have been closer to a factor of 3 during the Moscovian than a factor of 10.
However, it is well-known that dust cycle models tend to be stiff and require significant tuning, particularly to represent non-meteorological processes that affect dust transport such as glacial grinding of geologic materials, exposure of glacial flour, and perhaps other aspects of landscape evolution (e.g., Mahowald et al., 2006;Sugden et al., 2009). The tuning approach used by Marshall et al. (2016) (and adopted here) accounts for meteorological and vegetation change between interglacial and glacial intervals but has no way to account for changes in dust availability in dust source regions. Thus, the close agreement in simulated and measured dust deposition at the site during interglacial intervals suggests the coal-based deposition estimate used to tune the simulations of Marshall et al. (2016) was representative of interglacial conditions, while the discrepancy in reconstructing the amplitude of glacial-interglacial variability at this site suggests that dust sourced from dynamic glaciation on Milankovitch timescales may have contributed the vast majority of enhanced dust deposition recorded at glacial lowstands.

CONCLUSIONS
Dust extracted from shallow marine carbonate of the Central Persian Terranes of Iran along the northeastern Gondwanan margin archives changes in atmospheric dust loading and circulation during the middle Pennsylvanian-earliest Permian in a region relatively proximal to central Arabia and eastern African glacial regions. The carbonate strata occur in repetitive cycles with boundaries demarcated by either abrupt deepening or subtle signs of subaerial exposure and represent deposition in response to the waxing and waning of the Gondwanan ice sheets. The volume (wt%) and grain size distribution of dust particles vary systematically, and both increase during glacials. Orbital forcing is evident in the record of dust flux across the study interval. The estimated mass accumulation rate of dust particles records relatively minimal flux during interglacial highstands (0.19-0.27 g/cm 2 / kyr) and maximum flux during glacial lowstands (3.77-4.57 g/cm 2 /kyr), indicating that glacial dust flux exceeded that of interglacials by a factor of 16-19. Glacial phase increases in dust flux were likely driven primarily by (1) enhanced aridity and land exposure accompanying glacial conditions and (2) increased production of dust by mechanisms such as glacial grinding during glacial phases. Glacial phase aridity acted to mobilize dust from wider regions, including glacial outwash plains of the Arabian-Nubian Shield (proximally) and regions west in southwestern Gondwana (e.g., South America). The high amplitude of glacial-interglacial variability in dust deposition indicates the importance of the resupply of large amounts of dust to glacial outwash regions by dynamic advance and retreat of ice sheets on Milankovitch timescales. Finally, the greater, overall dust flux of the Asselian interval relative to that of the Moscovian interval indicates greater atmospheric dustiness during the Asselian, which is consistent with the interpretation of the Asselian as the most intense interval of the late Paleozoic ice age.