Mean annual temperature (MAT) on the USGP is typically 15° to 18°C and mean annual precipitation (MAP) ranges from <250 mm in the west to >1,500 mm in the east, spanning the transition between semi-arid to sub-humid climates. The latitudinal zonation of northern to southern ecoregions generally follows temperature trends (Figure 2A); the short-to-tall grass prairie transition parallels the strong west-to-east MAP gradient (Figures 2B,C). Grassland ANPP is highly correlated to precipitation, with peak values of >700 g m^-2 along the eastern margin of the Plains and decreasing to 80 g m^-2 in the fore of the Rocky Mountains (Sala et al., 1988). On the USGP, terrestrial carbon storage may range from ∼0.3 to 0.9 kg C m^-2 (Derner et al., 2006; Petrie et al., 2016), but can be a net atmosphere carbon source through diminished evapotranspiration and/or increased soil erosion (Meyers, 2001).

The USGP arid–humid divide roughly coincides with the 100th meridian, where precipitation decreases from ∼600 mm in the east to ∼400 mm in the west (Figure 2B) (Nielsen, 2018; Seager et al., 2018). This gradient shifted about 250 km eastward between 1931 and 1940 compared the prior MAP between 1895 and 1930, encompassing the largest extent of previously recognized severe wind erosion during the DBD (Figure 2D). Drought and land use directly impact grassland functioning via an increased risk of degradation, desertification, and subsequent reduction in stability and productivity, such as that observed during the DBD (Cook et al., 2009; Koerner and Collins, 2014; Ruppert et al., 2015; Lei et al., 2016; Byrne et al., 2017; Hu et al., 2018). However, the DBD was not a homogeneous event in time or space but consisted of several droughts and relative wet phases affecting different regions at different times (Figure 3) (Laird et al., 1998). Drought conditions were exacerbated by above-average temperatures exceeding 40°C and land–atmosphere interactions that increased the threshold for precipitable water (Cook et al., 2011a,b, 2014; Su et al., 2014; Donat et al., 2016). Climate modeling studies of the DBD have found that land cover changes increased the intensity and altered the spatial footprint of drought across the USGP (cf. Cook et al., 2008, 2009; Hu et al., 2018).

Reduced vegetation cover allows for enhanced wind erosion of soils, burying of adjacent grasses by the accumulation of eolian sediment and an increase in overland flow and riling from episodic and extreme rainfall events (Schlesinger et al., 1990). The subsequent recovery of plant communities may be delayed a decade or more, even with 5 years of precipitation above the historical average, as observed after the DBD (Weaver and Albertson, 1956, pp. 128–162). Ecological studies in the 20th century indicated that the diversity of grassland species decreased with drought, excessive grazing, and range fire (Albertson and Weaver, 1944; Tilman and Downing, 1994; Collins et al., 1998). The climax grassland species, such as Blue Grama (Bouteloua gracilis) and Buffalo Grass (Buchloe dactyloides), have decreased vitality with drought because of a shallow root system (<1 m). As drought persists for grasslands, the surficial heterogeneity increases with the dominance of bare surfaces with biomass and nutrients focused by deeply rooted (2–5 m) woody vegetation, forming “islands of fertility” (Weaver and Albertson, 1943, 1956, pp. 86–116; Schlesinger et al., 1990; Jurena and Archer, 2003).

In semi-arid regions such as the USGP, there are also apparent responses in eolian and fluvial geomorphic systems to extreme precipitation variability (Ewing et al., 2006; Derickson et al., 2008; Ewing and Kocurek, 2010; Ravi et al., 2010; Turnbull et al., 2010; Belnap et al., 2011; Liu and Coulthard, 2015). The pattern of stabilized dunes is often linked to bioclimatic variations. Wetter spring and summer conditions and heavy winter snowfall provide excess moisture for eolian systems to stabilize more effectively, with colonization by climax vegetation assemblages (Cordova et al., 2005; Ravi et al., 2010; Turnbull et al., 2010). Drier conditions are often correlated with decreased vegetation coverage and increased eolian activity, such as during the Medieval Climate Anomaly (MCA) when multiple, multi-decadal droughts propagated across the Americas with numerous presently stabilized eolian sand systems apparently reactivated, but to an unknown spatial extent (e.g., Mason et al., 2004; Forman et al., 2005, 2008; Lepper and Scott, 2005; Seifert et al., 2009; Halfen and Johnson, 2013). Many periods of eolian reactivation exhibit discontinuous blow-outs, parabolic dunes and sandsheet accretion, often leaving a spatially heterogeneous landform assemblage (Schlesinger et al., 1990; Hugenholtz and Wolfe, 2005, 2006). Questions remain if such degradation is in response to landscape-scale drier conditions, or if soil erosion is a stochastic process reflecting randomness or localized surficial disturbance of vegetation with drought, grazing, fire and/or pestilence (e.g., Schlesinger et al., 1990; Bel and Ashkenazy, 2014).

