Upper Devonian and Lower–Middle Mississippian strata of the North American midcontinent are ubiquitously fine-grained and silt-rich, comprising both so-called shale as well as argillaceous limestone (or calcareous siltstone) that accumulated in the Laurentian epeiric sea. Although long recognized as recording marine deposition, the origin and transport of the fine-grained siliciclastic material in these units remains enigmatic because they do not connect to any proximal deltaic feeder systems. Here, we present new data on grain size, whole-rock geochemistry, mineralogy, and U-Pb detrital-zircon geochronology from units across Oklahoma; we then integrate these data with models of surface wind circulation, refined paleogeographic reconstructions, and correlations from the greater midcontinent to test the hypothesis that wind transported the siliciclastic fraction to the marine system. The exclusively very fine silt to very fine sand grain size, clear detrital origin, widespread distribution over large regions of the epeiric sea, Appalachian sources, and paleogeographic setting in the subtropical arid belt far-removed from contemporaneous deltaic feeder systems are most consistent with eolian transport of dust lofted from subaerial delta plains of the greater Appalachian orogen and incorporated into subaqueous depositional systems. Delivery of dust that was minimally chemically weathered to Devono-Mississippian epeiric seas likely provided essential nutrients that stimulated organic productivity in these commonly organic-rich units.

