Introduction

Geohazards pose a significant threat to people and property and can cause significant damage to infrastructure and economic loss, often resulting in governments and non-governmental organizations expending significant resources to mitigate their impacts (Rogers 2004; Mahalingam et al. 2016). Hence, landslide susceptibility mapping is crucial for estimating loss potential for earth movement (Paudel et al. 2003). Rapid urbanization often erases the physical evidence of past landslides, which can only be evaluated by examining older aerial images or reports that notice changes in the geomorphology and environmental responses to anthropogenic alterations. Mass wasting features like shallow earthflows, rotational slumps, translational slides, sand boils, or lateral spreads are common features along erosional escarpments in the New Madrid Seismic Zone (NMSZ) (Lanbo and Yong 2001; Rogers 2004).

Missouri experienced an intense intraplate earthquake of Mw 7.2–8.2 on December 16, 1811, followed by aftershocks of Mw 7.4 that same day that caused multiple surface ruptures with mass wasting. Two other earthquakes occurred soon afterward: Mw 7.3 + − 0652.iu on January 23, 1812, and Mw 7.5 on February 7, 1812, along the Reelfoot fault in Missouri and Tennessee. Other moderate scenarios also correlated with substantial landslides and reports of sand blows following aftershocks of Mw 4.0–5.0, and about 2000 felt tremors between December 1811 and April 1812 (Reiss et al. 2005). The seismic energy emanating from the failed Precambrian rift triggered the 1811–1812 New Madrid earthquake sequence in Missouri at depths of 20 km to 30 km, accompanied by large-scale sliding of the Chickasaw Bluffs bordering the Upper Mississippi Embayment (Fuller 1912; Jibson 1985).

The potential for slope failure increases with ground shaking of saturated soils (Cruden and Varnes 1996). By 1985, Missouri had documented over 250 landslides in the NMSZ (Jibson and Keefer 1989) triggered by saturation from perched water tables and intense or semi-continuous precipitation. The United States Geological Survey noted a 35% probability that an Mw 6.0 quake would occur in the next fifty years (Petersen et al. 2014). We can reasonably expect that pre-existing landslides would be more susceptible to reactivating during earthquakes of Mw 6.0 or higher because of the number of equivalent load cycles. The safety factor would also be influenced by antecedent soil moisture from precipitation and natural or anthropogenic drainage, since soil moisture tends to be trapped within dormant landslides. This results in low permeability along their basal slip surface(s). For these reasons, the seismicity of the NMSZ would be expected to increase the risk of landslide recurrence.

Sustained ground shaking is another common trigger mechanism for slump-block slides, liquefaction-induced flow slides, and lateral spreads. Lateral spreads are commonplace during and after earthquakes (Honegger and Wijewickreme 2013). As ground shaking continues, the lateral flow movement of saturated cohesionless soils has finite expressions on gentle low-angle slopes induced by earthquakes (Varnes 1978; Rauch 1997; Kramer 2013). The movement is generally triggered by liquefaction of a saturated cohesionless horizon, usually composed of silt, sand, or occasionally gravel (Hansen 1966; Varnes 1978). Hansen (1966) identified slip surfaces within sensitive silty clay surrounded by stiff clay.

The dolomite–limestone formations in southeastern Missouri (Fig. 1) pre-dispose the upland slopes to landslides and slumps that tend to favor steep slopes underlain by unconsolidated material or thick soils. Karst features form in carbonate units, such as limestone and dolomite, and other soluble rocks like anhydrite, gypsum, and evaporated salts (Bansa 2018). Karst features are usually formed over long periods of time through cation exchange with oxygenated groundwater percolating along joints, fissures, and stratigraphic horizons. Mining of clay and crushed limestone requires frequent blasting and vibration of heavy equipment. Frequent vibration combined with small amounts of moisture (e.g. 10% by weight) can actually enhance the bulk density of cohesionless gravel-sand mixtures. The high solubility of karst terrain makes it of immense concern to authorities and residents.

Fig. 1
figure 1

Source: Missouri Department of Natural Resources)