A multi-criterion, spatially explicit weighted overlay is employed to identify USGP counties that are representative of the broader landscape for the time and place of negative acquisition and to minimize subjectivity in site selection (Table 1). We include three climatic parameters to estimate the magnitude of anomalies in each county between 1931 and 1940, with larger deficits given greater weight: percent change in maximum MAT and MAP relative to the preceding 35-year period (1895–1930), and the central tendency of evapotranspiration as expressed by the mid-range Palmer Drought Severity Index (PDSI_MR; Figures 5A–C). Three physical surficial parameters averaged at the county level were included to evaluate the similarity of landscape setting amongst counties within the same ecoregion, with values closest to the regional mean weighted highest: sediment texture (% sand, silt, and clay in upper 20 cm of soil), density of fluvial drainages expressed as a ratio of total channel length to county area, and percent area under cultivation per the Agricultural Census of 1935, identified in historical analyses as the year of maximal cropland extent on the USGP (Cunfer, 2005; Figures 5D–F). These criteria are standardized to a scale of 1–9 and given equal influence in an overlay to assign each county a final rating that reflects the intensity of drought conditions and similitude of surficial properties to the broader ecoregion (Figure 5G). A fishnet of 1,000 km² cells is created across counties within USGP ecoregions, and a random sample of cells distributed across the precipitation and temperature gradients is selected for acquisition of aerial photographs.

Access to the primary photo archive is through the National Archives and Record Administration[NARA] (2017) II facility in College Park, MD, United States, where much of early aerial imagery has been centralized, with limited holdings at smaller archives (Rango et al., 2008; Sylvester and Rupley, 2012). At the facility, photograph reproductions are available as hardcopy photomosaics of the original flight lines, organized by record group, state, county, and year. Once the appropriate survey symbol, flight line and photo number associated with an area of interest are identified, the specific canister of original reel film, typically held at off-site cold storage, is ordered to NARA for digitization by a vendor with approved scanning equipment. The original scans used in this study are archived by county and available from the Historical Air Photo Dataverse at the Texas Data Repository (Bolles, 2018a–i). A subset of 657 high resolution (≥1,200 dpi) photogrammetric scans of overlapping negatives from six counties is utilized in this study (Table 2 and Figure 6). These images were taken at 1:20,000 scale between November 1936 and November 1940, with 35% lateral overlap between images and forward overlap from ∼35 to 67%.

In lieu of coeval camera calibration reports, which are largely unavailable, the interior orientation of digitized images is corrected to pseudo-calibrated fiducial coordinates from a template image frame via an affine transformation in MATLAB (Leedy, 1948; Redweik et al., 2009; Gonçalves, 2016; Nagarajan and Schenk, 2016; Giordano et al., 2017). Images are then cropped to standard dimensions to remove outlying pixels from the frame, such as those obscured by the photograph collar and identification numbers (Heipke, 1997; Redweik et al., 2009; Doneus et al., 2016). Ground controls points (GCPs) are identified to correct for exterior orientation between image and coordinate space. Absent detailed flight logs containing roll, pitch, yaw, and position coordinates, GCPs are obtained by matching time-invariant features (TIFs) within the image to features in contemporary imagery with a known location (Nagarajan and Schenk, 2016). The base map within ERSI ArcGIS and vector data of roads from state GIS data repositories are used to derive 5–20 GCPs per frame to attain an RMSE ≤ 5 m.