Siltstone, argillaceous siltstone, silty carbonate, siliciclastic mudstone, shale, and chert are common lithofacies in the upper Devonian–Middle Mississippian strata in basins of the North American midcontinent that accumulated in epeiric seas of the southern Laurentian craton (Fig. 1, 2; Lowe 1976; Lane and De Keyser 1980; Gutschick and Sandberg 1983; Over 1990; Schwartzapfel 1990; Turner et al. 2015; Schieber 2016; Kondas et al. 2018; Godwin et al. 2019; Mazzullo et al. 2019; Miller et al. 2019b; Milad et al. 2020; Price et al. 2020). These strata are commonly organic-rich, and have attracted interest primarily as petroleum source rocks and unconventional reservoirs (Barrows and Cluff 1984; Smith and Bustin 1998; Smith and Bustin 2000; Slatt and O'Brien 2011; Angulo and Buatois 2012; Higley 2013; Wang and Philp 2019), but they also provide insights on climate, tectonics, and carbon sequestration of the time. The range of explanations for deposition of the siliciclastic fraction includes: 1) suspension settling (Kirkland et al. 1992; Smith and Bustin 1998), 2) a combination of suspension settling, density currents, and gravity flows (e.g., Wright and Friedrichs 2006; Loucks and Ruppel 2007; Macquaker et al. 2010; Milad et al. 2020; Price et al. 2020), and 3) wind- or tide-induced subaqueous currents that remobilized muds across ancient epeiric seas (Schieber 2016). Although all accept that these strata ultimately accumulated in marine environments, the provenance, generation, and transport of the significant volume of uniformly fine-grained siliciclastic material remain enigmatic and unresolved.

To better constrain the origin and transport of the siliciclastic fraction in these units, this work integrates new data from Upper Devonian and Lower–Middle Mississippian siltstone, mudrock, and argillaceous carbonate units across Oklahoma with compilations of coeval systems of the greater midcontinent region (Fig. 2). We use these data to test the hypothesis of eolian delivery to the marine environment for the vast volume of this material. An eolian-dust origin explains otherwise perplexing attributes of these units and has significant implications for interpretation of organic productivity in these and analogous systems.

The Anadarko Basin was initially part of the Oklahoma Basin (Johnson et al. 1988) which originated in the Early Cambrian in response to thermal subsidence associated with development of the Southern Oklahoma Aulacogen (Wickham 1978; Dalziel 1991; Meert and Torsvik 2003). Cambro-Ordovician carbonate strata accumulated in the axis of the failed rift (Perry 1989) as a broad epicontinental shelf developed on the flanks (Johnson et al. 1989). Carbonate and, subordinately, siliciclastic mudstone sedimentation continued through the Devono-Mississippian (e.g., Hunton Limestone, Woodford Shale, Sycamore Limestone, Caney Shale), a time of relative tectonic quiescence in this part of the Anadarko Basin and greater midcontinent region, before the tectonic disruptions of the Ouachita and Ancestral Rocky Mountains orogenies (Feinstein 1981; Johnson et al. 1988, 1989). Closure of the Rheic Ocean and associated contractional tectonics beginning in the latest Mississippian (Serpukhovian)–Early Pennsylvanian (Nance et al. 2010) ultimately produced the thin-skinned Ouachita fold-and-thrust belt and the associated Arkoma basin (Fig. 1). By Pennsylvanian time, a series of NW–SE-oriented basement-cored uplifts developed, termed the Ancestral Rocky Mountains, including the Wichita–Amarillo uplift and yoked Anadarko Basin of southwestern Oklahoma and northern Texas (Kluth and Coney 1981). The present configuration of the Anadarko Basin is a product of this compressional inversion of the earlier aulacogen along basement-involved reverse faults, and development of accommodation space along the mountain front (Johnson et al. 1989).

Beyond Oklahoma, Late Devonian orogenic events affected both the eastern and western margins of the craton, in the form of the Acadian–Neoacadian and Antler orogenies, respectively. The Acadian and Neoacadian orogenies produced a series of westward-prograding clastic wedges, the Catskill Delta complex during the Middle–Late Devonian and the Price–Pocono–Borden–Grainger delta complex during the Early–Middle Mississippian (Ettensohn and Woodrow 1985; Ettensohn 2004, 2008, 2022), whereas the Antler clastic wedge prograded eastward (Sandberg and Gutschick 1980; Johnson and Pendergast 1981) (Fig. 3). The Oklahoma Basin remained isolated and far removed from each of these clastic wedge complexes, ∼ 1800 km from the Antler and ∼ 1500 km from the Appalachian, throughout the Late Devonian to Middle Mississippian (Fig. 3). Note that we use North American stage names for the Mississippian, for consistency with the literature.

The eastern deltaic complexes are well exposed in outcrops and boreholes throughout the Appalachian Valley and Ridge and Appalachian Plateau provinces, and exhibit both marine and continental facies for nearly 1000 km along orogenic strike (Fig. 1). The Catskill Delta complex developed during the Middle to Late Devonian, transporting siliciclastics from the Acadian Orogen to marginal marine environments along the Appalachian foreland (Fig. 3). During the Famennian, this complex comprised an estimated four main delta systems, fed by low-gradient coastal rivers emanating from alluvial systems in the highlands—akin in scale to the Eocene and modern Texas–Louisiana coastal plain, or the northern Sumatran coast (Sevon 1985; Woodrow et al. 1985; Boswell and Donaldson 1988). Organic-rich black shales characterize the basinal facies of the Catskill Delta complex in the Appalachian foreland basin (Ettensohn and Woodrow 1985; Boswell and Donaldson 1988; Ettensohn 2008; Ettensohn et al. 2019).

Clastic-wedge deposition continued into the Early–Middle Mississippian (Kinderhookian–Osagean) in response to Neoacadian tectonism, such that the Appalachian foreland overfilled and deltaic deposition of the Borden and Grainger formations (Figs. 2, 3) prograded into the eastern and southeastern Illinois Basin (Ettensohn et al. 2022). Active denudation ceased in the Middle to Late Mississippian (Meramecian–middle Chesterian), resulting in extensive carbonate deposition in the Appalachian foreland (Greenbrier and Newman Limestones, Fig. 2; Ettensohn et al. 2022). Toward the eastern Appalachian foreland basin, carbonate strata intertongue with continental redbeds and eolianites throughout West Virginia, Virginia, and Maryland (Wynn and Read 2008; Ettensohn et al. 2022). Early development of the Mauch Chunk Delta began in eastern Pennsylvania in the Meramecian (Wynn and Read 2008; Ettensohn et al. 2022) (Fig. 3).

In the extensive cratonic region between the Acadian (east) and Antler (west) orogenic systems, broad epeiric seas, intracratonic basins, and emergent arches prevailed through the Late Devonian–Early Mississippian, with organic-rich mudstones, including the Chattanooga Shale, Ohio Shale, New Albany Shale, Antrim Shale, and Woodford Shale accumulating east of the Transcontinental Arch and correlative mudstones such as the Bakken Shale and Duvernay Shale accumulating west of the Transcontinental Arch (Figs. 2, 3). Carbonate-dominant systems with variably silty and cherty components prevailed through the Early–Middle Mississippian (e.g., the North Hill Limestone, Keokuk Limestone, St. Joe Group, Boone Group, Osage Limestone of northern Oklahoma, and Madison Limestone; Figs. 2, 3; Gutschick and Sandberg 1983, and others). By the Middle to Late Mississippian (Meramecian–late Chesterian), fine siliciclastic contents increased, with deposition of silty carbonates, marine siltstones, and silty organic-rich mudstones (e.g., Blue River Group, Warsaw Shale, Barnett Shale, Meramec and Sycamore formations, Caney Shale and Madison Limestone; Figs. 2, 3).

The midcontinent region ranged from southern subtropical latitudes (∼ 15° S–30° S) in the Late Devonian to near-tropical latitudes (∼ 15° S) by the Early–Middle Mississippian (e.g., Boucot et al. 2013; Domeier and Torsvik 2014; Scotese and Wright 2018) (Fig. 3). Paleoclimatic reconstructions document arid conditions in lowland regions throughout the midcontinent within these time intervals (Boucot et al. 2013) (Fig. 3). Additionally, a number of instances of ice-contact deposits document evidence for suspected icesheets and/or alpine and even tidewater glaciation within the central Appalachian orogenic system during the latest Devonian (approximately lower to middle Siphonodella praesulcata conodont zonation; LE, Retispora lepidophyta–Indotriradites explanatus to LN, Retispora lepidophyta–Verrucosisporites nitidus miospore zonation) (e.g., Spechty Kopf Formation, Isaacson et al. 2008; Brezinski et al. 2008, 2009, 2010; Lakin et al. 2016; Ettensohn et al. 2020a, 2020b) (Fig. 2).

Glaciation in the central Appalachians marks the apex of the late Devonian icehouse. Glaciation in Gondwanaland began in the late Famennian (approximately Palmatolepis postera conodont zonation; VCo, Diducites versabilis–Grandispora cornuta miospore zonation) and extended to the latest Famennian (Upper Siphonodella praesulcata conodont zonation; LN, Retispora lepidophyta–Verrucosisporites nitidus miospore zonation) (Fig. 2) (Caputo et al. 2008; Isaacson et al. 2008; Lakin et al. 2016). Additionally, minor evidence for Gondwanan glaciation exists for the Early and Middle Mississippian (Kinderhookian–Osagean and Meramecian–Chesterian) (Caputo et al. 2008; compilation in Soreghan et al. 2019).

In Table 1, we summarize the facies attributes and inferred depositional environments for Upper Devonian and Lower–Middle Mississippian strata of Oklahoma (Figs. 2, 4, 5). These units are well exposed and have been widely characterized in the Arbuckle Mountains (southern Oklahoma), the Ozark Plateau (northeastern Oklahoma), and the subsurface; correlatives occur throughout the greater midcontinent as demonstrated through hundreds of thousands of well penetrations (Craig et al. 1978; See Supplemental File 1).

Paleogeographic Reconstructions

To place the study units (Fig. 4) in a broader context (Figs. 2, 3), we used the paleogeographic maps of Blakey (2013), North American Key Time Slices ©2013 Colorado Plateau Geosystems, Inc, as base maps, and appended these with major climate belts following Boucot et al. (2013) and Scotese and Wright (2018). The Correlation of Stratigraphic Units of North America (COSUNA) project charts (Ballard et al. 1983; Hills and Kottlowski 1983; Patchen et al. 1984a, 1984b; Bergstrom and Morey 1985; Hintze 1985; Shaver 1985; Adler 1986; Mankin 1987; Kent et al. 1988) provided the initial framework for age constraints and correlations, but all unit ages were updated with the most recent available data (Cohen et al. 2013) (See Supplemental File 2 for references).

Grain-Size Analysis

To assess the nature of and possible transport modes for the siliciclastic material in the target strata, we collected samples for particle-size analysis from seven outcrops and ten cores distributed throughout Oklahoma (Fig. 4; See Supplemental Files 2 and 3). While carbonate grains and iron oxides, as well as grain coatings, may constitute a minor fraction of the original detrital material, we specifically target the silicate mineral fraction (SMF) as a proxy for detrital grain size. Therefore, all organic, carbonate, and iron oxide components were removed using a sequential procedure following standard methods (cf. Sur et al. 2010a, 2010b; Jiang and Liu 2011; Sweet et al. 2013; Foster et al. 2014), as outlined briefly here.

Samples were gently crushed with a ceramic mortar and pestle to pea-size gravel to accelerate chemical reactions, then rinsed with distilled water and sieved at 250 μm to remove any fines generated during crushing. Samples were subsequently treated with 2N HCl (minimum 24 hours) for carbonate removal. To remove organic matter and oxidize pyrite, samples were combusted at 500 °C for ∼ 12–20 hours. This does not affect the grain size of silicate minerals, since silicates do not begin to sinter until temperatures exceed 500 °C (Schomburg 1991; Rickard et al. 2016). Samples were then treated with sodium citrate–bicarbonate–dithionite (CBD) (Sur et al. 2010a; Rea and Janecek 1981) rinses until all iron oxides were removed.

Smear slides of residues were inspected using a polarizing reflected-light microscope under 10×, 20×, and 40× objectives to assess sample composition and ensure disaggregation, and check for the presence of remnant authigenic phases such as silicified microfossils (Fig. 6). If present, silicified microfossils were physically removed by sieving following visual inspection to determine appropriate sieve size needed to remove microfossils, microfossil fragments, and/or aggregates while retaining detrital material (See Supplemental Files 2 and 3 for additional information). Samples exhibiting a predominance of authigenic material were eliminated from analysis. Furthermore, both visual inspection and previous studies on these strata indicate no to minimal instances of quartz overgrowths (Shelley et al. 2019; Price et al. 2020; Milad et al. 2020; Duarte et al. 2021; Supplemental Files 3, 4), lending further confidence to the fidelity of subsequent particle-size analyses. Particle sizes of fully disaggregated samples were measured using a Malvern Mastersizer 3000 Laser Particle Size Analyzer (LPSA). Of 103 samples processed in this way, 82 disaggregated sufficiently to enable analysis by LPSA.

Geochemical Analysis

To address compositional characteristics of the silicate material in the Sycamore and Woodford Shale, nineteen ∼ 5 g samples from the Arbuckle Uplift transect (Figs. 2, 4) were fragmented to gravel size and submerged in a 1N HCl for 24 hours to eliminate carbonate before whole-rock and trace-element geochemical analysis (acquired by ALS Geochemical). Remnant CaO values and high loss-on-ignition (LOI) necessitated application of the correction by Taylor and McClennan (1985) to interpret the silicate fraction. Data were analyzed (and biplots were created) using principal-components analysis (PCA) in R.

Detrital-Zircon Geochronology

To assess provenance of the siliciclastic material, three 10–20 kg samples from each of the Arbuckle Uplift units (Figs. 2, 4; Woodford Shale, Sycamore Limestone, and Caney Shale) were processed for zircon separation and analyzed for detrital-zircon geochronology at the University of Arizona (UA) LaserChron Center. Analyses were conducted using the Element 2 laser ablation–single collector–inductively coupled plasma–mass spectrometer (LA-SC-ICP-MS) following the analytical and data processing methods described in Pullen et al. (2018) and Gehrels and Pecha (2014). Insufficient yield from the Woodford Shale sample (despite > 20 kg of processed material) precluded analysis. Best ages from the Sycamore Limestone and Caney Shale are reported as 206Pb/207Pb for ages greater than 900 Ma and 206Pb/238U for ages less than 900 Ma. These ages are displayed as normalized probability density plots using an Excel macro, provided by the UA LaserChron Center (Gehrels 2010). Analyses with > 10% discordance and 5% reverse discordance were eliminated before analysis.

General-Circulation Modeling

To address climatic attributes and potential atmospheric transport pathways, we reanalyzed results of a total of 12 coupled general-circulation-model (GCM) experiments of global climate from intervals spanning the Middle Devonian to the Late Mississippian. The simulations were performed with the HadCM3BL-M2.1aD atmosphere–ocean–vegetation model (Valdes et al. 2017) and are described in detail in Valdes et al. (2021). The model uses the same horizontal grid with a resolution of 3.75° × 2.5° in the atmosphere and ocean, with 19 and 20 unequally spaced vertical levels, respectively. The paleogeographies underlying the model simulations were derived from the PALEOMAP paleogeographic atlas (Scotese and Wright 2018). Current geographical coordinates of Oklahoma and the Appalachian clastic wedges were rotated using the PALEOMAP rotation model in the pyGplates software (Müller et al. 2018) to place these features within the appropriate paleogeographical positions. Locations for the Appalachian clastic wedges are from Ettensohn et al. (2019). Besides the model geography, the simulations differ only by a linearly increasing solar constant (0.8% per 100 Myr) and varying concentrations of atmospheric carbon dioxide. CO2 levels for the selected simulations were derived from the Foster et al. (2017) curve and range from a peak of 1377 ppmv in the Middle Devonian to a minimum of 233 ppmv in the Middle Mississippian. Model performance, the simulated climate evolution throughout the Phanerozoic and comparison to proxy data are presented in Valdes et al. (2021).

Analyses of Composition and Grain-Size

Results

Silicate material from the Woodford Shale contains silt-size quartz, mica, and feldspar (and trace accessory phases), and minor clay. Visual inspection of the coarsest (> ∼ 50 μm) fraction revealed silicified algal cysts (Tasmanites) and, in some cases, mudstone aggregates, both of which were removed before LPSA analysis. Grain-size modes typically range from ∼ 5 to 12 μm (n = 12). Spatially, the distribution of grain sizes appears generally uniform with a slight coarsening toward north-northwest Oklahoma.

Silicate material from the Mississippian units (Meramec, Sycamore Limestone, Caney Shale, Bachelor Formation, St. Joe Group, and Boone Group) includes quartz, feldspar, biotite, and metamorphic lithic fragments, as well as trace accessory phases (Fig. 6). In general, the older (Kinderhookian–Osagean) units (i.e., Welden Limestone, Bachelor Formation, Northview Formation, Compton Formation) are finer grained (modes of 5–17 μm; two samples from Skaggs core, 34–71 μm) relative to the younger (Meramecian–Chesterian) units (7–83 μm) (Fig. 7; Supplemental File 2).

Interpretation

The presence of quartz, feldspar, mica, lithic fragments, and accessory silicates (e.g., zircon) confirm the substantial detrital contributions to the study strata (Fig. 6). The Upper Devonian and Lower Mississippian units exhibit the finest modes, whereas Meramecian—Chesterian units are coarser (Fig. 7). The Meramec and Sycamore sections are the siltiest strata by volume of the studied units (Supplemental File 5). All units have modes falling between very fine silt and very fine sand fractions. Such uniformly fine modes are consistent with the hypothesis of eolian transport for this material, addressed further below.

Geochemistry

Results

Principal-component-analysis (PCA) biplots of major oxides and trace elements (Fig. 8) indicate the relative compositional variability within the dataset depicted by the vectors (Greenacre 2010). The Woodford Shale and Sycamore Limestone occupy distinctly different fields (Fig. 8A), with the Sycamore Limestone displaying significantly greater compositional variability than the silica-rich Woodford Shale. Relative to the North American Shale Composite (NASC) and contemporaneous Appalachian Shale (e.g., Sunbury Shale; Rimmer 2004), samples from both the Woodford Shale and the Sycamore Limestone exhibit depletion in Al2O3 but enrichment in SiO2 (Fig. 8B). From the data illustrated in Figure 8, the siliciclastic fraction of the Sycamore Limestone is compositionally more similar to loess (both modern loess and ancient loessite) than to the NASC (Fig. 8B), especially in siltier layers that compose the upper part of the formation (avg. 85% SiO2; 6% Al2O3; Chemical Index of Alteration of 57–62). Elemental ratios of generally insoluble elements (Zr and Ti, Al oxides) are plotted vs. depth (Fig. 8C) for the Devonian–Mississippian Arbuckle location (Fig. 2, 4). For all three elemental ratios (Zr/TiO2, Zr/Al2O3, and TiO2/ Al2O3), values increase abruptly in the lower Sycamore, then return to values comparable to the Woodford but gradually increase again through the middle–upper Sycamore (Fig. 8C).

Interpretation

The geochemical data confirm the mineralogical evidence for detrital contributions to both the Sycamore Limestone and Woodford Shale. Intriguingly, neither unit overlaps with the North American Shale Composite or the Appalachian Shale, both of which exhibit much higher values of Al more typical of source material that is highly chemically weathered. Rather, both units—especially the Sycamore Limestone—exhibit a signature more similar to loess (Fig. 8B). The inflection in Zr and TiO2 with respect to Al2O3 (Zr/Ti O2, Zr/Al2O3, TiO2/Al2O3) in the basal Sycamore Limestone (Fig. 8C) is interpreted to reflect a possible excursion in provenance. These compositional data are independent of grain size, as evinced by the weak relationship (r2 < 0.07) between grain size and these elements (Supplemental File 4), as well as the negligible grain-size variation within each of the Sycamore and Woodford formations in general (Fig. 7).

Detrital-Zircon Geochronology

Results

U-Pb detrital-zircon ages from the Sycamore Limestone (N = 1; n = 218) and Caney Shale (N = 1; n = 215) yielded a total of 433 concordant to slightly discordant Paleozoic to Archean (age) grains (Fig. 9). Detrital zircon grains recovered for analysis from the Sycamore Limestone are well sorted and considerably finer than 100 μm, while zircons recovered from the Caney Shale are moderately sorted with a majority of the grains finer than 100 μm. Overall, the age spectra for both the Sycamore and Caney units reveal similar age distributions, with detrital zircon ages ranging from 365 ± 5 Ma to 3022 ± 12 Ma. The Sycamore Limestone has a primary population of Mesoproterozoic ages (56%), with a primary peak at 1047 Ma and subordinate peaks at 1184, 1382, and 1463 Ma. Secondary populations in the Sycamore unit include Paleoproterozoic ages (16%) with peaks at 1640 Ma and 1739 Ma, as well as early–mid Paleozoic ages (14%) with primary peaks at 449 Ma and 473 Ma, and subordinate peaks at 363, 395, and 507 Ma. The Caney Shale also has a primary detrital-zircon population of Mesoproterozoic ages (38%) with primary peaks at 1062 Ma and 1356 Ma, and secondary populations of Paleoproterozoic ages (26%) with a primary peak at 1817 Ma subordinate peak at 1678 Ma, and Archean ages (23%) with a primary peak at 2715 Ma (Fig. 9).

Interpretation

Potential source terranes for the siliciclastic fractions of these units include the Archean Superior and Wyoming provinces (> 2.5 Ga) in northern Laurentia, the Trans-Hudson and Penokean provinces (1.8–1.9 Ga) in north-central Laurentia, the Mazatzal–Yavapai (1.6–1.8 Ga), Granite–Rhyolite (1.3–1.5 Ga), and Grenville provinces (0.9–1.3 Ga) throughout central Laurentia and Peri-Gondwanan basement rocks, Pan-African Gondwanan terranes, Taconic plutons, Acadian plutons, and associated synorogenic strata (365–820 Ma) along the Appalachian suture zone in eastern Laurentia (e.g., Hoffman 1989; Van Schmus et al. 1993; Becker et al. 2005; Whitmeyer and Karlstrom 2007; Thomas et al. 2004). The Sycamore Limestone is interpreted to have primary detrital-zircon contributions from the Grenville Province (41%), with secondary contributions from the Granite–Rhyolite province (18%) and Acadian and Taconic plutons from the Appalachian orogen (18%), with minor contributions from the Mazatzal–Yavapai province (13%), Archean provinces (7%), and Trans-Hudson–Penokean provinces (3%; Fig. 9). The Caney Shale is also interpreted to contain primary detrital-zircon contributions from the Grenville Province (33%), but it has secondary contributions from Archean age populations (Wyoming and Superior cratons, or recycled early Paleozoic midcontinent sandstones; 23%), and minor contributions from the Granite–Rhyolite province (12%), Mazatzal–Yavapai province (11%), Appalachian terranes (∼11%) and Trans-Hudson–Penokean provinces (∼ 9%).

Both the Sycamore and Caney units contain abundant Grenville-age grains (33–40% of total age population) with primary peaks at 1047 Ma (Sycamore) and 1062 Ma (Caney), which fall in the range of Appalachian Grenville sources (Thomas et al. 2017). Age-equivalent Appalachian foreland strata from Park et al. (2010) also yield abundant Grenville-age detrital zircons with a peak age of 1048 Ma (Fig. 10). In addition, both the Sycamore Limestone and the Caney Shale show contributions (∼ 11–18%) from the Taconic and Acadian orogens in the northern Appalachian suture zone, represented by detrital-zircon ages ranging ca. 490–440 Ma and 420–350 Ma (e.g., Thomas et al. 2004, Ettensohn et al. 2019). The Caney Shale did yield significant Archean-age grains, potentially sourced from recycled lower Paleozoic sandstone units from the northern midcontinent and northern Laurentia (Konstantinou et al. 2014).

Simulated Surface Winds

Results

The paleoenvironmental conditions of the Oklahoma and greater North American mid-continent region changed from predominantly subtropical and arid conditions during the Famennian to a rather monsoonal system exhibiting seasonal wind reversals and summer precipitation during the Early to Middle Mississippian (Kinderhookian through Meramecian). These seasonal wind reversals suggest that the greater midcontinent region was influenced by trade winds during Southern Hemisphere winter months (June through August) (Fig. 10). The simulated surface wind patterns are dominated by an interplay between the prevailing southeasterly trade winds, the land–sea thermal gradient, and local orography. As such, Acadian and Caledonian orogens tend to, in general, guide the flow of air around the landmasses, and absolute wind speeds are reduced due to the presence of these orogenic belts.

Summer winds in the Late Devonian (Fig. 10A) are shown to have come from the south and diverge across the Laurentian epeiric seas with a westerly jet flowing along the northern front of the Acadian orogen. During Late Devonian winter months, a more southeasterly trade-wind circulation developed across Laurentia (Fig. 10B). Westerlies along the southern front of the Acadian orogen show reduced wind values topping over the mountain range with acceleration on the lee side of the orogen. The Osagean (Fig. 10C, D) and Meramecian (Fig. 10E, F) surface-wind patterns are very similar. With the progressive northward movement of Laurentia, the mid-continent region came under the influence of cross-equatorial, northwesterly winds bringing moist air and summer precipitation toward emergent landmasses (Fig. 10E, F). Winter precipitation was very low over the adjacent subtropical land masses and limited to tropical regions northwards of 10° S throughout the Devonian to Mississippian (Supplemental File 6). Surface windflow patterns during the Osagean and Meramecian winter months (JJA) flowed westward along the southern front of the orogen, topped over the range, and swept across the Laurentian epeiric seas in a west to northwestern direction (Fig. 10E, F; Supplemental File 6).

Interpretation

Model simulations reveal that the Laurentian landmass and associated orogens strongly influenced surface wind patterns. In fact, the presence of a larger landmass at subtropical latitudes during the Early Mississippian promoted seasonal wind reversals and a southward migration of the Intertropical Convergence Zone (ITCZ). According to the model data, a northwestward transport of continental dust toward the midcontinent seems physically plausible, although such a transport would have been possible only during the austral winter months. These months are characterized by northwestward, offshore winds across the potential source regions of the Acadian–Neoacadian clastic wedges.

These wind patterns may be analogous to present winds along the Andean orogen, where westerlies from the South Pacific anticyclone are deflected northward creating a southerly jet along the Andean Pacific front. Air masses top over the orogen and descend as warm dry air off the lee side of the mountain range, known as Zonda winds (Puliafito et al. 2015; Zou and Xi 2021) which are similar in process to Foehn winds off New Zealand's Alps, Chinook winds off the North American Rocky Mountains, and Santa Ana winds over the mountains of southern California (McGowan et al. 1996; Hugenholtz et al. 2007; Álvarez et al. 2021).

Even though coastal and continental summer precipitation increased towards the Early Mississippian, winter months around the orogenic highlands were still dry in the simulations and a stronger easterly wind component may have transported material to the Oklahoma region of Laurentia. A comparable analog for the Early–Middle Mississippian model results, which sits at the same latitude, might be present-day northwestern Australia, where easterlies transport dust westward towards the Indian Ocean during the winter months following the northward migration of the ITCZ (Christensen et al. 2017; Keep et al. 2018).

The Enigma of the Missing Proximal Deltas

Delivery of siliciclastic material to epeiric seas occurs by either deltaic or eolian transport. Various processes can then redistribute this material, e.g., hyperpycnal flows, sediment gravity flows induced by slope failure or wave and/or current agitation, or contour currents (Stow et al. 2001; Walsh and Nittrouer 2009; Macquaker et al. 2010; Schieber 2016). Assuming an average wave base of 50 m, Schieber (2016) posited a maximum transport distance of 100 km for remobilized sediment, or up to 400 km for an extreme case of a 200 m wave base. Acknowledging the vast distribution of siliciclastic mud across the Devonian Laurentian sea, Schieber (2016) advocated transport induced by bottom currents, but these arguments demand lateral transport for thousands of kilometers in low-gradient epeiric seas with wave bases << 200 m. While we acknowledge these arguments as one hypothesis, we propose that eolian delivery is an additional and simpler explanation that avoids the need to call upon unknown and unusually powerful storms or tides, as well as flocculation, which is inconsistent with the minimal clay-mineral content of the study units.

The genesis and initial mobilization of the large volumes of ubiquitously fine siliciclastic material characterizing the study units remain unresolved, in large part owing to the absence of proximal feeder system(s). For the Devono-Mississippian Bakken Formation of the Williston Basin, for example, Egenhoff and Fishman (2013) inferred the presence of a deltaic fan that is either not preserved or not exposed. Similar inferences of deltaic feeder systems are applied to Devono-Mississippian mudstone units of the midcontinent (Price et al. 2020) yet without documentation of any remnants of the subaerial (delta plain) or proximal (subaqueous delta) depositional systems. It is highly unlikely that researchers have overlooked such delta systems, given the great density of subsurface data points across the midcontinent region (e.g., > 100,000 wells in Kansas and Oklahoma; Newell et al. 1987; Oklahoma Corporation Commission 2021). The coeval deltas nearest to the midcontinent region during Late Devonian time are the Catskill (Appalachian foreland) and Antler (Antler foreland) systems, and, for the early Mississippian, the Price–Pocono–Grainger–Borden (Appalachians) deltas, that are separated from the midcontinent study area by ∼ 1,000–1,500 km (Fig. 3) (Lane and De Keyser 1980; Ausich and Meyer 1990; Gutschick and Sandberg 1991; Richardson and Ausich 2005; Ettensohn 2004, 2008; Ettensohn et al. 2019, 2022).

In contrast to the Upper Devonian, Lower–Middle Mississippian (Kinderhookian–Meramecian) strata in northeastern Oklahoma are cast as a shallow-water carbonate ramp that yielded to the sediment-starved (during Kinderhookian–Osagean) Ouachita trough (southeast) and proto-Anadarko Basin (southwest; Lane and De Keyser 1980; Gutschick and Sandberg 1983; Godwin and Puckette 2019; Miller et al. 2019b). Yet the occurrence of partial Bouma sequences inferred from the Sycamore Limestone (southern Oklahoma) (Schwartzapfel 1990; Miller and Cullen 2018; Milad et al. 2020), as well as clinoform geometries inferred from subsurface data, led Price et al. (2020) to similarly infer a nearby fluvio-deltaic source despite the ubiquitous absence of evidence for proximal feeder systems. The Transcontinental Arch and the Ozark Dome were the most proximal exposed regions during this time but consist of correlative carbonates with limited exposures of Cambro-Ordovician sandstones and localized occurrences of upper Devonian sandstone—insufficient in volume, grain size, and provenance to supply the enormous amounts of exclusively fine-grained material to surrounding regions (Hills and Kottlowski 1983; Bergstrom and Morey 1985; Adler 1986).

If any of these Devono-Mississippian systems were deltaic, they are uniformly “headless”—missing all but the distal-most “prodelta” fines, even where the shelf margins are preserved. In other words, there were no deltaic systems located nearby to transport the siliciclastic material basinward. And if the siliciclastic units represent prodelta fines, they exhibit remarkably uniform grain sizes that fail to fine with distance from the inferred deltaic feeders, and exhibit an unusual dearth of true clay (minerals or grain sizes; Figs. 7, 8).

Origin of Voluminous Mid-Paleozoic Silt

In contrast to the prevalence of fluvio-deltaic inferences for the study units, few have mentioned eolian processes (e.g., Miller and Cullen 2018; Shelley et al. 2019), whilst others make a connection to the origin of chert in correlative units such as the Arkansas Novaculite, the Caballos Novaculite, and cherts of the Osage–Meramecian strata (Banks 1970; Lowe 1976; Cecil 2004, 2018b). Schieber (2016) explicitly dismissed an eolian origin for these mudrocks based on the presumed need for large erg systems to supply the fines via saltation abrasion of sand; yet eolian abrasion produces very little silt (Swet et al. 2019, 2020; Adams and Soreghan 2020). Dust-emitting regions today include both low-subtropical latitude arid regions such as dry lakes of the Bodélé Depression that supply diatomaceous-rich dust (West Africa; Washington et al. 2006; Chappell et al. 2008) and the Nile Delta plain, which fuels the dunes and loess of the Sinai and Negev (Muhs et al. 2013), as well as high-latitude proglacial and periglacial regions, such as the Copper River delta, which seasonally emits large volumes (60–160 ktons in 2006; Crusius et al. 2011) of glacially ground rock powder over hundreds of kilometers into the Gulf of Alaska (Fig. 11). In all of these cases, the fines are inherent to, and emanate from, the (subaerially exposed and desiccated) depositional systems, rather than from eolian saltation.

The critical attributes of the Devono-Mississippian mudrocks of Oklahoma and the greater midcontinent include the: 1) fine grain size restricted to very fine silt to very fine sand with little to no clay (mineral or grain size), 2) blanket-like deposition over large regions, including far from paleocoastlines, 3) deposition in both sediment-starved basins as organic-rich and chert-rich mudstone, and in carbonate ramps as dispersed fines, 4) ultimate sources in (primarily) the Appalachian orogen, and 5) enigmatic absence of proximal fluvio-deltaic feeder systems. We propose that delivery of the siliciclastic material to epeiric seas by eolian transport best explains this confluence of attributes.

