Predicting the location and timing of mudslides with adequate lead time is a scientifically challenging problem that is critical for mitigating landslide impacts. Here, a new dynamic modeling system is described for monitoring and predicting storm-triggered landslides and their ecosystem implications. The model ingests both conventional and remotely sensed topographic and geologic data, whereas outputs include diagnostics required for the assessment of the physical and societal impacts of landslides.
The system first was evaluated successfully in a series of experiments under idealized conditions. In the main study, under real conditions, the system was assessed over a mountainous region of China, the Yangjiashan Creeping (YC) slope. For this data-rich section of the Changjiang River, the model estimated creeping rates that had RMS errors of ∼0.5 mm yr−1 when compared with a dataset generated from borehole measurements. A prediction of the creeping curve for 2010 was made that suggested significant slope movement will occur in the next 5 years, without any change in the current precipitation morphology. However, sliding will become imminent if a storm occurs in that 5-yr period that produces over 150 mm of precipitation. A sensitivity experiment shows that the identified location fails first, triggering domino-effect slides that progress upslope. This system for predicting storm-triggered landslides is intended to improve upon present warning lead times to minimize the impacts of shallow, fast moving, and therefore hazardous landslides.
Because landslides are movements of Earth’s surface, they illustrate many of the interconnected processes that influence landscape stability. The locations of landslides depend on local topography, geological composition, forest cover, soil water content, precipitation morphology, and seismic activity. For the majority of landslides, heavy rainfall events and earthquakes are the main triggers. The former is far more frequent than the latter, so this study focuses mostly on storm-triggered shallow landslides (Iverson 1997; Iverson 2000). Fatalities from landslides that followed the 2008 Wenchuan earthquake demonstrate how increasing populations in high-risk settings have enhanced the threat of this natural hazard and accordingly have increased the urgency to develop and implement effective early prediction and warning systems.
Landslides occur irregularly, so landslide studies concentrate on predicting when and where landslides will occur; predicting how large they will be; and developing procedures that convey warnings to the public to mitigate damage to infrastructure and ecosystems. Empirical and descriptive landslide models have contributed much to the public awareness of landslide hazards and have led to valuable accumulated experience in identifying the key causal factors. A synthetic consideration of preparatory and triggering factors, however, demands a comprehensive modeling of the physical processes involved in landslides (Costa 1984; Iverson 1997). The overview by Iverson (Iverson 1997) suggested several criteria for dynamic landslide models, including that a model should be capable of simulating the full start–movement–spread–cessation cycle of the detached material and should cover a wide spectrum of debris flows. With continued growth and expansion of human population, rain-triggered shallow landslides increasingly result in loss of life and significant economic cost. From an ecological viewpoint, landslides are an important factor in desertification over mountainous regions, because they are very effective in transferring biomass from live to dead respiring pools (Ren et al. 2009).
Forecasting landslide timings and impacts has been a significant research topic in geotechnical engineering over the last four decades. However, because of the complexity of natural conditions (e.g., geometrical and geological variability, nonlinear time-displacements relationships, and the superposition of seasonal effects) and of ever-changing human activities, it is difficult to predict precisely the timing of slope failures (Bhandari 1988). Geoengineers typically use instruments such as the extensometer, distortometer, inclinometer, and GPS to monitor slope displacement. Based on these monitoring data, a mathematical model is applied to predict future slope conditions. The soil slope creep curve is usually divided into three separate regimes: strain hardening, steady, and accelerating stages (e.g., Saito 1965). The models typically employ strain displacement as a metric to describe slope behavior leading to failure (Fukuzono 1985; Voight 1988; Voight 1989). Similarly, Crosta and Agliardi (Crosta and Agliardi 2003) proposed a nonlinear viscoplastic bulk behavior slope model. Recently, gray models and artificial intelligence (AI) approaches, including neural networks, have also been used to forecast the stability and timing of slope failures (Sakellariou and Ferentinou 2005).
2. A numerical modeling system
As reviewed by van Asch et al. (van Asch et al. 2007), current methods used for landslide analysis are highly varied. Empirical or statistical techniques generally are applicable to predicting landslide susceptibility at the regional scale, whereas more process-based approaches are applied at the local scale. This study uses an advanced model, namely, a Scalable and Extensible Geo-fluid Model of the Environment (SEGMENT; see Figure 1), which accounts explicitly for soil mechanics and root distribution and has a complete cycle of hydrological processes and general applicability to regional topography, soil thickness profile, and vegetation coverage. The sliding material can be a mixture of mud and material it has entrained. This numerical modeling system provides prognostic fields of the driving and resistive forces and describes the flow fields and the dynamic evolution of thickness profiles of the medium under consideration.
As a model of mudflow, SEGMENT has been used successfully to simulate the sensitivity of storm-triggered landslides to deforestation and basal sliding (Ren et al. 2008b; Ren et al. 2009). The landslide components are outlined briefly in the appendix, which includes the governing dynamic and thermodynamic equations, some unique to SEGMENT, and the numerical solution techniques, to establish a basis for discussion of recent modifications and updates. In summary, the soil mechanics component employs a strain–stress relationship governing granular rheology. It is a set of conservation laws that determine the time evolution of near-saturated soil slab in response to gravitational driving stress, failed slab strength, changes due to soil moisture increase or decaying of the root system, and basal sliding initiated by earthquakes. Because it retains the inertial terms of the momentum equations, the landslide model simulates the full cycle of start–movement–spread–cessation of the detached material.
Some model components recently have been improved. These include updates that allow for vertical stretching of the model grids, augmentation of the terrain-following σ coordinates to facilitate imposing of lateral boundary conditions, and an alternative biogeochemical submodel [Terrestrial Ecosystem Modeling System (TECO)-Nitrogen; Weng and Luo 2008] has been added as an option to the Collaborative Adaptive Sensing of the Atmosphere (CASA; Potter et al. 1993; Ren et al. 2008a).
In the remainder of this section, two interrelated SEGMENT components are described: a sophisticated land surface model to simulate the soil moisture regime and a well-tested biogeochemical submodel that budgets carbon and nutrient cycling in the ecosystem (see Figure 1a). Both components provide estimates of the ecosystem consequences of landslides.
2.1. Land surface component
To understand changes in mudslide frequency, geographical distribution, and magnitude, the greatest influences are variations in precipitation amount and intensity and how much either is infiltrated into the soil or is run off. This requires a land surface scheme that balances computational load and accuracy. The land surface model component, which governs surface runoff routing and water infiltration, is adapted from the Simulator for Hydrology and Energy Exchange at the Land Surface (SHEELS) and is described by Smith et al. (Smith et al. 1993) and Crosson et al. (Crosson et al. 2002). The Richards’ equation (Ren 2005) for soil moisture redistribution in the unsaturated zone is modified according to Iverson (Iverson 2000) for applications to sloping ground surfaces. In addition, the modified SHEELS accounts for macropores in the infiltration–drainage scheme. The macropore features are specified according to detailed geological surveys. Thus, in the simulation domain, the infiltration coefficients vary three dimensionally. With the macropore scheme, the sliding material is not required to be saturated uniformly down to the sliding surface to trigger landslides. Only when the shear zone granular material is sufficiently wet will landslides occur in SEGMENT. In general, the macropores allow a portion of the runoff water to be channeled into the shear zone granular material and contribute to enhanced basal sliding.
Fully three-dimensional simulations of mudslides, especially fast-moving mudslides with Reynolds numbers up to 1000, require small integration time steps and fine spatial resolutions. Such simulations are computationally highly demanding. Some of the lower-order physics terms in the land surface model component therefore have been simplified (Ren et al. 2009). For example, a computationally inexpensive soil temperature solver is used, the additional stress on the surface exerted by falling raindrops is not included, and an economical simplification approximation is adopted for evapotranspiration (ET) processes. There are improved features in the land surface component, such as thermomechanical and hydraulic parameterizations of organic soil (Lawrence and Slater 2008; Letts et al. 2000). The drag from Horton overland flow also is implemented, because it is critical in determining the potential for shallow landslide/mudslides and whether they will develop into catastrophic landslides.
Once sliding is initiated, the inertial (local acceleration and advection) and the viscous forces act against the gravitational driving force. Regardless of the constitutive relationship present, the maximum attainable kinematic energy depends exponentially on the slab thickness. This explains why the spatial distribution of soil depth exerts a strong control over local landslides (Casadei et al. 2003). In mountainous areas, besides the infiltrated precipitation amount, surface soil availability limits the saturated slab thickness. In initializing the maximum sliding mass thickness in SEGMENT, both the soil and vegetation maps are referenced. This applies only at the initial time level, and the subsequent soil thickness evolution is fully governed by the model simulation [Equation (A.4)]. Hence, as the simulation proceeds, the soil thickness can become positive for the impermeable surfaces by advection from neighboring grids.
For plants, large canopy sizes are advantageous in competing for light and CO2 but are a liability for wind stress or ice loading, because they contribute to stem break and hence to landslides for vegetated slopes. Similarly, thick stems are more resistant to wind stress, but the additional loading contributes to slope instability. Plants balance their growth to defend against extreme events against their competition for nutrients and water. Natural selection ensures that existing species have strong or flexible stems and roots to cope with wind throw and landslides under a naturally varying climate. Root strength (>0.3 MPa for most chaparral shrubs) plays an important role in slope stability, especially for vegetated slopes composed of coarse colluviums with cohesion values normally less than 0.1 kPa (Table 1). The tree root geometry is assumed to be vertically downward. The power-law root tensile strength–root diameter relationship (Gray and Barker 2004) is implemented, with a plant-species-dependent power index and reference root strength.
2.2. Biogeochemical component
Soils form a thin film over Earth’s surface and are located where geological and biological processes intersect. As the medium in which most decomposition organisms live, soils play a key role in the global cycles of carbon and nitrogen. The physical soil matrix is a source of water and nutrients to plants and microbes and is the support system for terrestrial vegetation. The physical and chemical characteristics of soils influence all aspects of ecosystem functioning.
Landslides are natural disasters that displace the soil medium, leaving scar regions with near- or nonexistent soil. As disturbances to ecosystems, landslides transfer live biomass into dead organic pools, as a result of plant mortality. In SEGMENT, vegetation destruction and burial occur when maximum sliding speed and surface elevation changes exceed vegetation-dependent critical values. Normal decomposition breaks down dead organic matter, releasing carbon to the atmosphere and nutrients in forms usable in plant and microbial production. As dead biomass decomposes, different parts decompose at different rates. For example, leafy and fine roots take days to decompose under warm moist climate conditions, whereas the denser, woody parts take up to 20–30 years to diminish to insignificant amounts. As litter decomposes, its mass decreases approximately exponentially with time, so there is a turnover time scale associated with different biomass pools in the numerical model. Currently, nitrogen cycling remains an open question. Soil moisture, soil temperature, ventilation condition, and the activity of nitrifying bacteria control both the quantity and quality. Among them, nitric oxide and nitrous oxide, which are produced during nitrification, are gases that have important effects on atmospheric chemistry. TECO-Nitrogen was chosen primarily because it budgets carbon and nutrient cycling in the ecosystem in a more sophisticated manner than just using the C:N ratio.
Ren et al. (Ren et al. 2008b) have described the model’s governing equations and stress–strain decomposition. Details of the lower interface entrainment process include an improved treatment of the entrainment process that enables the model to describe the progressive increase in surface runoff (Cannon et al. 2001) for freshly littered ground, such as recently burnt areas. Analogous to Morton et al. (Morton et al. 1956) for jets and plumes, the rate of entrainment is assumed to be proportional to the mass mean flow speed and the effective density of the flow material and debris to be entrained. The entrainment constant is parameterized as a function of soil moisture content and shear magnitude near the lower boundary. Because entrainment accounts for only a small proportion of the deceleration, the low Reynolds number approximation reconciles several approaches found in the literature. The entrainment process is especially pertinent in describing rain-triggered landslides in high debris-accumulation circumstances (e.g., the fire scars left by Southern California wildfires).
The model is applicable to regional topography and ingests all data required by a sophisticated land surface scheme and, in addition, requires a vegetation weight mask and a root distribution profile. The key datasets for the development and verification of this distributed dynamic landslide model include
the Tropical Rainfall Measuring Mission (TRMM)–based Multisatellite Precipitation Analysis (TMPA3B42RT; Huffman et al. 2007; available online at http://trmm.gsfc.nasa.gov; of 0.25° horizontal resolution) real-time rainfall data;
SRTM-derived hydrological parameter files [Hydrological Data and Maps Based on Shuttle Elevation Derivatives at Multiple Scales (HydroSHEDS); Lehner et al. 2008; available online at http://hydrosheds.cr.usgs.gov];
soil parameters provided by the Food and Agricultural Organization (available online at http://www.fao.org/AG/agl/agll/dsmw.htm); and
the Moderate Resolution Imaging Spectroradiometer (MODIS) land classification map, which is used here as a surrogate for land use/cover, with 17 classes of land cover according to the International Geosphere-Biosphere Programme classification (Friedl et al. 2002).
Also, the University of Maryland 1-km global land cover (Hansen et al. 2000) can be used to characterize vegetation distribution. The morphological characteristics of the root system (e.g., root distribution with depth, root distribution over different diameter classes and root area ratio, root tensile strength and tensile modulus values) are assigned according to land-cover type and a remotely sensed above-ground biomass density map (Zhang and Kondragunta 2006). Below is a description of the datasets for generic, mountainous regions of Earth and also for data-rich regions.
3.1. Data for generic mountainous areas
Satellite-based geospatial datasets are valuable for cost-effective detection and early warning of storm-triggered mudslides, particularly for regions of the world where in situ measurements are sparse or nonexistent. As an alternative to the ground-based observations, timely satellite-based rainfall measurements now are available and are virtually uninterrupted. It is best to start the model integration with remotely sensed precipitation. Simulations with only observed precipitation may provide early warning of mudslides. However, for a truly predictive system, future precipitation usually is provided by a numerical weather prediction (NWP) model. For indications of imminent landslides, simulations with forecast precipitation amounts should be performed to provide more accurate predictions of the timing of storm-triggered landslides.
The landslide modeling system requires ground surface loading information. Ideally, this is the above-ground biomass. Unfortunately, in situ measurements are not globally available. Compiling and updating the traditional inventory data also is highly time and resource consuming (Brown et al. 1999). Zhang and Kondragunta (Zhang and Kondragunta 2006) estimated forest biomass for the United States using generalized allometric models and remotely sensed vegetation canopy foliage products from MODIS. The above-ground biomass map (Figure 2c of Zhang and Kondragunta 2006) can be used as a surrogate for the vegetation weight mask. This dataset has continental coverage and the potential to be extended to global coverage (X. Zhang 2009, personal communication).
For model verification, data-rich regions such as California are best. California has a forest biomass with a very high spatial resolution (∼30 m). This dataset was generated by combining the National Forest Inventory Analysis Program (FIA; available online at http://www.fia.fs.fed.us)–derived timber volume estimates with a forest cover map (Franklin et al. 2000). The timber volume data were then converted to biomass values using an expansion factor coefficient. For a small region along the Qingjiang River, a tributary of the Changjiang River in China, the so-called Yangjiashan Creeping (YC) slope, located at 29°50′N, 109°14′E (Figure 2), covers an area of only ∼20 km2, and a geological survey team collected the necessary biomass data for the study in section 4.2, following a request by the fifth author.
As input, SEGMENT requires geomechanical parameters related to soil/rock material property, including cohesion, dry repose angle, density, porosity, field capacity, saturated hydraulic conductivity, and initial bedrock topography. A global landslide hazard evaluation by Nadim et al. (Nadim et al. 2006) assigned material strength properties to the land surface. Larger-scale geologic mapping is readily available for the continental United States (King and Beikman 1974; Schruben et al. 1998) and Alaska (Beikman 1980). Outside the United States, the geological database was compiled from existing 1:5 000 000 scale maps from the United Nations Educational, Scientific and Cultural Organization (UNESCO; Hearn et al. 2003). The Nadim et al. (Nadim et al. 2006) ranking is modified to include several units in the highest susceptibility category and to assign friction angles and cohesion based on published values (Selby 1993; Jibson et al. 2000).
3.2. Data for a surveyed slope
Detailed geological surveys are conducted only for areas with important infrastructures in danger from natural hazards. Fortunately, one such geophysical dataset is available for YC (Figure 2b). In addition to a detailed, current geological map, an engineering study was performed of the soil and rock profiles by means of borehole (BH) drilling, extracting soil samples for laboratory testing (Table 1), and continuous measurements directly within the creeping zone (Figure 2b). The shearing strength properties of the various materials present also have been obtained. Both drained and nondrained shearing tests are carried out on the soil specimens (Table 1). In addition to recovery creeping zone rock and fluid for laboratory analyses, intensive down-hole geophysical measurements and long-term monitoring provide direct information on the composition and mechanical properties of active creeping zone rocks, the nature of stress responsible for sliding, and the role of pressurized water in controlling landslides recurrence for field conditions with a wide spectrum of dynamic ranges. Displacement surveys have been carried out continuously since July 2007, using the RST-IC3500 digital inclinometer. Displacement data from five boreholes (the labeled red circles in Figure 2b) are analyzed in this study: BH7, BH8, BH8–1, BH9, and BH10 (see Table 1).
The YC region has a subtropical moist climate with a typical monsoonal precipitation pattern. The annual-mean precipitation ranges between 1100 and 1700 mm. This section of the Yangjiashan Mountain, because of its proximity to a dam, is well instrumented and is an ideal region to verify the numerical model. For example, there are in situ measurements of material properties as well as the meteorological parameters required for the land surface model component. The model directly evaluates the roles of pore pressure, biomass loading and root distribution, and intrinsic friction in stress distribution. In this region, the major slides are preceded by creeping movements. Since the 1960s, the YC slope has had repeated failures, as a result of exceptionally heavy rainfall periods (e.g., in 1960, 1980, and 1998). It appears that after major rainstorms the toe-slope failure occurs first, reducing the stability of the upper portion, and the failure moves gradually to the upper slopes. The data collecting schemes therefore were designed to optimally extract information for a toewise wedge.
4. System verification
The model is verified first for idealized conditions and then for the YC region. For the data-rich YC, our modeling results are verified against borehole measurements of creeping rates. A prediction also was made of the next likely serial sliding cycle for this slope. The YC is not representative of most areas, because it has detailed in situ measurements of geological and meteorological parameters. Applications to most regional areas rely on the coarse-resolution meteorological, land surface vegetation, and geological data.
4.1. Verification with idealized environmental settings
The numerical model was verified using idealized topography and a soil thickness profile (Figure 3). A uniform sandy soil layer (d = 0.75 mm, repose angle ϕ = 30°) of 1.1-m depth is assumed, with a soil moisture content of 35% relative saturation and spread over two identical bell-shaped mesa domes located at x = 210 m, y = 210 m (M1 in Figure 3a) and at x = 840 m, y = 840 m (M2 in Figure 3a). The two domes have peak elevations of 2400 m and e-folding distances of 7500 m. A cylindrical boulder (x = 450 m, y = 840 m; see the hatched area in Figures 3a–f) of 300-m diameter is buried 0.5 m underground. The surface elevation changes were examined, which also is soil thickness change because the bedrock is assumed unchanged, at 15, 30, 180, 500, 1000, and 2000 s after the onset of sliding. Velocity fields in the midlevel (initially at depth 0.55 m) also are overlaid (vectors). As a result of the dynamic drag effect from the cylindrical obstacle, there is an upstream accumulation of mass and a downstream mass sink (Figure 3b). The ravine between the two tolls (the diagonal line from top left to bottom right) accumulates mass while the two domes lose mass. The soil located above the 2300-m elevation contour is depleted at 500 s after flow onset (Figure 3d). The flow fields are physically reasonable. A steady flow field is formed at 10 s. Because of the boulder’s obstruction, flow is concentrated in the lower-right section and is divergent west of M2. As mass is depleted, the flow weakens (cf. Figures 3c,d). This feature starts at the peaks and propagates down the slopes. As a result of the no-slip lower and zero-stress lateral and upper boundary conditions, the speeds at depth are weaker than at the surface. A zero-stress lateral boundary condition also prevents soil accumulation at the boundaries. Determined primarily by the shear zone depth (1.1 m) and the slope steepness of the bedrock, the slow slope movement lasted 50 min before ceasing, with the isolated areas of mud accumulation around the obstacle being depleted last (cf. Figures 3e,f). Over 80% of the total mass loss occurred during the first 10 min, though. The full life cycle of the idealized mudslide therefore has been well delineated.
With idealized slide material properties, basal sliding conditions, soil moisture content, and vegetation loading and root distribution, the adaptability of our modeling system was tested on a range of experiments including a slide–stop cycle for material scoured by a surface wave; simulation of a slide–slip–stop cycle for an originally stable but perturbed slope (finite period basal sliding with given magnitude); the effects of vegetation for various soil moisture regimes; and the effects of vegetation removal on slope stability.
Finally, in a series of sensitivity studies using the model, it was found that soil moisture direct loading from precipitation on the slope is miniscule. However, the soil moisture increments drastically change soil rheological properties. As the soil is saturated to a depth deeper than most roots, the preventive effects of the rooting system (Sidle 1992) are diminished. For all vegetated surfaces, canopy runoff is an effective preventive factor for landslides. For forest cover, biomass loading and wind stress on canopy are important factors for triggering slope instability, particularly when pore pressure builds up between sliding material and bedrock. In addition, also tested were various combinations of slope geometry, shear zone thicknesses, soil moisture conditions, and slide material compositions. SEGMENT can simulate the five principal types of slope movement, as categorized by Cruden and Varnes (Cruden and Varnes 1996).
4.2. A decadal prediction for YC slope
The YC community is located on hills composed of interbedded silt stones and sandstones, with occasional altered clay layers. The protruding weathered stones are from the Paleozoic era, the 200–900-m-thick yellowish sand stones/silt stones are from Silurian period, and the gray silt stone from the Triassic period is especially sensitive to pore pressure changes. Hence, infiltration of rainfall through macropores, which already are well developed in the area, plays a critical role in slope stability. The hills intersect with canyons in which increased erosion takes place during the spring and fall rainy seasons. Although many of the drainage patterns in this region have been altered by geoengineers to increase slope stability, some remain unchanged, even the inhabited areas east of the demarcation line in Figure 4a. The existence of the YC landslides has been known for decades. During a storm in July 1998, very heavy rainfall fell over several days and caused substantial erosion in one small canyon. A field study revealed that a toe-slope failure was the most likely triggering factor activating the slide in the upper slope area. Near the end of a rain storm, a block of material was undercut by the stream and moved into the canyon, its movement downhill leaving an unsupported upslope block, which followed its movement. This, in turn, was followed by a third block as the movement progressed uphill. Although there were no casualties, many homes were badly damaged. Because of the detailed data available for this region, the model can be verified over this small region with the borehole-monitored creeping rates.
The model was initiated using five borehole soil/rock profiles (18–50-m depths). A horizontal resolution of 10 m was used for the 3000 m by 5500 m simulation domain, while fully encompassing the region shown in Figure 4a, to reduce possible boundary effects. The sliding material is thicker downslope than upslope. A thin plate splines (TPS; Burrough and McDonnell 2004) interpolator was used to obtain a sliding material depth distribution over the grids. Three layers were used to delineate the surface soil layer, and one layer was used to represent the nonfractured montmorillonite layer and a nonslip bottom bedrock layer. From the profile data (Table 1), layers 9–12 form the slip surface, which is confirmed further by the creep monitoring data (Figure 4b). This zone has the same chemical composition as the underlying zone, but physically it is fractured. There are four sublayers delineating this shallow zone because, although only ∼1 m deep, it is the lubricating layer. Because it represents granular material under high confining pressure, the viscosity of this thin bed of finer-grained materials is smaller than for the adjacent layers. The aforementioned macropore scheme is activated for the YC sliding slope simulation. Ground surface biomass loading information for YC, because of its limited coverage (only ∼20 km2), was collected by a geological survey team.
Under mean annual soil moisture conditions from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis (available online at http://www.esrl.noaa.gov/psd/data/reanalysis/), a creeping rate was simulated near BH8 of 17, 17.2, 16.2, 8.5, 5.2, and 3.1 mm yr−1 at 1.5-, 3-, 10-, 30-, 40-, and 45-m depths, respectively, agreeing well with the measurements. The RMS error, when compared with the 108 measurement data grids, is only a 0.42 mm yr−1 overestimation, within the instrumental error range. The depth of the sliding surface, 49 m at this location, also is accurately delineated. Sliding eventually will develop on this plane of weakness composed of highly fractured sandy shale stones. At the time of writing, measurements are current only as of 27 March 2008. A prediction was made of the 1 January 2010 displacement curve based on the climatological precipitation over the region (Figure 4b). The basal sliding from the Wenchuan earthquake of 2008 also is taken into consideration. Compared with the initial measurements, the shallow level displacements are as large as 18 mm. The creeping soon evolves into the accelerating stage. The Wenchuan earthquake reduced the natural creeping period by at least 5 years. Even with little change in precipitation morphology, the next 5 years may see significant movement of this slope. If the region experiences a storm producing over 150 mm of precipitation, sliding is imminent.
Three significant historical YC slides are separated by about 20 years (Figure 4c). The years 1960 and 1998 both had well above average annual precipitation amounts of 1819 and 1771 mm, respectively. The year 1980 was relatively dry, with 1200 mm of total precipitation compared to 1600 and 1360 mm for the years 1979 and 1981, respectively. Fortunately, daily reanalyses of precipitation began in 1979. Examination of the daily precipitation series from 1979 and 1980 indicates that 1979 had over 90% of its precipitation in the latter half of the year (with no significant precipitation before June 1979). Although the total precipitation for 1980 was small, the heavy precipitation in January immediately followed the previous year’s heavy precipitation formed an extended wet period. A significant precipitation event occurred in January 1980 (Figure 4c), reaching 58.3 mm day−1 intensity. As a result, deep soil moisture (0.23 volume per total volume) was relatively high over the ensuing several months. The trigger for the landslide was a heavy precipitation event of 100-yr recurrence frequency, a single day precipitation of 230 mm on 11 June. Higher soil moisture was a shared feature of all three major slides (red arrows in Figure 4c). For YC, consecutive storms with heavy rainfall constitute dangerous scenarios for triggering landslides. However, the observational data are 18 months long, too short to unambiguously confirm this hypothesis. Model simulations indicate that, under long-term mean soil moisture conditions, the creeping will reach a significant rate (3 cm yr−1 near the surface near the toe of the slope) in about 37 years. Interestingly, this intrinsic frequency may not have been realized during the past half century. Because of frequent interruptions by storms, landslides have occurred ahead of the natural schedules.
For a storm of 50-yr recurrence frequency (∼170 mm day−1), more than 80% of the total water mass is channeled to the sliding surface through macropores and the sliding can become dramatic toward the end of the storm. The areas of significant deformation (i.e., elevation changes greater than 0.5 m) are shown in Figure 4a. The escarpments are distributed more or less parallel to the elevation contour lines. So far, the geological survey team has identified about 10 scars and slide terraces (green crosses in Figure 4a). For example, there are obvious scars near BH8 and another major one near the position labeled point A. A shear surface exists around point A (dimension ∼70 m). When water drains down this surface, the material strength is reduced to its residual value. The artesian pressures along the failure surface add to the instability of the sliding mass. Point C was inaccessible to the survey team, because of its steep topography. However, the results obtained from the model indicate that point C is highly unstable under such an extreme precipitation event and may act as the first domino in a major sliding event. The sliding material can accumulate a thickness up to 20 m at the lower edge. The lowest sliding mass spread about 50 m down slope and was brought to rest by the lateral stress from the walls of the V-shaped canyon (parallel to the 750-m elevation contour). The canyon usually is dry, but the sliding mass may block the discharge of rainfall and enhance the infiltration of rain into deep layer and cause pore-water pressure rise, a lubricating effect in the model parameterization. This is especially important for the land segment lying between the 700- and 600-m contours, which contains the primary residential areas. To obtain additional information, several additional boreholes are required in the rock slopes adjacent to the canyon on the west bank.
For the above scenario, a volume of about 6.3 × 107 m3 of soil and rock is estimated to have been involved in the slides. It is likely that landslides have occurred in this region at intervals through history, but only part of the material is involved in a particular landslide. Specifically, one portion may slide, causing a reduction in the stability of an adjoining portion. Then, years, decades or even centuries later, a land slide will occur. As a consequence, the topography of the area is hilly and uneven.
The creeping rate of a slope is determined by the material viscosity. The modeled viscosity change was examined near BH8 (not shown). Because the rocks are heavily weathered, the viscosity is on the order of 1016 Pa · s−1, about two orders of magnitude smaller than the unfractured same material. So, when dealing fractured rocks as a whole, they need to be viewed as granular material lumped together. The viscosity of granular material changes as strain accumulates. In this case, the viscosity reduced substantially with time at all levels, reducing by >40% from 2007 to 2008. This explains why the creeping tends to accelerate with time.
A scalable, extensible, and adaptable geofluid model (SEGMENT) is described here and was developed to improve the landslide disaster mitigation capability of relevant agencies. It can be used either as a warning system for possible landslide hazards or as a predictive tool when precipitation rates are provided by a weather prediction model. Future operational use of SEGMENT hopefully will aid interested agencies in decision making for landslide damage mitigation.
Different regions have different requirements for early warning systems. For regions in active earthquake zones, such as parts of China and the west coast of the United States, sensitivity experiments such as those with prescribed basal sliding velocities are desirable. For ENSO-affected areas such as California and Peru, sensitivity experiments with respect to soil moisture content may need to be performed, one small region at a time. The combined figure can be used as a likelihood of occurrence map or as a levels of safety map.
The numerical model was verified for the Yangjiashan Creeping (YC) slope in China. YC’s special geological constituents are the primary cause of the landslides, and storms are an enhancing factor. This study suggests an additional process that needs increased monitoring, because it is a likely cause of the cascading of storm-triggered landslides; it was found that increased lubrication effects from storms enhances creeping rates dramatically for weathered slopes. For slopes experiencing repeated failure–restore cycles, increased precipitation amounts and intensities under a future warming climate are the two most important factors determining the long-term increase in landslides frequency.
The modeling system presented here was designed to estimate the potential for landslides over a regional area, rather than for a single slope. It was applied initially to an idealized situation for calibration and verification, then to the YC slope region of China. The promising performance of the model in the applications is attributed to the use of a new, fully three-dimensional modeling framework based on a newly proposed granular rheology and to the use of a land surface scheme that explicitly parameterizes the hydrological characteristics of macropores. Some requirements of the model, such as vegetation loading and root distribution in soils and weathered rocks, generally are not available, even in modern geological maps. Extension of our study to other regions is limited primarily by the lack of high-resolution input datasets.
In summary, the landslide component of the SEGMENT explicitly accounts for soil mechanics, vegetation transpiration and root mechanical reinforcement, and relevant hydrological processes. It has a nonlocal dynamic balance of the three-dimensional topography, soil thickness profile, basal conditions, and vegetation coverage in determining the prognostic fields of the driving and resistive forces, and it describes the flow fields and the dynamic evolution of thickness profiles of the medium being considered. It therefore can answer the following critical questions for storm-triggered landslides: How will a particular landslide be initiated? How fast will it move? How far will it travel? How large will it be? Together with a land-use map, the model can predict how much sliding material will be shed onto nearby infrastructures.
This study is supported by start-up funding provided to the second author. Shaikh Muhammad, Marc Wiederspahn, and Effie Jarrett from Jackson School of Geosciences, the University of Texas at Austin, helped with technical aspects of the study.
Partial Differential Equations of the Modeling System
Developing a continuum mechanics framework based on the laws of conservation of mass, momentum, and constitutive laws has proved to be challenging for natural soil slabs (T. W. J. van Asch 2008, personal communication). However, it is now possible to pursue a detailed landslide tendency study using a physically sound model applicable to regional scales rather than just individual hill slopes. Ren et al. (Ren et al. 2008b) recently proposed a new dynamical model for mudslides that is one integral component of a larger, more sophisticated multiple-phase scalable and extensible geofluid modeling system (Figure 1a). In the following, the focus is on its landslide component.
A.1. Landslide model component of SEGMENT
Soil rheology in general is complicated. Under different soil moisture conditions, the same soil is brittle, elastic–plastic, or viscous. For a storm-triggered landslide study, the situation is simplified because only the (near-) saturated moisture range is of concern. Under this soil moisture range, most soils behave viscoelasticoplastically. The near-saturated soil constitutive relationship of the SEGMENT is based on Terzaghi’s effective stress principle (Fredlund and Rahardjo 1993). The uniqueness of the critical state line in the deviatoric stress-effective mean stress plane for saturated/near-saturated soil is the theoretical basis of the SEGMENT-landslide model. The tensile strength and shear resistance of wet granular materials are closely related to soil moisture content. For swelling clay and quick clay, dampening or enhancement factors need to be further applied.
A soil dynamics model is proposed, based on the incompressible Navier–Stokes formulation. In standard notation, the mass conservation equation is
and the momentum conservation equation is
where ρ is bulk density, V is velocity vector, σ is internal stress tensor, and F is the body force (e.g., gravity ρg). Instead of the standard decomposition of the full stress tensor σ into static and dynamic stress parts, here we decompose it into lithostatic (L) and resistive (R) parts (Van der Veen 1999): σij = Rij + δijL, where Rij denotes the components of the resistive tensor and δ is the Kronecker operator.
Applying rheological relationship, the internal stress tensor term is expressed explicitly in velocity. Among several options, we choose the constitutive relationship proposed by Jop et al. (Jop et al. 2006) because it has proven to be acceptably accurate for studying the flow characteristics of real landslides (Ren et al. 2009) and therefore should provide good estimates of total mass scoured out of the scar and overall changes in the surface elevations. According to our decomposition of the stress, the granular viscosity is parameterized as
where ν is viscosity; S = [Rkk − ρg(h − z)]/3 is the spherical part of the stress tensor σ; μ0 and μ1 are the limiting values for the friction coefficient μ; is the effective strain rate and ; I0 is a constant depending on the local slope of the footing bed as well as the material properties; and I is inertial number defined as , where d is particle diameter and ρs is the particle density. The particle size is assigned according to soil type. Soil moisture enhancement factor on viscosity is assumed varying according a sigmoid curve formally as Equation (9) of Sidle (Sidle 1992) but with the time decay term replaced by relative saturation. We further specify c = 0, a = 1, and b as a soil type–dependent parameter varying between 0.3 and 0.95. The mechanical parameters of the bedrock are assigned according to the geological map (Nadim et al. 2006; Schruben et al. 1998; King and Beikman 1974).
The prognostic equation for surface elevation [h(x, y)] is
where X|top indicates evaluation at the free surface elevation. In the case with slope movements, Equation (A.4) is solved regularly to update the sliding material geometry.
Sliding can cause intense heating inside the material, especially near the bottom. This can be a positive feedback mechanism that accelerates sliding through further reduction of the shear band of loose granular media material strength. The thermal equation used is
where c is heat capacity (J kg−1 K−1), T is temperature (K), κ is thermal conductivity (W K−1 m−1), and σeff is effective stress (Pa). The last term is “strain heating,” converting of work done by gravity into heat affecting the sliding material by changing viscosity or causing a phase change.
The dynamic [(A.2)] and thermodynamic [(A.5)] equations are coupled in SEGMENT-landslide. The upper boundary condition for Equation (A.2) is that of free stress. The zero-velocity gradient is applied as lateral boundary condition. The boundary conditions for Equation (A.5) are a surface temperature at the top of each grid [T|z=zt = Ts(x, y)] and a geothermal flux to the base of the lithosphere layer (k(∂T/∂z)|z= zb= G0, where G0 is geothermal heat flux in W m−2). The above governing equations are capable to describe the episodic movement of unstable slopes, how large a storm-triggered landslide will be, and the deposition thickness and range.
Numerical solving of this system is not trivial. The finite differences technique is used to discretize the governing equations. The grid stencil is staggered so that mass and other scalars such as temperature and the velocity components are interlaced. The vertical grid is a terrain-following σ coordinate (Greve and Blatter 2009, chapter 5). For our time frame of interest, the most salient atmospheric forcing fluctuation is near the upper surface. The stain rate, however, usually is large near the bottom. Instead of uniformly increasing the model layer number to better represent vertical simulation domain, SEGMENT, to be efficient, uses a vertically stretched grid stencil that better differentiates both bottom and near-surface ice domain (Figure 1b). In the illustration, we use 31 vertically stretched layers in delineating the sliding material thickness, including one lithosphere layer at the bottom. The vertical layers are finer near the upper and lower boundaries but coarser in between, according to a hyperbolic tangent stretching. As a result, the uppermost layer is usually less than 5 cm thick, fine enough to accurately simulate the upper surface moisture and energy states on hourly time scale. Right now, the modeling frame allows arbitrary refining of the vertical layer division, with vertical stretching being the default setting. The u- and υ-component velocities are solved for explicitly, and the w component is obtained from the continuity equation [Equation (A.1)]. It is noted that, even for swelling clay, Equation (A.1) is applicable. The volume of swelling clays changes with soil moisture content, so the vertical regridding takes the expansion ratio into consideration in the time marching. To effectively reduce noise in the flow caused by the model numerics and machine accuracy, an Asselin (Asselin 1972) time filter is applied with a coefficient of 0.15.
The horizontal grid spacing is the digital elevation map (DEM) resolution. The SEGMENT required static input data are of very different horizontal resolutions. They are linearly downscaled to the DEM resolution before being ingested. In addition, the SEGMENT has an adjoint-based (e.g., Ren 2005; Ren 2010) four-dimensional variational data assimilation system. Except for the adjoint code, which is designed strictly corresponding to the forward landslide model component, all components such as optimization and adjoint code verification are taken from Ren (Ren 2004).
A.2. Seismic considerations
Earthquakes are well known to be difficult to predict, so it is prudent to use prescribed earthquake conditions to determine the stability of a slope. SEGMENT incorporates lateral and vertical accelerations of a dynamic nature generated by an assumed imaginary earthquake. At the time of writing, these include changes in soil strength properties, pore pressure changes, and alteration of the macropores. A static seismic coefficient is used to assign extra basal velocity and loading to the sliding material for the period of the earthquake.
A.3. Time stepping within the landslide model
When using the model in real landslides, high-resolution real elevation and soil vertical datasets are used. The time step for integration is 5 × 10−4 s (δts) for soil/rock mechanics and half an hour (δtl) for the soil hydrological component. There is a critical soil moisture lookup table using stand-alone runs of the landslide model for a range of slope angle, slope length, soil thickness, granular size distribution, and bedrock material strength. When the flow speed reaches 0.005 mm day−1, the soil moisture value (critical soil moisture value) is recorded for the combination of the factors in a static array. The land surface model is run over the entire region and the soil moisture of each grid is compared with its critical soil moisture value. If any grid point is unstable, the mudslide model is triggered on for the entire region and integrated for a time span of δtl using integration time step δts. The surface elevation is updated, and so are the soil thickness profiles. If flow speeds are significant (>0.005 mm day−1) somewhere, the flow velocity field for the entire region is recorded for the last three integration steps (measured in δts) for a restart integration following the next land surface model integration. The land surface model integration then continues, with the updated topography and soil thickness profiles and evaluates the slope stability. If there are unstable slopes or the previous velocity fields are significant, the mudslide model will be integrated, initialized with the recorded three steps of velocity fields, forming a complete model cycle.
Corresponding author address: Dr. Diandong Ren, Department of Geological Sciences, The University of Texas at Austin, Austin, TX 78712. email@example.com