Color correction is concentrated on contrast enhancement to distinguish the intergrades in brightness from bare to fully vegetated surfaces, and between the textural properties of those surfaces (Kadmon and Harari-Kremer, 1999; Murray et al., 2013; Liu and Mason, 2016). The open-source, python-based GNU Image Manipulation Program (GIMP) is utilized to normalize the histogram of each frame individually, then to gamma-balance across frames so overlapping features are displayed along the same radiometric range (Kadmon and Harari-Kremer, 1999; Murray et al., 2013). Images are monochromatically balanced to correct for light fall-off introduced by lens distortion, spatially variable characteristics of the original film, and/or improper exposure of the photograph (Figures 4B,C; Redecker, 2008; Morgan et al., 2010). Frames are then further color balanced in an ESRI ArcGIS mosaic dataset using a first-order smoothing algorithm to minimize merging of pixels that do not match. Frames are blended together along seamlines generated from spectral patterns of overlapping features; the frames with greatest spatial accuracy (i.e., lowest RMSE) are weighted forward in seamline delineation and blending.

Multi-view stereoscopy of overlapping aerial photography is a rapidly developing method for landscape-level analyses of historical surfaces (Gomez, 2012; Nebiker et al., 2014; Doneus et al., 2016; Mertes et al., 2017; Mölg and Bolch, 2017; Pacina and Popelka, 2017; Sevara et al., 2018). The Agisoft Photoscan (APS) software package is used for automated three-dimensional reconstruction of the surface; the APS algorithm generates dense point clouds on par with airborne LiDAR data to derive a triangulated surface representing a DSM (Figure 4D) (Nebiker et al., 2014; Pacina and Popelka, 2017). Subsequently the DSM is used to calculate slope and elevation as covariate parameters for image classification. The centroid coordinates of each frame footprint created in the mosaic dataset is utilized to approximate the camera position at the time of acquisition to increase accuracy in photogrammetric alignment within APS.

The Iterative Self-Organizing Data Analysis (ISODATA) classification in ESRI ArcGIS is a robust clustering algorithm for areas with sparse ground-truth data and is used initially to understand the distribution of pixels DNs within a mosaic (Ball and Hall, 1967; Liu and Mason, 2016). An exaggerated number of clusters (≥15) is established across the range of spectral values, and neighboring clusters are merged until the overlap between cluster quartiles is ≤1% (Figure 4E) (Anderson and Cobb, 2004; Skirvin et al., 2004; Okeke and Karnieli, 2006; Bolles et al., 2017). A majority filter is applied to reclassify pixels where it can be reasonably assumed the pixel belongs to the surrounding cluster, producing a smoothing effect (Mather and Koch, 2011; Liu and Mason, 2016). Whereas this is a powerful tool to assess vegetation cover and soil texture, distinctions in land use and geomorphic form are often obscured. For example, an unpaved road, eroding field, and migrating dune are clustered together given similarities in surficial brightness. Thus, more information about patterns between adjacent pixels is needed to classify surficial processes and land use.

Using a MATLAB script, we calculate gray-level co-occurrence matrices for several textural properties within a 3 × 3 moving window (Haralick et al., 1973; Laliberte et al., 2004; Morgan and Gergel, 2010, 2013; Morgan et al., 2010; Mather and Koch, 2011; Schowengerdt, 2012; Vogels et al., 2017). This returns the probabilities of occurrence of a specific pixel pairing in eight directions, the average of which is taken to define the final texture parameter for each texel (i.e., the block of pixels within the moving window). The resulting texture rasters are used as inputs for principal component analysis (PCA) to identify statistically significant relationships between texture parameters and group pixels for object-based image segmentation (Figure 4F). Segmented objects are overlaid with the DSM to categorize landform and land use, then combined with surficial soil texture and tonal class to obtain the final classification of unique surface types (cf. Morgan and Gergel, 2010, 2013; Blaschke et al., 2014; Parajuli et al., 2014; Baddock et al., 2016; Parajuli and Zender, 2017; Vogels et al., 2017).