The predominantly arid and seasonally dry climate in the mostly subtropical latitudes of Devono-Mississippian Laurentia would have favored dust emission. In this model, the Neoacadian Orogen clastic-wedge system (Catskill delta for the Devonian, and Price–Pocono–Borden delta and Mauch Chunk delta for the Mississippian) then provided the source of the material, which was derived from the orogenic highlands. Notably, the Catskill delta system drained (at times) regions thought to have been glaciated during the latest Devonian (Brezinski et al. 2008, 2010; Ettensohn et al. 2020a, 2020b), and thus would have delivered a ready and voluminous source of Al-poor rock flour to the alluvial plain, analogous to the Copper River Delta today (Fig. 11). Although no such glaciation is recorded outside of Gondwanaland for the Mississippian, the subaerial delta plains of the Price–Pocono–Borden–Grainger and Mauch Chunk deltas off the Appalachian Orogen would have supplied abundant silt, analogous to the Nile Delta today. Mass-balance arguments suggest that only 1% (363 km3) of the total quartz eroded from the Neoacadian orogen and delivered to the Mississippian fluvio-deltaic plains is needed to supply the volume of quartz in the Meramec and Sycamore systems (259 km3) (Supplemental File 5).

The eolian-delivery model holds that austral winter winds along the Appalachian orogen front and Zonda–Foehn–Chinook-like leeward winds translating into easterlies (Figs. 10, 12) deflated these subaerial delta plains and lofted siliciclastic fines for ultimate deposition into the Laurentian epeiric seas and emergent land regions to the (paleo) north and west (Devonian), and (paleo) west (Mississippian) (Figs. 10, 12). The strength of leeward winds can impose gusty conditions to generate intense dust storms, such as the Zonda winds off the Andes (Norte 2015; Puliafito et al. 2015; Isla et al. 2021), and Chinook winds that rework proglacial strata across the northern Great Plains today (Hugenholtz et al. 2007). This model also explains the dearth of clay minerals (and thus Al), given the lack of intense chemical weathering expected in this arid climate region and—for the latest Devonian—the source in a proglacial delta system.

