Changes in storm-triggered landslide activity for Southern California in a future warming climate are estimated using an advanced, fully three-dimensional, process-based landslide model, the Scalable and Extensible Geofluid Modeling System for landslides (SEGMENT-Landslide). SEGMENT-Landslide is driven by extreme rainfall projections from the Geophysical Fluid Dynamics Laboratory High Resolution Atmospheric Model (GFDL-HIRAM). Landslide changes are derived from GFDL-HIRAM forcing for two periods: 1) the twentieth century (CNTRL) and 2) the twenty-first century under the moderate Intergovernmental Panel on Climate Change Special Report on Emissions Scenarios A1B enhanced greenhouse gas emissions scenario (EGHG). Here, differences are calculated in landslide frequency and magnitude between the CNTRL and EGHG projections; kernel density estimation (KDE) is used to determine differences in projected landslide locations. This study also reveals that extreme precipitation events in Southern California are strongly correlated with several climate drivers and that GFDL-HIRAM simulates well the southern (relative to Aleutian synoptic systems) storm tracks in El Niño years and the rare (~27-yr recurrence period) hurricane-landfalling events. GFDL-HIRAM therefore can provide satisfactory projections of the geographical distribution, seasonal cycle, and interannual variability of future extreme precipitation events (>50 mm) that have possible landslide consequences for Southern California. Although relatively infrequent, extreme precipitation events contribute most of the annual total precipitation in Southern California. Two findings of this study have major implications for Southern California. First is a possible increase in landslide frequency and areal distribution during the twenty-first century. Second, the KDE reveals three clusters in both the CNTRL and EGHG model mean scarp positions, with a future eastward (inland) shift of ~0.5° and a northward shift of ~1°. These findings suggest that previously stable areas might become susceptible to storm-triggered landslides in the twenty-first century.
Trends in extreme rainfall over the twenty-first century in mountainous areas have great relevance for storm-triggered landslides, particularly when the slopes are precipitous and prone to slippage following heavy rain events. For example, Southern California experiences a range of landslide types, including shallow, fast-moving landslides triggered by rainstorms on the steep, unstable slopes of the region. Almost all of Southern California’s rain falls in winter (November–April) and the most common Southern California landslides that are triggered by winter rainstorms are debris flows, which often are referred to as mudslides. Debris flows usually occur soon after a heavy rainstorm event, and during wet winters with seasonal rainfall totals exceeding 250 mm. Such conditions are most likely in southern California during the warm (El Niño) phase of El Niño–Southern Oscillation (ENSO) (see, e.g., Neelin et al. 1998; McPhaden et al. 2006; Luo et al. 2011; Bretherton et al. 2004).
To explore future changes of rainfall over Southern California and their implications for landslide trends in the twenty-first century, high-resolution general circulation models can provide simulations of future rainfall (e.g., Knutson et al. 2010) Here we use the Geophysical Fluid Dynamics Laboratory High Resolution Atmospheric Model (GFDL-HIRAM; Putman and Lin 2007). GFDL-HIRAM is verified against observed precipitation over Southern California to validate whether the physical processes in GFDL-HIRAM are the same as those known to control rainfall changes on interannual, decadal, or longer time scales. If GFDL-HIRAM reproduces the major climate drivers of precipitation over Southern California, the credibility of the GFDL-HIRAM twenty-first-century rainfall projections is increased (Palmer et al. 2008). The Scalable and Extensible Geofluid Modeling System for landslides (SEGMENT-Landslide; Ren et al. 2007, 2009, 2012), driven by atmospheric fields from GFDL-HIRAM, then can justifiably be used to extend our projection ability of future landslides changes in California, especially Southern California, through the twenty-first century. As described in detail below, SEGMENT-Landslide is one of few fully 3D, process-based landslide models available and has been used in a range of landslide studies.
The study region (Fig. 1, ~390 000 km2), which includes Southern California, comprises several low mountain ranges in the southwest corner of California, near Los Angeles and San Diego, and the higher, sparsely inhabited Sierra Nevada range to the east. The complex geological properties of the region support essentially all type of landslides, as defined in Hutchinson (1968) and Varnes (1978). For example, when the coastal bluffs (of Cenozoic marine sedimentary rocks) experience storms in El Niño years, deep-seated and rotational landslides are likely to occur. For the scars left by wild fires or after timber harvesting, shallow landslides such as the progressive bulking type of debris flows are common for intense rainfalls. Land cover and land use are a mix of urban, shrub, and forest environments. The dominant vegetation types of Southern California are winter deciduous and evergreen trees and shrubs, annual grasslands, and some drought-resistant deciduous plants. Material properties vary widely, as do local surface relief, vegetation, precipitation climatology, and seismic activity. These factors determine where and when landslides will occur. Previous landslides leave a “cold trace legacy” (Dietrich and Dunne 1978; Reneau and Dietrich 1987; Dietrich et al. 1995; Crozier and Preston 1999) on future slope stability (i.e., reducing the possibility of another similar occurrence). Detailed process modeling is needed to identify future unstable regions.
In this study, SEGMENT-Landslide provides twenty-first-century projections of trends in storm-triggered landslides over California. This study region is data rich and is used to demonstrate that the GFDL-HIRAM model used can reproduce most climate aspects of rainfall over Southern California. Moreover, SEGMENT-Landslide is capable of investigating all major landslide types that occur in Southern California. Section 2 provides details of SEGMENT-Landslide, and its application methodology. Section 3 describes the available data for both input to the model and for verification of the model simulations. Section 4 presents the results of the model projections, and section 5 discusses the results and draws conclusions arising from the findings of the study.
2. The SEGMENT-Landslide model and methodology
A range of techniques is employed to analyze Southern California rainfall, with emphasis on trends in extreme precipitation, as landslide activity is defined here as that triggered by extreme rainfall, with a landslide mass of 103–106 m3. Statistical methods for the rainfall estimations include Monte Carlo–based composite analyses (Hendon et al. 2007), scattering analysis (Scott 1992), and nonparametric trend analysis (e.g., Mann–Kendall test; Kendall 1975). These statistical techniques are described briefly in appendix A.
a. Landslide models
Present and future landslide research is concerned with generating more accurate predictions about landslide timing, location, and size, and about how the sliding material is redistributed. There are three categories of landslide models available at present: rainfall intensity–duration (ID) threshold line models (Caine 1980; Campbell, 1974, 1975); 2D slope stability models [specific examples include Transient Rainfall Infiltration and Grid-based Regional Slope-Stability (TRIGRS; Baum et al. 2008), the Shallow Slope Stability Model (SHALSTAB; Montgomery and Dietrich 1994), and stability index mapping (SINMAP; Pack and Tarboton 2004)] and full 3D, process-based landslide models (e.g., SEGMENT-Landslide; Ren et al. 2007).
In ID methods, a dataset of rainfall intensity (I; mm day−1) and rainfall duration (D; h) of landslide events is prepared. A scatterplot then is generated with rainfall intensity plotted against rainfall duration. The rainfall threshold equation is a power-law curve that fits the points in the scatterplot (it is a lower envelope, because if a point lies to the upper right of the curve, landslides are possible) and usually is of the form , where a and b are positive constants of soil, vegetation, and land use. Importantly, ID thresholds are only one of the many kinds of rainfall thresholds proposed in the literature. Some researchers have demonstrated that deep-seated landslides are triggered by prolonged rainfall, so thresholds based on total rainfall amounts are usually preferred to ID thresholds (Lagomarsino et al. 2013).
The slope stability models (SSMs), like the IDs, are concerned only with where landslides will occur. Slope stability models are based on the following reasoning. On a sloping surface, the gravitational force can be partitioned into a component normal to the slope (Fn), contributing to friction that resists sliding erosion, and a component parallel to the slope (Fp) that promotes sliding. A stability parameter S is defined as , where is the internal friction coefficient. Note that S is also referred to as the stability factor. This is the general form of S, but there are many specific forms, based on the fact that landslides, as downslope movements of rock and debris or earth, occur when the shear stress exceeds the shear strength (i.e., when S < 1).
To systematically model a landslide life cycle, a complete landslide model must solve the governing equations of mass, momentum, and energy of sliding material for a specific granular rheology. The following description therefore focuses on 3D, physically based process modeling, using SEGMENT-Landslide as an example. It can be shown that both IDs and SSMs are reduced forms of the full 3D model, with different emphases on the failure criteria.
For the sliding material, a coupled system is solved for conservation of mass,
under the multiphase rheological relationships, for water, condensed mud, and wet sliding granular materials, with granular viscosity parameterized as
where V is velocity, with horizontal components u and υ and vertical component w (as in atmospheric sciences), h is surface elevation of slope, F is external force (e.g., gravity), is the full second-order internal stress tensor, and indicates evaluation at the free surface elevation. In the case of slope movements, Eq. (1) is solved regularly to update the sliding material geometry. The w terms may include sedimentation rate and entrainment rate for clear debris flows. After the sliding flow entrains the material, it changes the rheological property of the medium and is advected within the sliding flow. Here is bulk density; V is velocity vector; F is the body force (e.g., gravity ); υ is viscosity; is the spherical part of the stress tensor (R is the resistive tensor); μ0 and μ1 are the limiting values for the friction coefficient μ; is the effective strain rate defined as ; I0 is a constant depending on the local slope of the footing bed as well as the material properties; and I is the inertial number defined as , where d is particle diameter and is particle density. A soil moisture factor enhancing viscosity is included and is assumed to vary according to a sigmoid curve, expressed formally as Eq. (9) of Sidle (1992), but with the time decay term replaced by relative saturation. The transient soil moisture regime in SEGMENT-Landslide is obtained by solving the 3D Richards equation, using a land surface model ( appendix B).
The viscous term in Eq. (2) implies an energy conversion from kinetic energy to heat. To ensure full energy closure, the following thermal equation is required:
where c is heat capacity (J kg−1 K−1), T is temperature (K), k is thermal conductivity (W K−1 m−1), and σeff is effective stress (Pa). The last term is “strain heating,” which is the converting of work done by gravity into heat and affects the sliding material by altering its viscosity, or by causing a phase change. The above equations form the landslide component of SEGMENT-Landslide. Other components are shown in Fig. 2. That 2D slope stability models are simplified forms of this modeling system is easy to demonstrate. Omitting the time dependence (because SSMs do not include the sliding material’s ensuing motion and mass redistribution) and working only in the x–z plane for simplicity, the governing equations can be written in component form as
Further assuming that the slope geometry has a constant bedrock slope (α), a uniform sliding material thickness (which implies that the surface slope also is α), and a stable bed, the volume integration of Eqs. (5) and (6) yields
where is the hydrostatic pressure from groundwater, is the net vertical support at the top and the toe (determined by boundary conditions), is the horizontal support evaluated at the top and the bottom (or “toe”), and is the horizontal force exerted during an earthquake. Equation (7) indicates that slope stability is influenced by various factors, such as the slope gradient (), soil properties (implicit in G, Ceff, and ), groundwater table, and geomorphology ( and the boundary conditions). The triggering mechanisms therefore formally include not only the excessive rainfall, which is the focus of this study, but also possible earthquake and volcanic activities (). If the full stress is defined generically, then Eq. (7) is the Coulomb–Navier fracture criterion for internal stress [RHS of Eq. (8)]:
where is the gravitational potential, C is effective cohesion, also called apparent cohesion or interlocking cohesion (i.e., the shear stress that the material can sustain at zero normal stress), and is the pressure perturbation caused by earthquake-, volcano-, or human-induced disturbances at that location. Note that μ < μ1 = tgΦ, where Φ is the granular repose angle (the angle of friction). It is found experimentally that the behavior of soil materials above the water table can be adequately represented by the Coulomb–Navier fracture criterion. Hydraulic properties and the root reinforcements (detailed below) all are included in the apparent cohesion C in Eq. (8). As soil moisture plays a critical role in slope stability, the SSMs often have a prognostic soil moisture profile in the sliding material. However, unlike the full 3D landslide models, the SSMs typically use analytical solutions to the Richards equation (see appendix B). For example, with some extension, TRIGRS uses the analytical solution described in Srivastava and Yeh (1991) to link with a stability factor, enabling an assessment of landslide potential over broad geographical areas.
For landslides of different origins, the factors in Eq. (7) contribute disproportionally. For storm-triggered landslides, the focus of this study, soil moisture changes are the leading cause. The surface runoff has a positive feedback in the magnitude of storm-triggered landslides. Rainfall morphology also is critical in triggering landslides (Ren et al. 2010, 2011). The ICs, as first-order approximations, only address the vital role played by precipitation, which is the source of soil moisture, and use a rainfall intensity–duration power-law relationship, ascribing all geological, climatological, and geographical controls to just two empirical parameters.
Soil slips, debris flows, shallow landslides, and slumps (Varnes 1978) are seemingly different in their sliding material collection region, the abruptness of the movements, and the resultant flow and redistribution of the sliding material. In SEGMENT-Landslide, neighboring grids transfer information through the conservation equations of mass, momentum, and energy. Provided the simulation area is sufficiently large, several forms of landslides can coexist and be accurately modeled by SEGMENT-Landslide. In contrast, slope stability models are more suitable for the case where sliding material is locally scraped out of the residing slope, whereas full 3D models can treat cases where the sliding material is transported from elsewhere. For example, in “progressive bulking” debris flows there usually are large upstream collecting areas, where infiltration–excess runoff entrains debris and transports it to creeks/valleys downstream. Then the accumulated material slides down as visible debris flows. For the upstream collecting area, the slopes usually are stable (safety factors S > 1). The progressive bulking debris flows have been shown to be tractable in SEGMENT-Landslide (Ren et al. 2013).
As a fully 3D model, SEGMENT-Landslide uses a complete form of energy bundle (i.e., the grouping energy terms; e.g., the moist static energy in atmospheric science is cpT + lq + gz, which is a form of energy bundle), allowing investigation of a wide spectrum of landslides, from shallow, storm-triggered landslides and debris flows to deep-seated rotational landslides. The mechanical, chemical, and hydrological effects of vegetation roots are implemented. For example, a scheme is proposed in this study for implementing the reinforcement of soil strength from live vegetation roots. The uniqueness of this scheme is that remotely sensed biomass is used to infer root mechanical properties. The governing equations for mass, momentum, and energy are documented in Ren et al. (2008, 2012) and need not be repeated here. Root mechanics is implemented because of its relevance for the slope stability of the region of interest. In many of these slopes, several kilopascals of stress can determine the shallow stability of materials.
The role of vegetation in increasing slope stability is well known (Wu 1995; Coder 2009; Morgan and Rickson 1995, and references therein; Wu and Watson 1998) and shear tests have been carried out in the laboratory (Gray and Ohashi 1983) and in the field (Abe and Iwamoto 1988) and have even been retrieved from failed slopes containing a root membrane (Riestenberg 1994). The presence of above-ground vegetation introduces the following effects: above-ground biomass loading (gravitational), growing season soil moisture extraction by live roots (hydrological), fortification of the soil within its extension range (mechanical; Fig. 3a), and changing chemical environment of the soils through life processes (e.g., respiration, absorption of minerals selectively, and secretion of organic substances) and therefore the bond strength among unit cells (chemical), and wind stress loading (meteorological). The overall effect is from the interaction of the above factors, and it is difficult to generalize before a detailed analysis is carried out that is specific to a certain situation. For example, the fortifying roots have yield strengths larger than dry soils and the existence of roots is commonly thought to increase the resistance of soils. However, the presence of roots, especially when there is precipitation, also facilitates water channeling into deeper depths. After the soil is moistened, the cohesion between soil and the root surface is reduced greatly, to negligible strengths < 0.001 MPa (Lawrance et al. 1996) and the root strength cannot be effectively exerted (i.e., tensile strength is only “partially mobilized”). Also, the effect of roots is to “unify” the soil particles within root distribution range (Fig. 3a). Once the entire rhyzosphere soil layer, which is that part of the soil layer affected by the plant roots, is saturated, the fortifying effects are lost. Thus, the reinforcement effects from roots are an effective mitigating factor for shallow storm-triggered debris flows.
It is known that shallow interlocking root networks can contribute to mechanical reinforcement to soils (Sidle et al. 1985; Selby 1993; Lawrance et al. 1996). For a pasture species, Selby (1993) estimated the “additional” cohesion as ranging from 0.1 to 9.8 kPa, with changes in soil moisture. There also are experimental tests of root strength for a variety of species (Table 1). The parameters obtained from in situ sampling usually are inadequate as input for 3D mechanistic models that run on regional scales.
Coder (2009) tabulated root strength for a variety of plants. In addition, roots from the same plant can behave differently: they behave elastically for thinner roots but elastoplastically for thicker roots (Pollen and Simon 2005). The maximum strength of roots depends primarily on the fiber content in the xylem. Thus, ephemeral absorbing roots have different strengths than structural roots. The same types of roots of different thickness have similar tensile strengths per unit cross-sectional area (Table 1). For example, the structural roots of the Swiss stone pine of thickness 12.6 and 13.9 mm have similar tensile strengths of 12.9 and 13.5 MPa, respectively. Because different vegetation, or the same vegetation in different growing stages or living in different climate zones, can have different root mass densities or even orientations within their ranges of influence in the soil, generic schemes are still lacking for estimating root strengths for regional or global scales. Because roots occupy a small fraction of the soil volume, the root reinforcement can be factored in as an added stress over the case of no root presence (Fig. 3). As not only tensile strength but also the cross-sectional area (thickness of roots) is critical, the following “allometric” approach is proposed that uses “root weight density” and vegetation type to characterize the added tensile strength to the soil medium shear strength and elevated yielding criterion:
where is the root mechanical reinforcement [kPa, as an additional term within C of Eq. (8)]; Csmc is the soil moisture control (0–1); index i differentiates woody transport roots (i = 1), woody supporting roots (i = 2), and ephemeral absorbing roots (i = 3); (kPa) is a species-dependent reference value (i.e., the tensile strength of xylem of roots varies from 10 to 30 MPa for most plants); NPP is net primary production (kg m−2 yr−1); Pa is annual precipitation (mm); Tair is annual mean air temperature (K); and Iveg is the species-dependent reference NPP value. Also, F varies from 0 to 1 and represents the weight fraction of roots in a unit volume of soil, within the range of influence of the roots. The value of F typically is close to 0.2%. The functional form of Csmc can be tabulated according to soil type. The modifications resulting from climate conditions are necessary because the same plants have very different strategies in allocating biomass when they are located in different climate zones. For example, plant leaves tend to be thin and wide in moist tropical rain forests. The same plant, growing in drier locations, tends to have smaller leaves and a more extensive root system to better use the limited water resources. A positive side of this formula is that root volumetric (or mass) density can be obtained remotely using satellite data (Zhao and Running 2010). Compared with existing upscaling schemes for root-introduced extra cohesion (e.g., Pollen and Simon 2005), this scheme emphasizes that roots with similar ecosystem functions have similar tensile strength, rather than simply assuming that roots of similar thickness have similar strength.
Mechanical properties of fractured bedrock are quantitatively formulated. The applications of this model to several recent cases in various climate zones and different triggering mechanisms were successful. Quantifying the damage capacity of storm-triggered landslides is critical for evaluating their deforestation effects. To this end, SEGMENT-Landslide has a sound speed–based index. In SEGMENT-Landslide, a metric for estimating deforestation effects is presented in terms of ecosystem parameters (Fig. 2). High-resolution simulations of the region require very large in-memory computer processing. The SEGMENT-Landslide code is fully parallelized and all simulations in this study were carried out on the iVEC supercomputer at Curtin University, Perth, Australia. The lakes, open waters, and high-elevation glaciated parts of the Sierra Nevada are automatically filtered out, based on invalid land cover types.
The model simulation domain contains all possible sliding material (the entire soil layer) and one bedrock layer as the bottom model sublayer. SEGMENT-Landslide requires the following input parameters: initial moisture content of the failed slab, soil material properties, vegetation root depth, and distribution morphology such as root length density and root strength and initial bedrock topography (high-resolution elevation data). High-resolution elevation data are needed to identify unstable slopes because it is the surface elevation gradients that drive the downward motion. This study uses the high-resolution digital elevation map (DEM) from the National Aeronautics and Space Administration Shuttle Radar Topography Mission (SRTM; http://www2.jpl.nasa.gov/srtm/). The dataset is ~90-m resolution. Higher-resolution DEMs are available, but for this study 90 m was chosen as it enabled the computations to be made in reasonable time without compromising the aims of this study. The soil characteristics are from the Food and Agriculture Organization (FAO) of the United Nations (http://www.fao.org/geonetwork/srv/en/metadata.show?id=14116). These soil mechanical properties (Table 2) are used to parameterize the granular viscosity of the sliding material according to Eq. (3). The setting of the bedrock mechanical properties references to the geological map of California. At the fracture stage (velocity is zero in the governing equations), the static mechanical strength is set at 30 MPa for the intrusive igneous rocks and other crystalline rocks of the Sierra Nevada range and only 5.6 MPa for the sedimentary rocks at the coasts (e.g., shale and clay stones). Weakly cemented sandstones have an intermediate value of 11.1 MPa. The compressive strengths of rocks are orders of magnitude larger than their tensile strengths and rarely are involved in the landsliding events. At flow stages (with apparent velocities of sliding material), the viscosities of the bedrocks also vary from 1015 to 109 Pa s, depending on lithology and degree of fracture. In setting up the bedrock mechanical strength, the burying depth, in addition to the geologic unit reflecting lithology, is taken into consideration to reflect increasing rock strength with confining pressure.
Compared with the physical properties of hill slopes that affect the potential for landslides, vegetation and rainfall are far more temporally variable. The U.S. Geological Survey National Land Cover Dataset 2001 (http://www.epa.gov/mrlc/nlcd-2001.html; Homer et al. 2004) provides land cover information (of reliable kilometer spatial resolution). For vegetated surfaces, especially those affected by fires, representing soils as purely mineral (e.g., sand, silt, and clay classifications) is inadequate (DeBano 2000; Moody and Martin 2001; Smith et al. 2002). For previous fire sites within the last decade, a parameter fi is introduced to represent the fraction of soil sublayer (i) that is organic matter. The thermal and hydraulic parameters for organic soil are weighted averages of the corresponding mineral soil and those of humus, using an empirical relationship (Lawrence and Slater 2008; Letts et al. 2000). A unique feature of SEGMENT-Landslide is its treatment of root properties relevant for storm-triggered landslides. Bedrock rheological properties such as viscosity are set according to published geological maps (http://mrdata.usgs.gov/sgmc/ca.html). The high-resolution climate model GFDL-HIRAM provides the atmospheric parameters.
As a result of the complex geological properties underneath Californian slopes, soil type may vary within the same slope (Segoni et al. 2012; Mahmood et al. 2013). For storm-triggered landslides the focus is on the mechanical properties of the soil particles. Many rheological features and hydromechanical properties of granular material originate from particle interlocking, which relies on soil particle size distribution. Ideally, the soil classification should classify soils such that their hydromechanical properties are similar within the same classification. A number of soil classification systems are available, depending on the disciplines and their general engineering or scientific objectives. Some of the most widely used classification systems in the United States are the U.S. Department of Agriculture (USDA) system (Eswaran et al. 2002), the international system (Buol et al. 2003), the Unified Soil Classification System (ASTM 1985), and American Association of State Highways and Transportation Officials system (Hogentogler and Terzaghi 1929). The commonality among these systems lies in the use of soil particle sizes in dividing soils into different categories of clay, silt, sand, gravel, and larger materials, although the boundaries for each of these systems are somewhat different (see Soil Survey Division Staff 1993 for more details). Here, the USDA classification system is used because SEGMENT-Landslide’s hydrological submodel component is parameterized according to this classification system. Soil chemical composition also is important for hydromechanical properties. Unfortunately, none of the four classifications is based on soil chemical compositions. In setting the spatial variable soil mechanical properties, the State Soil Geographic (STATSGO; http://water.usgs.gov/GIS/metadata/usgswrd/XML/ussoils.xml) soil data were used, as well as the FAO soil maps (reclassified using USDA soil types over the region of interest), for its wealth of vertical profile information. The general range of soil mechanical properties is listed in Table 2.
To study the daily variation of California rainfall associated with climate drivers such as Pacific–North American pattern (PNA), Pacific decadal oscillation, North Atlantic Oscillation (NAO), and Southern Oscillation index (SOI), the North American Regional Reanalysis daily composite precipitation is used for the period 1979–2010 (http://www.esrl.noaa.gov/psd/data/gridded/data.narr.html). The indices are obtained from the National Oceanic and Atmospheric Administration National Climatic Data Center website (http://www.ncdc.noaa.gov/teleconnections/) For the projection of future global warming consequences on storm-triggered landslides, the GFDL-HIRAM, specifically the C180 HIRAM2.1 simulations, are chosen over other Program for Climate Model Diagnosis and Intercomparison climate model simulations for AR4 primarily because of its finer horizontal (~50 km) and vertical grid spacing (32 vertical levels) in the version described by Zhao and Held (2010). The twentieth-century run is referred to as the control run (CNTRL). For the twenty-first century projection, only the Special Report on Emissions Scenarios (SRES) A1B emission scenario runs are used. The simulations used here are for two continuous 20-yr periods: from 1979 to 2000 and 2079–99. The SRES A1B scenario run is referred to as the enhanced greenhouse gas run (EGHG).
4. Results and discussion
The primary sources of precipitation over the west coast of the United States are the synoptic troughs traceable to parent Aleutian lows. The GFDL-HIRAM provides very satisfactory simulations of the occurrence and tracks of these precipitation events (not shown). However, the precipitation amount and intensity are found to be underestimated when compared with observations and there usually is little damage, because the troughs weaken or diminish as they pass over the colder waters of the California Current flowing southward along the California coast. The longer-lasting rainfall events are usually associated with ENSO events, which have about a 4–7-yr recurrence frequency (McPhaden et al. 2006). These rainfall events (usually >50-mm rainfall) cause most of the observed debris flows in Southern California (Ren et al. 2011). Closely related to ENSO events, GFDL-HIRAM occasionally produces a southern, subtropical branch of storm tracks in addition to that originating from the Aleutians (not shown). This is a response to the unusually warm sea surface temperature (SST) in the eastern tropical and subtropical Pacific, and the storm tracks bring tropical moisture to Southern California. In this study of the GFDL-HIRAM output it is shown that the most severe rainfall events over Southern California result from midlatitude troughs trapping the moisture remnants of hurricanes that formed in the eastern North Pacific Ocean and traveled sufficiently far northwest. This situation is found to be rare (~27-yr recurrence period) but can produce rainfall totals exceeding 200 mm within several days, consistent with the studies by Blake (1935), Smith (1986), and Landsea and Chenoweth (2004). The climatological statistics of hurricanes simulated in GFDL-HIRAM, as well as the interannual variability in the model forced by observed SSTs, indicate that it can reasonably simulate the geographical distribution of storm genesis locations, as well as the seasonal cycles and interannual variability of hurricane frequencies for the major basins (Zhao et al. 2009). The deficiency of GFDL-HIRAM in simulating the >100-mm category rainfall events over Southern California is due primarily to an inability in simulating tropical cyclone intensity (Zhao and Held 2010). Thus, GFDL-HIRAM simulates well the occurrence of this type of event but the simulated intensity and longevity were found to be inadequate at present. Rescaling the model precipitation empirically (e.g., Zhao et al. 2009) can better fit the observed rainfall morphology. In this study, we do not make such an empirical correction because this deficiency should not be allowed to influence GFDL-HIRAM’s ability to indicate changes of extreme rainfall events as the climate warms. That is, GFDL-HIRAM projections are highly valuable since they indicate whether there are significant changes in occurrence frequencies of extreme rainfalls. Instead, we reclassify the precipitation bins as 0–5–20–50–100 mm and then the correspondence between the CNTRL and reality is improved significantly (figure not shown).
The mountainous terrain south of the Peninsular Range and Transverse Ranges also receives heavy rain from the North American monsoon. The GFDL-HIRAM provides satisfactory timing and location of this monsoonal precipitation. Note that the coupled models cannot produce precise timing in their simulations, but do provide statistically acceptable time series to assess the probability of changes in extreme rain events with landslide consequences.
To investigate the relationship between the daily precipitation variations and the teleconnection patterns, composite maps for the period 1980–2010 were constructed for all days in high index polarity and low index polarity (i.e., periods with indices beyond one standard deviation). Composite differences are defined as the composite high index minus the composite low index. The natural variability of composite differences or, if the index is completely unrelated to the precipitation variations, is estimated from a Monte Carlo method (Hendon et al. 2007). For the SOI except for spotty areas north of California’s Central Valley, all were found to be significant at the 95% level. The correlations depend on geographical location and can explain up to 17% of the total variance (e.g., for the region containing Yosemite National Park, Mandeville Canyon, and Oroville). Similar analyses carried out in this study for the Arctic Oscillation and NAO indices show only modest correlations. From the composite analyses it is apparent that the most notable precipitation signals on an annual time scale across California are from El Niño events, which show an overall positive anomaly for all seven subclimate zones, especially Southern California. The composite analysis adapted from Hendon et al. (2007) reveals that as much as 18% of the variance can be explained by the SOI. On the longer decadal time scale, the fluctuations in precipitation are related to the PNA patterns. For Northern California, it can explain as much as 20% of the annual variance. Thus, whether GFDL-HIRAM can successfully simulate the El Niño–related precipitation is critical in estimating changes in storm-triggered landslides.
Figure 4 shows that, characteristic of a Mediterranean climate, the regional average rainfall for Southern California is concentrated in the 5-month period between November and the following March. Although the annual mean precipitation is accurately simulated, the error in February precipitation can reach 18%. The light precipitation (<25 mm), although occurring most frequently (~70% of rainfall occurrence), comprised less than 40% of the annual total rainfall. However, the extreme rainfall (100–250-mm category) occurred 20 times less frequently but contributed about 25% of the total precipitation. The model significantly underestimates this rainfall category.
Unlike typical precipitation events with lifespans of less than a day, extreme precipitation, with total amounts greater than 100 mm in the southern California region, usually are distributed over consecutive days, or several nonconsecutive days separated by no more than two days. For storm-triggered landslides, those precipitation events separated by less than two dry days (a “break”) are considered as one “super rain” event. Thus, unlike many previous studies (e.g., O’Gorman and Schneider 2009, and references therein), which count daily precipitation one day at a time, this study counts those extended super rain events, defined as a near-continuous rainfall period nowhere separated by more than two consecutive dry days (Ren et al. 2011). Figure 5 shows changes in rainfall morphology as the climate warms. Rainfall categories are set according to the extreme rainfall events defined in Ren et al. (2011). For the control period (x axis), except for the <5-mm category, all other categories have apparent (2–6 times) errors of persistent underestimation of rainfall occurrence frequency (climate model produces rainfall less frequently than reality). The GFDL-HIRAM deficiency in simulating the extreme categories (>50 mm) arises mostly from the weaker simulated hurricanes formed in the eastern North Pacific Ocean. The deficiency can be corrected by reclassification or by proportionally increasing the modeled precipitation. For example, suppose the 5–10-mm category corresponds to the 50–100-mm observed category and so on. The large rainfall categories (50–100- and 100–250-mm categories) have a tendency to increase in occurrence frequency (the figure lying above the diagonal line). By examining three-dimensional atmospheric flow fields, it is noteworthy that there is a broadening of the tropical meridional circulation.
It is of interest to note that in this study, to identify a tropical cyclone (TC)-like vortex in GFDL-HIRAM, a number of criteria were used. These criteria are required because the resolution of climate models still is too coarse to simulate TC vortices in detail. The criteria used here are similar to those employed by Oouchi et al. (2006), designed specifically to extract TCs with genesis in the tropics and to avoid the extraction of midlatitude cyclones. They are as follows:
The 850-hPa relative vorticity exceeds a threshold value of 3 × 10−5 s−1. In the Southern Hemisphere a negative threshold of the same magnitude must be used.
The maximum near-surface (10 m) wind speed in a centered 7° × 7° box exceeds the wind speed threshold of 10 m s−1.
The sea level pressure is the minimum in a centered 7° × 7° box.
The temperature anomaly averaged over the 7° × 7° box and three pressure levels (300, 500, and 700 hPa) exceeds the temperature anomaly threshold.
The local temperature anomaly averaged over the 7° × 7° box is positive at all three pressure levels (300, 500, and 700 hPa).
The local temperature anomaly, averaged over the 7° × 7° at 300 hPa is greater than at 850 hPa.
The mean speed averaged over a 7° × 7° grid box is larger at 850 hPa than at 300 hPa.
The grid points representing the center of storms that obeyed all the criteria above are connected if they are less than a certain distance from the center of the previous time step analyzed. This distance is defined by the frequency of the output of the model. For 6-hourly outputs, two grid points are used in longitude and/or latitude, while for daily output three grid points is the maximum distance.
The storm lasts for at least 36 h.
Using a hurricane-tracking program (Walsh 2004), the landfalling tropical cyclone numbers are found to increase significantly, reaching ~20-yr recurrence period. A clustering analysis ( appendix A) also indicates that the landfalling locations of tropical cyclones move poleward during the twenty-first century. This factor alone can explain over 50% of the extreme rainfall increase.
From the GFDL-HIRAM atmospheric fields for the twenty-first century, SEGMENT-Landslide can estimate changes in landslide occurrence over Southern California (Fig. 6). Because single landslides cover only a very limited area, plotting the individual scarp is difficult. The markers of different sizes in Fig. 6 indicate the average magnitudes of landslides within a grid box of ~2000 km2, defined as the amount of sliding material involved. The model grid is 90 m × 90 m but, because of the limited figure size, a thinning ratio of 50 (i.e., using every 50th model grid value) is applied in each direction. Markers of different types indicate climatological occurrences of storm-triggered landslides, with plus signs specifying 0.05 yr−1, circles indicating 0.1 yr−1, and barred circles representing 0.25 yr−1. Contour lines are increases of occurrence frequencies relative to climatological occurrence frequency. There are two “hot spots” along the coastal area: the La Conchita area (~34.36°N, 119.4°W) to the south, and the Oakland area (~37.8°N, 122.3°W) to the north. Note that these two hot spots also are populated areas (Fig. 1). These two areas have the most frequent occurrence of mudslides but the magnitudes are moderate, as <105 m3 of sliding material is involved for a single event. Precipitation over the Sierra Nevada, especially the western slope, increases significantly during the twenty-first century. However, because the base rocks are primarily least weathered–fragmented granites, the occurrence of debris flow is infrequent. However, for scattered areas covered by volcanic materials, such as the northernmost cross marker in Fig. 6, storm-triggered landslides are possible. The increase in occurrence frequency is significant (reaching once per 20 yr at some locations), relative to the base occurrence frequency. For these less populated regions the landslide magnitudes, which primarily are from deposits of rockslides, generally are larger than those occurring along the coast. The hydrological processes, the geotechnical properties of the sliding material, and even the thickness of the sliding material control the size of a landslide. Malamud et al. (2004) proposed an inverse gamma distribution to describe the size distribution for the landslide events. The three-parameter formula [their Eq. (3)] seems to be globally applicable. As the climate warms, the parameter , which controls the power-law decay of occurrence frequency with increasing landslide areas, tends to decrease. That means the landslide size spectrum has a shift toward the large landslide categories. For our study region, decreases from 1.7 to 1.5 for 2079–2100 compared with the 1979–2000. Potential landslide locations with areas >500 m2 increased by 4%.
The kernel density estimation (KDE) values (see appendix A) were calculated for the CNTRL run landslide locations and the EGHG run landslide locations, for the normal reference 2D bandwidths of (3.95, 1.3) as defined in Eq. (A5). The KDEs suggest that three clusters are present, one in the Sierra Nevada, a second in Southern California centered near Malibu, and a third near Oakland in central California. Note that the three clusters possibly can be interpreted as one large disjoint cluster. A spot is defined as the scattering center of landslide scarps if the sum total of the Euclidean spherical distances between this point and all the scarps as a whole is a minimum. The scatter centers of the CNTRL and the EGHG run periods are compared and it is found that the center is shifted northward by 1° and shifted eastward by 0.5°, from 35°N, 120°W to 35.8°N, 119.5°W, as shown by the stars in Fig. 6. Extreme precipitation events associated with ENSO events and, less important, the northward-shifted hurricane tracks increase landslide numbers in northern coastal California. This also is reflected in the shape of the contour lines in Fig. 6. Except for the limited region parallel to 34.6°N, and the strip portion of the San Fernando Valley, there generally is no landslide-free region (i.e., green line circled areas). More important, very few regions show a decreasing trend of occurrence of landslides.
5. Discussion and conclusions
Landslides in the United States, and particularly in Southern California, are severe hazards that cause substantial human and economic losses, estimated at an annual average of 25–50 deaths and a financial cost of approximately $1–3 billion per year (National Research Council 2008). Simulations of future changes in storm-triggered landslides, also commonly referred to as debris flows, were carried out over Southern California for the period 2079–2100, using projected precipitation from a coupled climate model, the GFDL-HIRAM. Two model projections were assessed. One was a control run (CNTRL); the other was a model run (EGHG) with enhanced greenhouse gas levels, from the moderate Intergovernmental Panel on Climate Change Special Reports on Emission Scenarios A1B future scenario. The results of this study are significant because any significant trends in landslide activity, and/or changes in the landslide genesis locations, have major socioeconomic and environmental consequences for Southern California.
First, the projected trends in occurrence frequencies and intensities were compared. It was found that there are significant differences in the mean numbers of landslides between the CNTRL and EGHG projections, as well as a notable difference in the mean magnitudes, with the EGHG landslides being more frequent and larger, at the 95% statistical significance level. The conclusion is that for Southern California, landslides are likely to become more frequent and more intense during the twenty-first century, relative to the twentieth century. Second, the locations of the landslides from the CNTRL and EGHG projections were compared. The results from a nonparametric permutation test revealed statistically significant differences between the occurrence locations. The importance of this finding is that it demonstrated that previously stable areas could become susceptible to landslides. This partly is caused because our analysis of the GFDL-HIRAM projected an increased numbers of storms along the southern track, as well as more tropical cyclones that bring moisture to higher latitudes, although landfalling tropical cyclones likely will remain rare. Further analysis indicated that the landslide activity (geographical distribution patterns) caused by extreme rainfall events associated with ENSO has not changed significantly. However, the most uncommon category of extreme rainfall, namely that caused by hurricane landfall, was found to have a tendency to extend farther north as the climate warms.
It is emphasized that this study region is very complex and, in particular, is affected by triggering factors other than storms. These include, but are not limited to, wildfire burning of surface vegetation cover, earthquakes, and human activities such as irrigation and construction of hydraulic projects. All are potential contributors to increased landslide occurrences but are not part of this study.
The authors thank Dr. Ming Zhao of GFDL for providing the high-resolution meteorological data for precipitation analysis and for driving the SEGMENT-Landslide. We thank M. Cartina for providing the root strength data and Drs. Matthew Wheeler and Harry Hendon, Bureau of Meteorology, Melbourne, Australia, for providing technical assistance in the composite analyses and for discussions of the theory and application of the nonparametric tests of significance.
The cluster analysis of the landslide positions involved a multivariate kernel density estimation method, with a Gaussian kernel. The univariate KDE approximates the probability density function of a point x as
where n is the number of observations, h is a smoothing variable (the bandwidth), and is a kernel function. There are a number of possible choices of kernel function. For this study the Gaussian (normal) kernel is used, defined as
A two-dimensional KDE can be represented by
where x = [x(1), x(2)], and with separate bandwidths for each dimension. Here, multiple bandwidth pairs (h1, h2) were used for comparison, as well as the bandwidth pair determined by the normal reference rule for a multivariate KDE given in Scott (1992):
where d is the number of dimensions. In this study, d = 2 and thus
Full Form of the Richards Equation
Fluid flow and moisture content in hilly environments vary spatially and temporally due to time-dependent environmental changes (e.g., precipitation, solar radiation, or transpiration of vegetation through distributed roots) and the storage capacity of porous granular soils. In modeling, the time-dependent environmental changes are usually accounted for by boundary conditions of distributed forcing terms, whereas the soil hydrologic property determined soil storage is usually cast into the flow laws. Considering a representative soil volume (RSV), this RSV is of arbitrary volume and arbitrary surface area with an out-normal of n. We assume the soil porosity is P (maximum of the volumetric water content wc) and moisture flux is q.
From Green’s first theorem and mass continuity,
Equation (B1) holds everywhere inside the porous media. From Darcy’s law, soil water flux (kg s−1) is proportional to hydraulic conductivity (K) and total hydraulic gradient (ht):
where ht usually is assumed to be the sum total of matric potential (hm, head due to pore water pressure), the head due to gravity (hg), the head due to osmosis (ho), and the head due to kinetic energy (hυ, negligible). Darcy first realized that the total energy of the system stored in the form of the total water potential is always lost during flow through the porous material and the amount of energy loss depends on the material. This is a manifestation of the second law of thermodynamics. In general, energy stored or released by an element of soil due to change in the fluid mass can be related to the hydraulic head (Freeze and Cherry 1979, p. 51):
where is the compressibility of bulk soil and is compressibility of water voids. The term is also the mean compressibility of wet soil (1 Pa−1). Whereas the compressibility of water under earth surface environment is close to a constant value of about 4.4 × 10−10 Pa−1, the compressibility of clay soil ranges from 10−6 to 10−8, of sandy soil from 10−7 to 10−9, and of pebbles–gravels from 10−8 to 10−10. Assuming constant water density in Eq. (B1) and substituting in Eqs. (B2) and (B3) gives
Equation (B4) is applicable to saturated and unsaturated flow and is the most general form of the Richards equation. For unsaturated situations, the compressibility of the bulk soil usually is smaller than the compressibility of the fluid voids, so the void compressibility dominates the average compressibility. By definition
In the absence of osmotic head, total head composes only matric suction head and geopotential, gives the most common form of the Richards equation (3D hydraulic conductivity in the component form):
For 3D landslide models, the land surface model subcomponent must solve the Richards equation in a full 3D setting [Eq. (B7)]. For example, SEGMENT-Landslide includes the Simulator for Hydrology and Energy Exchange at Land Surface (SHEELS; Ren et al. 2004) as a subscheme to simulate the soil moisture content. For application in 2D (x–z plane) slope models (e.g., TRIGRS; Baum et al. 2008), the Richards equation in hill slope settings is proposed based on Philip’s (1991) assumption, although several recent theoretical and experimental studies question the validity of this assumption (e.g., Jackson 1992; Philip 1993; Sinai and Dirksen 2006). Following are two analytical solutions to the Richards equation because of their value in dealing with flux boundary conditions. The Green and Ampt (1911) infiltration model is an analytical solution for flow into an initially dry, uniform column of soil. Although not a direct analytical solution of the Richards equation, it has the essence by two assumptions arising from observed wetting front during infiltration: 1) the soil suction beyond (ahead of) the wetting front is constant in space and time, and 2) the water content and the corresponding hydraulic conductivity of the soil behind the wetting front (the wet portion of soil column) also are constant in space and time. For some efficient land surface schemes, this relation is utilized for unsteady rainfall flux upper boundary conditions (e.g., Chu 1978). The infiltration scheme inside TRIGRS is also a unique analytical solution of Eq. (B7). It is based on a Gardner (1958) linearization of a 1D Richards equation applied to a water table. An analytical solution for transient infiltration above the groundwater table is thus obtained (Srivastava and Yeh 1991; Baum et al. 2008). Suction head and soil moisture profiles, as a function of vertical elevation above water table and of time, are standardized output of this type of model.