There is a wealth of observations on land surface conditions during the 1930s that has previously provided the basis for environmental interpretations (Browning et al., 2009; Skaggs et al., 2011; Williamson et al., 2011). Data sources are from contemporaneous peer-reviewed publications, government records from the Soil Conservation Service (SCS), ancillary experiment stations, and ground-based photographs held at NARA II facility in College Park, MD, United States. Items reviewed and digitized include: handwritten field notes, field photographs, planimetric maps, documentation of tillage operations, and correspondence between land owners, scientists and government officials. A complementary data source is the study of near-surface eolian deposits and landforms that are associated with landscape degradation during the DBD (e.g., Forman et al., 2008; Bolles et al., 2017). Specifically, the dry Munsell Color of bare sand and buried soils is a useful corollary to connect image grayscale and soil texture. This index designates the hue, chroma, and value of a given material, where value (or brightness) is the first parameter established and ranges in equal intervals from black to white. The Web Soil Survey from the National Resource Conservation Service (NRCS; Soil Survey Staff [SSS] et al., 2017) is also used to delineate the lateral variability in soil texture. Finally, an Abrams stereoscope is used with hardcopy image stereo-pairs to randomly “spot check” the sense of slope and elevation measurements taken from the DSMs.

Bivariate and multivariate statistical approaches are applied to parse the potential relationships between a matrix of 52 variables measured from the six study sites. Variables include criteria from the site selection scheme (Table 3), classified land surface parameters (Table 4), and summer PDSI from 1931 to 1940 interpolated from the North American Drought Atlas (Figure 3) (Cook et al., 2010). Data are standardized, and a comparison of skewness and standard error values indicates about half of variables measured are not normally distributed. Therefore, Spearman’s rank correlation coefficient (R_S) is calculated to explore the strength and direction of dependence between the rank order of variable pairs (Borradaile, 2003). Pairs with an R_S over ±0.8 and a p-value ≤0.05 are considered as significant variables with a strong monotonic relation. A PCA of the same matrix reveals linear combinations of variables with the maximum amount of covariance, and the subsequent principal components (PCs) are mutually uncorrelated (Wilks, 2011). A threshold is applied to extract coefficients in the upper or lower quintile of PC eigenvalues and identify the dominant variables contributing to intra- and inter-site variance. Finally, two pairwise distance metrics are computed between sites to create hierarchical cluster trees (Wilks, 2011): cosine (one minus the cosine of the angle between points) and hamming (percentage of coordinates that differ).

Ward County, ND (Figure 7A) occurs at the transition from northern mixed-grass prairie to steppe vegetation and by 1935 was ∼60% cultivated. This county ranked 7, experiencing maximum temperatures 11% higher in the 1930s than in the prior three decades, and with a 13% precipitation deficit. At the time of negative acquisition in the fall of 1938, most of the study area (∼70%) is an uncultivated mix of grasslands over alluvial floodplains and lacustrine basins, with 71% vegetation cover irrespective of land use or landform. Nearly 40% of uncultivated areas with plant cover are classified as patchy weeds and grasses, 25% as denser grasses and forbs, and ∼5% as woody vegetation. Cultivated fields appear mostly vegetated (73%) with severe cover loss when bare. Of the most extreme denudation, 80% is associated with natural landforms, at least half of which are identified as dried-out lake beds exposing silty clay loam. Nearly all lacustrine basins lack standing water; by 1939 total lake surface area was about one-third of the recorded high from 1868 to 1877 (Shapley et al., 2005). In total, 242 surface types are identified, with over half the study area described by six types (Table 5). The most common surface type is grassland dominated by weeds and perennial grasses colonizing loamy soils developed on glacial till.

McPherson County, NE (Figure 7B), ranked 5, falls within the central mixed-grass prairie ecoregion, was 23% cultivated by 1935, and during the DBD experienced a 7% increase in maximum temperatures and a 19% deficit in precipitation. The soils in this county, at the margin of the Nebraska Sand Hills, are characterized by ≥90% sand content (Figure 5D). The earliest available aerial photographic coverage for this county was July 1939, with 22% of the surface in cropland, ∼40% uncultivated undulating sandsheets and 35% interdunal grassland. Nearly three-quarters of the surveyed area (74%) is vegetated, though severe vegetation loss is frequently mapped in cultivated fields. Surfaces classified with the greatest degradation occur 64% on cultivated fields and 31% on blowout dunes, nascent parabolic dunes, and sandsheets. The crests of sandsheet ridges are most frequently disturbed. Blowouts are trough-shaped, sometimes reaching 0.5 km in length along the dominant axis, oriented NNW/SSE. In total, 105 surface types are identified, with >50% of the study area described by three classes (Table 5). This area is ∼23% vegetated sandsheet, and nearly equally common are patchy-to-dense weeds, grasses, and forbs in interdunal grasslands.