This model predicts the occurrence of similarly uniform fines distributed in coeval units across a large region of the Laurentian epeiric sea, as indeed occurs (Figs. 2, 3). Additionally, it builds on several other examples of eolian–marine units posited for similarly enigmatic, fine-grained (and generally clay-poor) siliciclastic systems that lack fluvio-deltaic components, such as those of the Permian Delaware Mountain Group (Fischer and Sarnthein 1988), Phosphoria Formation (Carroll et al. 1998), various units of the Permo-Pennsylvanian (Soreghan 1992; Soreghan et al. 2008), and Lower Mississippian units of the Appalachian Basin (Cecil et al. 2018a). Similarly, it builds on the eolian-dust model for the origin of chert in Paleozoic units of Laurentia proposed by Lowe (1976) and Cecil et al. (2018b).

The Woodford Shale consistently exhibits the finest modes of the study units and the highest silica content (Figs. 7, 8). Its uniquely high silica content is commonly attributed to the influence of biogenic silica from nutrient-rich upwelling (Turner et al. 2016), but this explanation evades the root issue of an ultimate Si source. Cecil et al. (2018b) posited that much Phanerozoic chert was likely sourced by eolian dust; this model seems viable for the Woodford Shale given the paleogeographic context of arid conditions across the midcontinent of Laurentia during the Late Devonian, conducive to the availability and transport of ultrafine siliceous dust. We propose that the fine detrital modes reflect: 1) the distance to source and 2) the direction of surface wind circulation patterns, given the significant paleogeographic shifts from Late Devonian to Early–Middle Mississippian times (Figs. 10, 12). During the Late Devonian, Oklahoma and the greater midcontinent seaway were situated west-northwest of the Catskill delta at ∼ 30° S. Prevailing winds off the leeward front of the Appalachian orogen would have transported eolian dust in a north and west-northwesterly direction, across the North American midcontinent to Oklahoma and beyond. The emergent transcontinental arch may have provided temporary storage areas along this path.