The political map of conterminous United States (top left), the State of Missouri (top right), and the geologic map of Cape Girardeau and Bollinger counties and parts of Stoddard and Scott counties in southeast Missouri showing, the regional controlling faults trending northwest–southeast in the study area in the New Madrid Seismic Zone (

An even more significant problem might be the reactivation of dormant landslides or lateral spread features triggered by earthquake shaking in 1811–1812 in the NMSZ; approximately 100 earthquakes occur annually in the Wabash Valley Fault Zone (WVFZ), 220 km to the northeast (McBride et al. 2002; Frankel et al. 2009). An extensive lateral spread feature occupying 67 km2 has been identified near Advance, Missouri, placing it among the most extensive of such features globally (Watkins and Rogers 2009a, b). This lateral spread feature formed as part of the series of structural alterations that occurred during the 1811- 1812 New Madrid earthquakes. Much of the NMSZ lies beneath the Mississippi Embayment, a trough filled with sediments that may greatly amplify earthquake shaking (USGS 2002) (Fig. 2). Some sliding occurred in the lower Paleozoic interbedded shale and limestone in the dissected area near the Mississippi River (Radbruch-Hall et al. 1982; Thompson 1991).

Fig. 2
figure 2

Modified USGS Earthquake Shaking Intensity map for the NMSZ showing an insert of the study area. The 2017 intensity shaking map projects that areas along the Mississippi River are most at risk of earthquake-induced ground shaking (Rogers 2010)

Disasters can occur due to unmonitored hazards, highlighting the need to monitor threats to lives and properties (Davies 2015). Aerial photos of remotely sensed Landsat imagery and Light Detection and Ranging (LiDAR) datasets are commonly employed in a geospatial environment to evaluate anomalies in slope morphology commonly associated with mass wasting. These datasets have often been used to compare pre- and post-event imageries (Liang 1952; Balazy et al. 2019). Detailed landslide susceptibility mapping has historically focused on identifying various types of landslide features, whether active, dormant, or old (Fig. 3) (Fuller 1912; McGill 1959). Their inventory usually denotes at-risk areas by the density and scale of past land slippage. Workers in the field have long assumed that landslides seldom occur as isolated phenomena. They almost always occur in irregular coalescing complexes that result in infrastructure, human life, and socio-economic losses. LiDAR technology has proven effective in investigating landslides by creating accurate and precise representations of the topography (Jaboyedoff et al. 2010; Rogers and Chung 2016a, b). Geospatial technology has supplanted previous technologies as the preferred medium for landslide susceptibility mapping. Most of these Geographic Information Systems (GIS)-driven efforts have relied on the application of quantitative or "expert" techniques (Carrar et al. 1995). Studies elsewhere have compared the coefficient of global logistic regression and geographically weighted logistic regression models (GWLR), which noted that the GWLR models are more easily equipped to classify mass movements of varying scale (Zhang et al. 2016). Slope stability algorithms have also employed the simplified infinite slope equation to evaluate surficial slope stability relative to safety factors based on the presumption of full saturation (Palenzuela et al. 2015). In comparing the GWLR models' results with traditional susceptibility maps, GWLR models tend to predict more slope failures.

Fig. 3
figure 3

Map of Earthquake Features of the New Madrid District (1912) in Parts of Missouri, Arkansas, Illinois, Kentucky and Tennessee (Fuller 1912). The subset is a detailed map of the study area showing the Paleozoic uplands and paleo-landslides of Cenozoic deposits

This study employed LiDAR and GIS for landslide susceptibility mapping in Cape Girardeau and Bollinger counties and portions of Scott and Stoddard counties in southeastern Missouri. The main aim was to analyze the relationship between landslides and susceptibility to other geohazards by evaluating surficial morphology using statistical computation in spatial analyses, focusing on the impacts of geomorphic alteration by evaluating the observed changes (e.g. before vs. after). The overarching objectives are to ascertain the susceptibility of a study area to reactivation or generation of new mass wasting activities and identify the hazard potential as a function of land use practices for regulatory agencies.

Study area

The geo-political extent of the study area covers Cape Girardeau and Bollinger counties and parts of Stoddard and Scott counties in southeastern Missouri (Fig. 1). Greater Cape Girardeau is a key commercial center of the New Madrid Seismic Zone, which includes portions of Missouri, Illinois, Arkansas, and Kentucky. The topology of the area has changed significantly in the last 50 years, especially in pockets of expanding urbanization. The Mississippi Embayment (Fig. 1, in yellow) borders the study area in the southwestern section, while gentle hills covered by Holocene-Quaternary alluvium and wind-blown loess typify the northwestern section. The bedrock geology of the Cape Girardeau and Bollinger counties reveals northwest-southeast trending limestone, dolomite, and highly fractured sandstone. It is dominated in the north by rocks in the Ordovician Series, namely the St. Peter Sandstone and Everton Formation, Joachim Dolomite, Dutchtown Formation, and Decorah and Plattin Groups (Satterfield 1993). The significant faults and rivers in the Paleozoic units trend northwest-southeast, while the NMSZ faults trend southwest to northeast. The average annual rainfall in the study area is about 63–122 cm, with snow precipitation between zero and 25 cm. The precipitation drains into the Mississippi and St Francis Rivers. Floods periodically inundate the Mississippi River floodplain south of Cape Girardeau and Bollinger counties, flowing through the Bell City Gap and Upper Whitewater Creek. The northern flanks of Scott and Stoddard counties are covered by Quaternary alluvium of the Mississippi Embayment. The major faults trend northeast-southwest (Gomberg and Ellis 1994).

Methods

Hazard mapping was conducted in a phased approach: (1) desktop studies of historic resources supplied by the Missouri Geological Survey; (2) site reconnaissance and fieldwork; and (3) data validation and analyses. These steps enabled the construction of a landslide hazard map in 1:24,000 scale (1 in. = 2000 ft).

Desktop studies

The survey entailed the collection of historical maps dating back to 1875. Potential landslide areas were identified on historical topographic maps dating back to 1912 and 1963, followed by scanning and georeferencing images, converting all the scanned data into raster files, importing, and overlaying the digital information layers, and interpretation. The output generated maps of the study area at scales of 1:600,000 to 1:24,000. The data was compared to the 2016 USGS Landsat orthoimagery map to screen for landscape alterations. Manually identified landslides were digitized in the ArcGIS environment and overlain on the topographic map, along with lithology, faults, earth dams, and depth to bedrock maps (Fig. 4). The watershed layers include locations of ponds and watersheds, and direction of flows were annotated and imported into ArcGIS. Scanned files were georeferenced to the Missouri base map on a projection of North American Datum (NAD 83). Clusters of landslides were drawn by combining the external boundaries of the adjoining landslides under a single feature identity (ID). Line point features were created for the landslides, which allowed the individual landslide morphology to be drawn as a single identity (Fig. 4). Clusters of landslides were subsequently described by combining the external boundaries of the adjoining landslides under a single ID raster image generated from the LiDAR cloud points to provide a high-resolution hillshade raster images. Locations of dams, bridges, open-pit mines, and gravel quarries were also added to the base map.

Fig. 4
figure 4

The geological map of the Cape Girardeau and Gordonville Quadrangles showing mapped landslides (in gray) and the expected distribution of ground accelerations emanating from the active fault strands (in red)

Site reconnaissance and fieldwork

Site visits were carried out in December 2017 and May 2018 to determine the compatibility of the mapped landslides to current terrain, outcrop observations, and channel changes. Field work was used to evaluate various types of landslide clusters, as well as active sinkholes and flood-prone terrain. The mapped landslide clusters were surveyed.

Data validation and analyses

Several analytical tools were used to interpret the data. Landsat Digital Elevation Models (DEMs) were downloaded from the USGS website, and an orthophoto from the Google website. The DEMs were combined with LiDAR data downloaded from the Missouri Geospatial Information System (MOGIS) and processed to produce key analyses including slope, hillshade, GWLR and co-seismic hazard evaluation.

  1. a.

    Hillshade analysis The hillshade map provides the input for delineating the DEM, azimuth altitude, and scaling for the z-factor, which serve as a distinct expression of landslide morphology. The visualization of DEM allows the interpolation of the z-factor for image analysis. The hillshade extraction exposes the unique expressions of landslides as they exist. The image analysis simplifies the identification of anomalous topographic patterns that are commonly related to landslides and other mass movements. Hillshade measures the values for a topographic surface described as the illumination angle, and shadows in each raster, with the raster output values ranging from 0 to 255. Raster cells that are in shadows are assigned a value of zero, while fully illuminated raster cells are assigned a value of 255. The x and y units are expressed in decimal degrees, while the z units are in meters.

  2. b.

    Slope analysis The Landsat Digital Elevation Models (DEMs) downloaded for the study were combined and processed to enable analysis of natural slopes. Ashraf et al. (2012) described the slope function as a measure of hydrological processes such as overland flow, soil erosion, and sediment transport along slope fall lines (90° to slope contours). Slope was expressed as a function of the maximum change in height “z-value” in degrees (from 0 to 90) over a defined space or a raster pixel. The ArcGIS Slope tool calculates the maximum rate of change between each cell and its eight neighbors (Khadri et al. 2015). Consequently, every cell in the output raster has a slope value. The lower the slope value, the flatter the terrain, and; the higher the slope value, the steeper the terrain (the slope tool also provides critical corrections to slope calculations when the surface z units are expressed in different units from the ground x, y units). The slope can also be measured as a percentage value:

    $${\text{Percent}}\,{\text{of}}\,{\text{slope}} = \left( {{\text{rise}}/{\text{run}}} \right) \times {1}00\% .$$
    (1)

Given that the rise can be any value range from zero to near infinity, the slope value of a flat surface is 0%. When the slope angle equals 45°, the rise is equal to the run, and the slope value is 100%. Consequently, the slope angle measured in degrees (θ) is:

$${\text{Tan}}\left( \theta \right) = {\text{rise}}/{\text{run}}{.}$$
(2)

As the slope angle approaches vertical (90°), the slope percentage approaches infinity (Cao 2016). For surfaces tending towards vertical (higher than 45°), the percent rise becomes increasingly more significant, which makes the slope unstable for weathered rocks and cohesionless soils. The optical values were descriptive of the steepness of the slopes between 0° and 62°.

  1. iii.

    Geographically weighted logistic regression (GWLR) The GWLR allows the assignment of high statistical value to high-risk factors and low statistical value to factors considered to exhibit low risk. Logistic regression analysis is based on the premise that certain geologic factors must exist for failure to occur. The following landslide conditioning variables were used for analyzing the study area: precipitation and hydrology, lithology, gravity control, structurally controlled faults, slope steepness control, dam failure potential, sediment thickness, liquefaction, and earthquake shaking potential. The dependent variables like faults, lithology, and slope angles that might be expected to precondition a hillside for slope failure are normally assigned higher weight on a scale of 1–10 (where 1 is the minimum, and 10 is the maximum). Other independent variables, such as precipitation or dam failure, are assigned less weight if their impact or the likelihood of occurrence is low.

  2. iv.

    Co-seismic hazard evaluation For hazard evaluations, credible sources of ground shaking in the area were used. This process usually begins with a simple evaluation of the historicity of earthquakes in an area (Figs. 2, 3) using the map of Earthquake Features of the New Madrid District (1912), the New Madrid Seismic Zone map and the Earthquake Hazards Map of Southeast Missouri (Marcus 1993).

Results

Field observations

The site observations in some areas indicate steep slopes that experience repetitive mass wasting. Changes in vegetation were observed from the changes in soil characteristics and moisture resulting from shallow mass wasting on slope sides (Fig. 5). Figures 5a, b illustrate plants with yellow flowers and grasses with a deeper green color that is distinctly different from the surrounding grasses, respectively. In Figs. 5c, d, hillsides were observed to exhibit vegetated patches indicative of different adiabatic conditions. The verified landslides are characterized by shallow earth flows, rotational slumps, and block slides, while other areas were altered by building construction and infrastructure. Mass earthwork grading for site development was observed in several rural locations. Many hills mapped as landslide terrains have recently been excavated, re-worked, and re-compacted for residential developments with paved streets, lot pads for houses, man-made lakes and/or ponds, and golf courses within gated communities. Although some of the grading appeared to be in areas prone to earth movement, we were unable to discern any evidence of the nominal age of such features. Some of the hillside developments exhibited pavement cracks (Figs. 6a, b), which may be ascribable to differential settlement of compacted subgrade soils. Mass wasting was also observed along reaches of improved water courses (Figs. 6c, d). The subsidence of topsoil occasionally exposed underlying limestone where active sinkholes were impacting channels and public streets (Fig. 6e). Desiccation cracks prevalent in the southeast portion of the study area may be influenced by seasons of intermittent flooding and drying (Figs. 6f). The field observations were then integrated with the results obtained from the earthquake, precipitation, and seismic data map layers (landslide inventory). In Figs. 6c, d, the sides of the western bank of the Mississippi River exhibit the most severe angular differences in elevations.

Fig. 5
figure 5

Field pictures taken in May 2018 showing changes in vegetation. a Plants with yellow flowers due to the change in soil characteristics and moisture resulting from shallow mass wasting on slope sides. b Grasses with a deeper green coloration that is distinctly different from the surrounding grasses. c, d Hillsides with vegetated patches indicative of different adiabatic conditions

Fig. 6
figure 6

Field photos showing mass wasting in limestone terrain. a, b Cracks on the pavement as a result of differential settlement of the compacted soils. c, d  Channel bank undercut by active sinkholes. e Sink hole shows subsidence in the topsoil, which is indicative of an incompetent subsurface layer. f Mud cracks along the Buzzi Unicem rail tracks showing evidence of flooding of the tracks

Data interpretation and analyses

Hillshade analysis provided information on the bare-earth topology of the study area. The relief expressed in the hillshade shows the shape of hills and mountains illuminated by the sun’s relative position (levels of gray) on the map. The light gray shades indicate relative slope angles, mountain heights and peaks (Fig. 7). Dark gray shades indicate areas where there are relative steep angle changes that are unstable with high potential for slope failures, and landslides (Fig. 7b). In the lighter areas, changes are reflected by the difference in the illumination of a surface. Consequently, the variation in the intensity of light can expose the patterns and intricate changes in morphology over a larger area. The Advance, MO lateral spread feature is expressed in Fig. 7c as a flat, pancake-like feature straddled by two parallel faults. The intricate changes allow inferences based on geomorphic expressions typical of lateral spreads, and striations consistent with parallel faulting signatures (Figs. 7c, d).

Fig. 7
figure 7

Hillshade maps of the study area and snap shots of the areas reflecting the different types of landslides produced from LiDAR-derived DEM maps. The red lines are mapped faults, while the light blue lines are the mapped landslide and lateral spread features

The slope analyses describe the steepness and azimuth of the slopes (Fig. 8). The lower slope values (green to yellow) indicate relatively gentle undulating terrains, while the higher slope values (orange to red) indicate high angle slopes with steeply dipping terrains. The steepest areas in red have high landslide risks and are coincidental with riverbanks and watersheds. The lowest areas (green) are flat-lying alluvial plains currently cultivated for agriculture.

Fig. 8
figure 8

The slope analysis for the study area. The steepest areas (red) are the areas with high landslide risks, which are coincidental with river banks and watersheds. The lowest areas (green) are flat-lying alluvial plains currently cultivated for agriculture. The slope percentage reflects the steepness of the slope value

Furthermore, the GWLR analysis enabled a hierarchical evaluation of the most influential variables on slope failures. Gravity influences the migration of cohesionless earth materials pulled down slope by the force of gravity and is dependent on elevation, the severity of slope, and moisture cycles (wetting and drying). The influence of gravity is noticeable in areas of significant differential elevation and steep slopes. Slope steepness (Fig. 8) tends to control the topography, and the slope value of the study area reveals that hills that are sloped between 15° and 40° (from horizontal) are more at-risk than slopes inclined at lower angles. For this reason, the slope steepness factor was weighted as high risk. The thickness of overburden sediment averages 15–30 m. Surficial site reconnaissance suggests that the alluvial cover on the southwest flank of the Gordonville Quadrangle is > 7.6 m thick over Paleozoic age bedrock. However, the least sediment thickness occurs in dolomite, but the porosity of the crystalline bedrock increases because of systematic regional joints, which tend to be mutually orthogonal to bedding. It appears that the northwest-southeast trending faults in the study area exert a strong influence on the flow paths and geometry of the local watersheds and are considered high risk because of downslope flow.

Google orthoimagery was used to evaluate the current extent of mining activities in the southeastern portion of Cape Girardeau County (Figs. 9a, b). The hydrology of the Cape La Croix Creek has allowed local runoff to recharge the river channel as it flows into the Mississippi River channel (Fig. 9a). During the spring, the area receives more precipitation and runoff. Whenever precipitation continues unabated for days at a time, runoff is absorbed into joint-controlled “cutters” or soil-filled voids in the formation, which are observed to extend up to > 60 m deep. These provide natural conduits for selective point recharge, which led to a series of sinkholes around the northeastern corner of Buzzi Unicem’s open pit limestone quarry in 2007–2008, which extended to depths of > 91 m (Fig. 8b). The presence of linear soil-filled cutters can increase the risk of sinkholes developing if a low head outlet is established with an underground opening or with an open pit mine.

Fig. 9
figure 9

a The Google map imagery of the Cape La Croix River by Buzzi Unicem accessed October 16, 2020. b Field photo of the northeastern corner of the Buzzi Unicem mine pit showing the natural conduits for selective point recharge leading to a series of sinkholes

The co-seismic hazard potential susceptibility map combines all the risks various geohazards present in the study area and prioritizes the highest and/or combined threats. The co-seismic map derived from the analysis of GWLR risks, and Earthquake Hazards Map of Southeastern Missouri (Marcus 1993) (Fig. 10) evaluated the consequences of failure on the communities that may be affected. The gray zone in the earthquake hazard map highlights the erosional escarpments and slopes that are most susceptible to seismically induced landslides. Dam failure potential is considered a risk factor to the lives and infrastructure that are situated down-gradient. Seventeen dams regulated by the State of Missouri are located in the study area [regulated dams have a total height of > 10.7 m (35 ft) or volumes > 15 acre-feet]. Nine of the dams are classified as high risk, based on storage capacity, elevation, and location on highly fractured terrain. The dam heights range between 3.65 m and 13.10 m, with only one dam above the state-regulated height of 10.7 m. The maximum storage capacity of the high-risk dams falls between 100,000 and 580,000 gallons of water. These dams present both individual and cumulative hazard potentials to the population and infrastructure that could be impacted by their sudden failure.

Fig. 10
figure 10

Earthquake Hazards Map of Southeast Missouri (Marcus 1993) with an inset of a detailed map of the study area extracted for analysis. The study area shows severe liquefaction in the Mississippi embayment, while the collapse potential is projected along Interstate Highway 55. This zone falls within the high-volume Ordovician limestone, dolomite, and sandstone perturbed by high-angle, interconnected faults (see Fig. 1). Slope instability locally increases along channels where steep-sided banks are vertically undercut along tight radius bends

Sustained earthquake vibration can trigger liquefaction of saturated cohesionless soil materials, like rock flour, silt, and sand. Of great concern are the high impact areas where 10 oil and gas pipelines (green lines on Figs. 10 and 11) as well as 14 highway and railroad bridges (red squares on Fig. 11) (red squares on Fig. 11) founded on unconsolidated sediments comprising the Upper Mississippi Embayment. Therefore, hazard evaluations should consider credible sources of ground shaking in the area. Since the local network of seismographs were installed in the Upper Mississippi Embayment in 1973–1974 (including one in Cape Girardeau), only 29 earthquakes of between Mw 1.5 and Mw 3.1 have been recorded. The earthquake energy released increases 1024 times/Mw; therefore, the 29 events that occurred between 1996 and 2017 were capable of causing weak vibrations. While the mining and quarry activities from crushed stone, limestone, sand, and gravel trigger regular vibrations capable of producing topsoil liquefaction over time, the risk factor appears to be rather low.

Fig. 11
figure 11

Left—New Madrid Seismic Zone map showing major pipelines and bridges crossing within the Mississippi Embayment and adjacent to the Reelfoot Rift, compiled by David Hoffman at Missouri University of Science and Technology. Right—Map of the study area showing rock bearing faults (lines in red) that run almost perpendicular to the major oil pipelines. The southwest portions of the map in the Mississippi Embayment are at risk of liquefaction when high seismic loading occurs in the NMSZ (Rogers 2010)

Discussion

The vegetation categories in the study area were interpreted based on a few fundamental principles and rules of practice in remote sensing. The green-colored linear arrangement of vegetation patterns is usually interpreted as cultivated vegetation. The exposure of curvilinear bare-ground scars on hillsides and tilted trees are common manifestations of shallow earthflow landslides vegetated areas on Landsat orthoimages. We observed that between periods of active raveling and mass movements, the vegetation re-established itself with a significantly different vegetation signature—the native plants in landside areas with ephemeral springs and a high-water table were a deep green color. The study area receives between 63 and 122 cm of annual precipitation that can percolate through fractures, joints, or highly porous horizons, such as limestone with interbeds of shale and sandstone. Sustained precipitation can often result in flooding, causing liquefaction and slope stability issues that can culminate in shallow earth flows. These earth flows occur much more frequently than earthquake-related slides. When the slope analysis map was evaluated in conjunction with the USGS Earthquake Intensity Map and the geologic map of the study area, the areas most at-risk to slope failure also fall within the Mississippi Embayment, which is susceptible to ground movement during earthquakes. The cumulative effects of the presence of interconnected faults serve to cause local fracturing and brecciation of the fault zones, as well as local perching of groundwater compartments bounded by the old faults.

The limestone bedrock in southern Cape Girardeau exhibits a rapid response to flooding and supports a high-water table, which enhances the dissolution of the rocks. Figures 12 and 13a–f show the schematics of the stages of water-induced dissolution in a carbonate environment and examples of sinkholes in the Missouri Bootheel (Roberts 1964; Lolcama 2003, 2009; Palmer 2004; Lolcama and Gauffreau 2008; Bansa 2018). In 2007, the North–South BNSF Railroad line paralleling the Mississippi River, south of Cape Girardeau, was threatened by multiple sinkhole collapses adjacent to the Buzzi-Unicem Limestone Quarry. The impacted zone has been fitted with telemetered downhole monitors to measure settlement and provide real-time warnings to the railroad and local officials (Figs. 6e, 13f).

Fig. 12
figure 12

Source: Missouri Geological Survey)