Hamilton County, KS (Figure 7C) yielded a ranking of 8, and experienced significant drought conditions during the DBD, with 5% higher maximum temperatures and an extreme precipitation deficit (-25%). Shortgrass prairie dominates this county, with largely sandy soils and smaller areas of clay loam soils. In 1935, agricultural census data indicates that 35% of the county was under cultivation, though only ∼5% is observed as actively tilled in the 1939 imagery used for this study, suggestive of farmland abandonment. Much of the vegetation present are grasses on alluvial surfaces and sparse cover on sandsheets, with the densest vegetation along fluvial channels. However, nearly half of the study area is devoid of vegetation at the time of negative acquisition. Almost 33% of bare areas are severely denuded with the presence of barchanoid-type ridges palimpsest with sand sheet deposits and occasional large blowouts (∼1 km). Slight-to-severe degradation is typically associated with sandsheets, nascent barchan dunes, and devegetated alluvial floodplains. Approximately 25% of cultivated surfaces are eroded, with sand apparently saltating across the landscape from the northern edges of fields. Cultivated fields and pastures exhibiting degradation are most often associated with loamy sand or clay loam soils, whereas nearly all bare uncultivated landforms are composed of sand or silty sand. In total, 241 surface types are identified, with over half the study area described by six classes (Table 5). Most common in this study area is vegetated sandsheet on sandy soil, followed by active sand dunes, and severely-to-extremely denuded sandsheet surfaces.

The southern shortgrass prairie of Quay County, NM (Figure 7D) experienced slightly warmer maximum temperatures than previous decades but a large precipitation deficit during the DBD and is ranked a 6 in the site selection scheme (Table 3). A maximum of 15% of the county was under cultivation in 1935, and few cultivated surfaces are observed in the 1936 photographs analyzed in this study. Forty-five percent of the study area is bare; where vegetation is present, it is typically grasses and shrubs. Extreme degradation is largely associated with sand dunes (>75%) which most frequently are parabolic or partially reactivated barchanoid ridges, oriented NE–SW. Severe degradation is also observed for some fluvial surfaces, such as point bars and tributaries of the Canadian River, and sheet and rill erosion above incised channels. In total, 232 surface types are identified, with over half the study area described by nine classes (Table 5). The largest proportion of the study area is mapped as vegetated sandsheets developed on antecedent eolian sands, followed by severely degraded sandsheets.

Tillman County, OK (Figure 7E) falls at the confluence of the three dominant ecoregions of the Great Plains: short-, mixed-, and tallgrass prairie. This county, rated 7, experienced an 2% increase in maximum temperatures and 11% precipitation deficit during the DBD, and in spring of 1937 exhibited ∼40% bare surfaces. Vegetated surfaces are evenly distributed across grasses, woody vegetation, and cropland. Cultivated surfaces are frequent; over half of fields are bare, with 25–30% classified as severely denuded. However, of the unvegetated areas ∼60% is associated with fluvial deposits in the main channel and tributaries of the Red River, with dunes superimposed on terraces and partially active, reflecting a complex interaction between eolian and fluvial geomorphology and potentially influencing the evolution of both features (Liu and Coulthard, 2017). The majority of moderate-to-severe degradation occurs on cultivated fields, with the remainder common to older, topographically higher fluvial terraces with active sandsheets. This county exhibits a marked increase in diversity of soil type, land use and landforms than those counties west of the 500 mm precipitation isohyet. In total, 452 surface types are identified, with over half the study area described by 14 classes (Table 5). The greatest portion of denuded area is mapped as bare fluvial deposits of loamy sand, followed by cultivated fields with sandy loam soils.

Lamar County, TX (Figure 7F) covers the transition between southern tallgrass prairie to crosstimbers ecoregion and experienced largely neutral conditions during the DBD, receiving a rating of four. Approximately 45% of the study area is classified as cultivated fields and ∼50% associated with the fluvial landforms and floodplain of the Red River. Less than 30% of the surveyed area is bare, and a robust population of woody vegetation colonizes uncultivated surfaces. Generally dense vegetation is absent where land was cleared for cultivation. The majority (50–60%) of bare surfaces occurs in cultivated fields, and many fields relict show nebkha dunes, ∼10 to 25 m in diameter and up to 2 m in relief (cf. Seifert et al., 2009). Among natural surfaces, deposits within the active channel of the Red River lack plant cover. Much of the moderate-to-severe degradation is associated with gullies on older fluvial terraces, whereas incipient vegetation loss is most frequently observed across the active floodplain. In total, 343 surface types are identified, with over half the study area described by sixteen classes (Table 5). The largest percentage of the study area is mapped as vegetated cultivated fields on clay soils, and dense vegetation on uncultivated loamy sand, silt loam, and clay soils.