Ultimately, the finest fraction accumulated in the Woodford Shale; by the Mississippian, shifts occurred in both the location of the primary source (now the Price–Pocono–Borden–Grainger, the Mauch Chunk deltas, and coastal clastic systems) and the paleolatitude (∼ 15° northward). By this time, the Price–Pocono–Borden–Grainger and Mauch Chunk deltas were east-southeast of Oklahoma, enabling prevailing easterlies to loft fines from these systems directly into the epeiric seas across the greater midcontinent. With the northward drift of the Laurentian continent, Southern Hemisphere winter winds off the Appalachian orogen shifted to blow more directly toward the west and west-northwest (less circuitous than those of the Late Devonian). This shift might explain the elevated grain sizes, higher volume of quartz, and more “loess-like” signature in the Meramecian samples of the Meramec and Sycamore (Fig. 8). Additionally, orographic effects of the greater Appalachian orogen presumably further exacerbated the aridity across Laurentia, leading to enhanced eolian emissions.

Finally, widespread delivery of only minimally (chemically) weathered eolian fines to Devono-Mississippian epeiric seas may have had auxiliary effects on organic carbon in these systems. The Mississippian (and Upper Devonian) of North America contain anomalously high organic carbon, with the Woodford Shale and correlatives considered world-class source rocks (Ulmishek and Klemme 1990; Kuuskraa 2011; Sonnenberg 2011). We speculate that this remarkable organic richness might reflect eolian delivery of abundant nutrients (notably iron) to these distal epeiric seas, stimulating productivity and associated preservation of organic carbon, as has been observed in the Bering Sea (Koffman et al. 2021) and suggested in other deep-time, dust-sourced systems (e.g., Carroll et al. 1998; Gabbott et al. 2010; Sur et al. 2015; Abadi et al. 2020).

  • Siliciclastic silt and mudstone are common throughout Upper Devonian–Lower Mississippian epeiric-sea strata of Laurentia, but no proximal fluvial–deltaic feeder systems exist to explain transport of this material.

  • Mineralogical data indicate a clear detrital source. Grain-size modes range from very fine silt to very fine sand, with minimal lateral variation in units across broad swaths of Oklahoma and the greater midcontinent region. Geochemical signatures are Si-rich and Al-poor, reflective a dearth of clay minerals (particularly Woodford Shale and Sycamore systems). Detrital-zircon geochronology indicates primary sources in the Appalachian orogen.

  • Paleogeographic reconstructions show that the midcontinent occupied the subtropical arid belt, in the orographic rain shadow of the growing Appalachian orogen. Consequent aridity would have promoted dust emission from vast delta plains draining the Appalachian highlands. Model simulations of surface wind circulations indicate that Southern Hemisphere winter months were conducive for eolian delivery of fine siliciclastic material toward the west and northwest, derived from Appalachian clastic wedges. Alpine glaciation in these highlands during at least the latest Devonian produced large volumes of physically weathered fines.

  • Siliciclastic silt of the Upper Devonian Woodford Shale most likely experienced a circuitous eolian trajectory and greater transport distances, resulting in the predominance of very fine silt modes. In contrast, Mississippian units were directly downwind of deltaic plains of the northern Appalachian system, resulting in generally coarser silt to very fine sand modes.

  • Coeval units located most distally contain abundant chert, consistent with sourcing from very fine silica dust. Eolian transport of material that was minimally chemically weathered to Devono-Mississippian epeiric seas likely delivered essential nutrients that stimulated organic productivity in these commonly organic-rich units—a hypothesis that merits further investigation.

Supplemental material is available from the SEPM Data Archive: https://www.sepm.org/supplemental-materials.

This paper grew out of a graduate student seminar in the School of Geosciences at the University of Oklahoma led by Drs. Lynn Soreghan and Andrew Cullen. We thank the University of Oklahoma School of Geosciences, the Eberly Family Chair (Soreghan), and Warwick Energy (Cullen) for funding, Vyetta Jordan and the Oklahoma Petroleum Information Center (OPIC) for access to core, Cody Totten for contributions to initial data collection, the University of Arizona LaserChron Center NSF grant (NSF-EAR 1649254) for detrital-zircon acquisition and Mehrdad Sardar Abadi for discussions. We thank the reviewers and editors for their constructive comments on earlier versions of this manuscript. We also thank Patrick Kelly and Oluwatobiloba Francis Oyebanji for assistance with field data acquisition.

Open Access CC-BY 4.0.