Stages of sinkhole formation in carbonate environment; a Dissolution by groundwater circulating through the rock. b As the rock dissolves, spaces and caverns develop underground. c After the underground support for the topsoil has being dissolved, it loses support, and sudden collapse occurs. d Multiple collapse can result in shapes like shallow bowls or saucers, while others have vertical walls, and some hold water, forming natural ponds. e Inverse filter used in reversing the hydraulic piping in sinkholes (

Fig. 13
figure 13

a, b Sinkhole formation in the Missouri Bootheel region. c Sinkhole adjacent to the Buzzi Unicem Quarry and the BNSF Railway line. df Vertical section showing the early stages of sinkhole formation, as solution widening occurs

The spate of sinkholes was clustered around the northeastern corner of the Buzzi-Unicem Quarry (Figs. 13a–f), undermining the channel of lower Cape La Croix Creek, necessitating the demolition and replacement of the highway bridge carrying South Spring Street across the channel, and severing buried utilities along Spring Street. The City of Cape Girardeau also undertook a series of engineered mitigation measures to reopen South Spring Street and restore the operative discharge along the Cape La Croix Creek. In order to achieve this goal, the City undertook an extensive renovation of the creek channel and the South Spring Street Bridge. These corrective measures included reconstruction of 100 m of the channel’s right bank, just upstream of South Spring Street and closest to the northeastern corner of the Buzzi-Unicem Quarry.

An embankment of well-graded rockfill was constructed over a series of sinkholes that were over-excavated and backfilled with an inverse filter, to retard hydraulic piping along joint-bordered fissures (see Fig. 14). A more comprehensive correction was undertaken along South Spring Street, which began with the removal of the damaged bridge across the Cape La Croix Creek. In late 2007 to early 2008, the sinkholes in this area also began diverting seepage into the quarry floor, welling upward from the Cambrian-Precambrian contact at the base of the Plattin Limestone. The seepage increased to a level necessitating the backfilling of 15.3 vertical meters (50 feet) of graded rockfill to increase the confining pressure, which served to retard the flow from the base of the Plattin Limestone.

Fig. 14
figure 14

Modified from Deere and Patton 1971). a Collapsed sinkhole which has been infilled with residuum, forming a basin with perched water. This phenomenon is common in areas with slopes covered with unconsolidated materials. b A fresh sinkhole collapse which has not been filled. c An incipient sinkhole filled with layers of soft clay, stiff clay, and colluvium. The proposed base of cut is a line that shows optimum depth of excavation of the weakened carbonate layers in order to stabilize the failed karst zone. An inverse filter array is typically used to retard hydraulic piping through the open jointed rock