A statistical analysis was undertaken to evaluate the relation amongst land cover classes, landscape attributes and climatic variables. The extent of extremely denuded surfaces is positively correlated to loamy sand soils and fluvial deposits (e.g., point bars, alluvial islands; Figure 8A). Areas of severely denuded surfaces are negatively monotonic with silty loam soils, cultivated fields, and MAP (Figure 8B). Patchy cover of weeds and grasses increases where the summer PDSI of 1932, 1937, and the mean PDSI from 1931 to 1940 are lower (Figure 8C). Dense, woody vegetation coverage increases with area of sandy loam soils, the density of fluvial drainages, the mid-range PDSI value during the DBD, and the number of surface types explaining at least 50% of site surficial diversity (Figure 8D). The association between woody vegetation and greater diversity of surface types is further reinforced by shared monotonic trends with the summer PDSI for 1934 and 1935 (R_S = 0.943, p = 0.017), the change in maximum MAT (R_S = -0.841, p = 0.044), the average percentage of clay in county soils (R_S = 0.886, p = 0.033), and dominant ecoregion (R_S = 0.882, p = 0.033). The total count of surface types classified is positively correlated with area of silty loam soils (R_S = 0.88, p = 0.05) and the proportion of area under cultivation by 1935 (R_S = 0.886, p = 0.033). There is a positive correlation between areas of silty loam soils and cultivated fields (R_S = 0.941, p = 0.017). Cultivated fields are negatively correlated with sandsheets (R_S = -0.841, p = 0.044), though sandsheets trend positively with active dunes (R_S = 0.941, p = 0.022).

Principal component analysis reduces the 52 variables included to three PCs that explain >80% of the variance, and >90% of the variance is accounted for by including the fourth PC (Figures 9A,B). In component space, Ward County, ND plots in the lower quintile of the first three PCs, most closely associated with silty loam and silty clay loam soils, grassland, lacustrine basins, and change in maximum MAT. McPherson County, NE ranks in the lower quintile of PC1 and the upper quintile of PC3, with weed and grass cover and the average percentage of sand in county soils. Hamilton County, KS falls in the lower quintile of PC1 and the upper quintile of PC2, most closely associated with areas of mapped sand, clay loam soils, sandsheets, and the percentage of area uncultivated by 1935. Quay County, NM is in the upper quintile of PC2, generally correlated with areas of exposed bedrock, dunes, and moderately-to-extremely denuded surfaces. Tillman County, OK is in the upper quintile of PC1, most closely associated with: the total count of surface types classified, summer PDSI of 1934, 1935, and mid-range during the DBD, drainage density, and occurrence of sandy loam soils and woody vegetation. Lamar County, TX plots in the upper quintiles of PC1 and PC3, and in addition to the variables associated with Tillman County, is also correlated with 50% and 90% diversity thresholds in surface type, clay soils, percentage of land cultivated by 1935, and MAP.

As the PCA indicates four components explain >90% of the variance between sites, four hierarchical clusters elucidate the covariance in surficial response across sites. Cosine distance of site linkages provides a relatively reliable cluster tree, with a cophenetic correlation coefficient (c) of 0.676 (Figure 9C), though hamming distance determines a somewhat stronger correlation (c = 0.748, Figure 9D). The former measure demarcates four clusters as: (1) Tillman County, OK and Lamar County, TX, (2) Quay County, NM, (3) Ward County, ND, and (4) McPherson County, NE and Hamilton County, KS. For both distance measures, the largest break occurs between cluster 1 and clusters 2–4; however, cosine distance closely links Tillman and Lamar counties where hamming distance creates a discrete separation. Likewise, cosine distance calculates close linkages amongst clusters 2, 3, and 4 where hamming distance estimates more distinct, gradational steps between sites.

The DBD is a vivid example for an annual record of drought on the USGP, though is modest in intensity and duration compared to droughts in the past 1,000 years (Cook et al., 2007, 2010; Woodhouse et al., 2010) and forecasted droughts for the 21st and 22nd centuries (Cook et al., 2015). Spatial analysis of countywide anomalies in maximum MAT, MAP, and PDSI from 1931 to 1940, concomitant with average soil texture, density of linear water features, and extent cultivated by 1935, illuminates the magnitude of DBD aridity in USGP counties and surficial coherence within ecoregions. The distribution of the highest ranked counties (8) reveals a cluster of impacted counties on the northern USGP in the Dakotas (Figure 5G). The greater clastic influx between 1930 and 1950 observed in lake cores from the Dakotas and Minnesota underscores this regional aridity (Dean, 1997; Shapley et al., 2005; Wigdahl et al., 2014). Dust storms, like that observed November 12, 1933, moved progressively over Bismarck, ND and Omaha, NE, with suspended particles eventually raining out across southern and eastern states (Hovde, 1934; Miller, 1934). A less continuous cluster of southern counties occurs within the classically defined DBD area of severe wind erosion, constrained by a maximum MAT range of 20 to 25°C, MAP of >500 mm, and the ∼250 km eastward shift of this isohyet compared to the previous 35-year period (Figure 2). Interestingly, the highest ranked counties within this region are also where private hospitals were at capacity and Red Cross emergency hospitals were set-up to deal with a conspicuous increase in “dust pneumonia” (Chicago Tribune, 1935).

The apparent exception to this ranking system are the counties comprising the Nebraska Sand Hills, which ranked 4 and 5 despite persistent drought conditions. This suggests grassland response observed in this area may not be correlatable to other counties with respect to the surficial criteria included in the site selection overlay. The high infiltration capacity of sandy soils in tandem with the high levels of regional groundwater tables contribute to the dominance of baseflow control on local response, and baseflow is generally resistant to interannual climate variability (Wang et al., 2009). On annual time-scales the Sand Hills generate higher rates of mean runoff and recharge, correlated with lower actual evapotranspiration, than in adjacent areas with siltier soils and relative to the amount of precipitation received (Gosselin et al., 2006; Harvey et al., 2007; Wang et al., 2009). The consistency of groundwater levels since AD 1700 appears to buffer the impact of the DBD in this region of Nebraska. However, multi-decadal droughts during the MCA (∼AD 800–1200) were sufficiently severe for a regional decline in the water table, with concomitant landscape denudation and increase of eolian transport and deposition (e.g., Harvey et al., 2007; Miao et al., 2007; Halfen and Johnson, 2013).

The spread of landscape degradation identified reflects ephemeral and cumulative effects of multiple controls that vary with antecedent landscape conditions and persistence of climatic anomalies (Turner, 1989; Schlesinger et al., 1990; Peters et al., 2007; He, 2014; Petrie et al., 2016). PCA indicates that site variance largely parallels location east or west of the 100th meridian, coupled with the severity of the precipitation deficit during the 1930s. The secondary component is related to the latitudinal temperature gradient and predominant land use and geomorphology, with a tertiary influence from soil texture and/or the annual PDSI during the DBD. Similar studies of USGP grasslands find that many ecosystem processes are primarily driven by precipitation and modified by temperature (Petrie et al., 2016; Nielsen, 2018; Seager et al., 2018). This analysis is also consistent with a linear relationship between drought intensity and ecosystem resistance observed in a global analysis of grassland response to ∼500 drought events (Ruppert et al., 2015).

Clustering analysis illustrates the differences in landscape stability across the arid-to-humid divide. Biotic controls tend to drive ecosystem stability east of the 100th meridian, with increased plant functional diversity and greater below-ground biomass of woody species (He, 2014; Lei et al., 2016). Both Tillman County, OK and Lamar County, TX, are characterized by mostly continuous dense or woody vegetation increasing with areas of greater soil clay content, near-river landforms, and reduced temperature anomalies, likely reflecting water availability. Tallgrass prairie is generally resistant to inter-annual precipitation variability and resilient to daily or weekly extremes (Jones et al., 2016), and more sensitive to monthly and seasonal precipitation distribution (Petrie et al., 2016). The distribution of the densest vegetation cover in this region is most strongly correlated with the summer PDSI of 1934 and 1935. This association would suggest ecosystem degradation in the eastern USGP was largely driven by plant response to the magnitude shift in seasonal conditions from 1934, recognized as the most severe drought year of the last millennium (Cook et al., 2014), to 1935, a wet period for the southeastern USGP (Figure 3E). Noteworthy in Lamar County, TX are the relict nebkha dunes revealed at the surface in areas cleared for cultivation. These features are similar to sites dated to 700–1,000 year BP on the Ozark Plateau and underscore the spatial extent and intensity of MCA droughts in the late Holocene (e.g., Seifert et al., 2009).