Schematic view of the reconstruction of a series of sinkholes in the same area (

Implications of failures

The present work evaluated the susceptibility of the study area under various scenarios and identified areas where risks of various geohazards could be mitigated to reduce the likelihood of future damage in those situations where avoidance was not realistic. The probability that a hazard may occur evaluates the outcomes of the minimum and maximum consequences of any such events on lives, health, and social infrastructure. The risk increases as the probability of occurrence is multiplied by the volume of damages (loss) if the anticipated hazard or failure mode occurs.

$${\text{Risk}} = {\text{probability}} \times {\text{loss}}{.}$$
(3)

Thus, evaluation of failure susceptibility measures the probability and consequence of failure of a given hazard or failure mode that poses a credible threat to lives or protective infrastructure, which often has the potential to quickly cascade the resulting damage by several orders of magnitude. For instance, a single landslide could sever linear infrastructure (highways, railroads, buried utilities, oil and gas pipelines, transmission lines, fiber optic cables, etc.) that could be injurious to sustainability of developed areas at considerable distance from where the slope failure occurred.

For the susceptibility to be insignificant, a landslide must have little or no impact on the surrounding population, which, depending on the density or development and size of the mass wasting, could impact anyone whose sustainability hinges on the interrupted arterial (e.g. loss of electricity, potable water, or internet connectivity). Conversely, extreme outcomes may be experienced in densely populated areas, where lives and sustainability are reduced as a result of landslides. In the past 15 years, much more attention has been paid to probabilistic hazard assessments that consider all the various failure modes that might impact the operability and sustainability of several engineered systems, like highways, bridges, dams, levees, and associated control structures. The number of “failure modes” (scenario events) has been increasing as each failure is critically evaluated (e.g. the rapid retrogressive erosion of the Oroville Dam spillways in 2017). By considering the entire panoply of potential failure mechanisms, a community’s disaster resilience can be bolstered substantially. A community’s resilience is dependent on its capacity to recover quickly following any sort of geohazard or disaster. It is worth noting that the majority of the seismic events recorded in the study area range between Mw 1.5 and Mw 3.6. However, Petersen et al. (2014) has warned that the  likelihood of a Mw 6.0 earthquake has a 35% probability of occurrence over the next 50 years.

One of the indirect benefits of increasing resiliency in one hazard area is often an unanticipated benefit in another area. A prime example of such crossover benefits would be the parallel aspects of seismic and wind loads, as could be expected from earthquakes or tornadoes. The engineering mitigations for both hazards are virtually identical (the only difference being storm shelters in light structures bereft of basements). The common mitigation techniques for slope stability (increasing subdrainage and flattening slope inclinations) increase the safety factors for intense precipitation (increased runoff or debris torrents), flooding (which includes driver safety and commercial riverine navigation), channel erosion and deposition, and earthquakes.

Co-seismic hazard potential susceptibility map

The gray zone in Fig. 10 highlights the erosional escarpments and slopes most susceptible to seismically induced landslides. The bedrock is characterized by high volume Ordovician limestone, dolomite, and sandstone perturbed by high-angle, interconnected faults. Accordingly, the factors that increase geohazard risks in the study area to damaging impact were the clusters of landslides in places that identified the potential for smaller earth flows to coalesce into larger landslide complexes (Frankel et al. 2009). The preponderance of crystalline limestone bedrock with interconnected fault patterns and much younger unconsolidated materials like wind-blown loess and alluvial cover leads to the susceptibility of the northwest quadrant. This area has 17 embankment dams that fall just below the jurisdiction of Missouri’s dam safety statutes, so the state does not exercise any oversight of their design or maintenance. As these structures age, they will likely pose increasing risks to downstream inhabitants and infrastructure. The failure consequences arise from dam storage capacities and locations along the same water courses. Agricultural dams are typically situated in areas of low population density. The largest threat category in Southeast Missouri appears to be liquefaction potential, which encompasses the entire Missouri Bootheel, except for portions of Crowley’s Ridge (the Benton and Bloomfield Hills).

The liquefaction potential is moderate to high because of the high-water table and preponderance of cohesionless sands beneath a surficial layer of overbank silts that create a semi-confined aquifer at shallow depth (Chung and Rogers 2015). The region also exhibits evidence of lateral spread features triggered by the 1811–1812 earthquakes (Fig. 2). These have been identified on the eastern margins of Holly Ridge near Bloomfield and just east of Advance, where Cape Girardeau, Bollinger, Stoddard, and Scott counties converge near the Bell City Gap (Watkins and Rogers 2009a, b).

One of the significant risks to flood control and agribusiness (mostly soybeans, wheat, corn, and rice) are the protective levees along the western side of the Mississippi River, which form Missouri’s boundary with Illinois, Kentucky, and Tennessee. With the exception of Crowley’s Ridge, the entirety of the Bootheel region lies below the flood water surface of the Mississippi River. The St. Francis River forms the western boundary of the Bootheel with the State of Arkansas, and its flow level at the Arkansas border is 70.1 m (230 ft) above sea level. This elevation is about 9.9 m (32.5 ft) below the Mississippi River’s normal flow level, along an east–west line. This means that a levee breach or failure along the right bank of the Mississippi River in this area has the potential to flood an enormous land area. There is, therefore, considerable risk in any sort of breach along the right bank levees, whether by earthquake damage or any alternative failure mechanism. The earthquake Factor of Safety depends on a number of variables, such as the river’s flow level (freeboard below the levee crest) at the time of the earthquake, surface fault rupture, as occurred along the Reelfoot Rift reverse fault segment on February 12, 1812, as well as the magnitude, duration, and areal distance from the earthquake hypocenter.

The other areas noted in Fig. 9 are considered as lower risk with less threat to the population, private properties, and public infrastructures. However, the preexisting landslides can be triggered by earthquake shaking, precipitation, or flood-induced erosion, either along channels, within alluvial valleys, along levees, dams, spillways, man-made embankments, or natural hillsides. Over-steepened erosional escarpments are also subject to undercutting, slope creep, and dilation, which can hasten slope instability, such as rockslides and debris torrents. The porous loess blanketing upland areas like Crowley’s Ridge are highly erodible once saturated. These conditions cause earth materials to reach the limiting slope angle for cohesionless materials like loose rock and sandy soils.

Conclusions

LiDAR and geospatial analyses of the northern Missouri Bootheel in the NMSZ reveal its susceptibility to a moderate risk of experiencing an array of geohazards, especially when sustained precipitation and earthquake strong motion cycles accompany high intensity shaking, as occurred during four significant events in 1811–1812.

Landslide susceptibility mapping was used to identify the most common mass wasting conditioning factors and to highlight locations most vulnerable to landsliding given the activation of one or a combination of factors. The study explored the extent of spatial, temporal, and geomorphological alterations measured through image differencing in zones of multiple landslides. Primary physical conditions, such as precipitation and hydrology, appear to increase the susceptibility of a hillside area to reactivation or generation of new landslides, especially along the southeastern erosional escarpment of Crowley’s Ridge within the NMSZ.

Native vegetation observed during field visits exhibited preferred associative responses to slope creep and/or fluidization of unconsolidated topsoil and colluvium (slope wash). These observations suggest that water is often trapped within colluvial-filled bedrock ravines infilled with rock dust (loess), which appear prone to absorption to free moisture and softening that promote the establishment of perennially wet swales that are generally unsuitable for agriculture. The permeability contrast between the highly porous loess and the underlying crystalline bedrock may aid in the sustenance of such landslide features. The larger landslides leave substantial erosional escarpments, and can threaten downslope infrastructure with rockfalls.

GWLR was then performed, and the relative risks were ranked based on potential impacts on people, property, and infrastructure. In general, the geohazard risks in most parts of the greater Cape Girardeau area are moderate, except for the risk of flooding in the watershed areas, and the steep hillside slopes within the sandstone, limestone and dolomite complexes that exhibit closely spaced faults and systematic joints.

Hydrological conditions can also play a crucial role in the formation of sinkholes in areas underlain by karst features that can serve as flow conduits when carbonate rock has been excavated through multiple semi-horizontal horizons. Flow conduits can allow seepage to descend with little loss on hydraulic head, allowing significant back-pressure to develop above semi-impervious horizons. These situations can hasten the formation of new sinkholes in locations that have not previously been prone to developing sinkholes.