Abiotic controls appear to more strongly influence grassland degradation west of the 100th meridian and are associated with latitudinal temperature effects (Nielsen, 2018). The northern ecoregion of Ward County, ND experienced some of the most pronounced increases in maximum temperature during the DBD but is largely detached from other western sites in component space related to the overall lower MAT and dominance of siltier soils. Quay County, NM is likewise clustered individually but is near in component space to Tillman County, OK, as both are southern ecoregions with milder DBD conditions, and shows a greater diversity of surface types than other western counties. Nonetheless, this area of eastern New Mexico is statistically linked to Hamilton County, KS and McPherson County, NE, and correlated with the distribution of increasingly sandy soils and large patches of eolian landforms indicative of soil erosion. These latter two counties are both central ecoregions with a similar distribution of land use, high soil sand content, and exhibit less surficial diversity than eastern sites as sandsheets and dunes cascade across the landscape. McPherson County is influenced by groundwater levels at the margins of the Sand Hills and classified as mostly vegetated sandsheet and interdunal grassland with little change to extent of cultivation. Whereas Hamilton County shows ∼50% vegetation loss of shortgrass prairie, associated with active eolian landforms and a sharp increase in apparently fallow fields.

Broadly, the distribution of denuded surfaces during the DBD exhibits stronger non-parametric correlations with surficial properties than climatic factors. Land use interacts with degradation regardless of setting but is more significant to explain site variance in eastern counties, whereas uncultivated landforms dominate the drier, western areas. Sand-dominated soil textures on landscapes with maximum relief are correlated with more frequent denudation for western counties. The areal extent of blow-outs, parabolic dunes, sinuous ridges, and sandsheets positively covary, and this heterogenous landform assemblage correlates to intensity and persistence of 1930s aridity (Hugenholtz and Wolfe, 2005, 2006). Particle sources for the DBD have been attributed solely to areas of agricultural disturbance in the south-central Great Plains (cf. Bennett and Fowler, 1936; Johnson, 1947; Worster, 1982; Hansen and Libecap, 2004; Schubert et al., 2004; Cook et al., 2008, 2009; Lee and Gill, 2015). This analysis identifies large uncultivated areas that were equally or more a source of dust during the DBD than plowed agricultural fields. On the Northern Great Plains (North Dakota), one significant source of dust is with the drop of lake level and the exposure of lake sediments and adjacent eolian deposits. In the classic DBD area, the many stabilized, high-relief sand-rich dune systems adjacent to major rivers were partially to wholly denuded and were significant sources of dust.

A similar shift in the arid-to-humid transition as observed during the 1930s is projected for the future, which can strongly influence settlement sustainability and agricultural development (Nielsen, 2018; Seager et al., 2018). The overuse of groundwater coupled with rising costs of irrigation could cause groundwater tables to lower and increase susceptibility of previously stable systems (Parton et al., 2007). Future warming may reduce vegetation cover and soil moisture over large areas of USGP grasslands, exacerbating soil erosion, degradation feedbacks, and decline of terrestrial carbon sinks (Schlesinger et al., 1990; Lei et al., 2016; Hu et al., 2018). Regional hydroclimatic variability and carbon sequestration is therefore dependent on how land surface changes modulate climatic conditions (Cook et al., 2013; Petrie et al., 2016; Hu et al., 2018). Given the potential speed and irreversibility of such changes (Breshears and Barnes, 1999; Collins et al., 2014; Moran et al., 2014; Svejcar et al., 2015), characterizing USGP grassland response during severe drought is significant to constrain the large variability in landscape degradation (Petrie et al., 2016), and to better guide adaptable land and carbon management (Basara et al., 2013; Ruppert et al., 2015; Lei et al., 2016; Petrie et al., 2016; Byrne et al., 2017).