The Gravity Environment of Zhouqu Debris Flow of August 2010 and Its Implication for Future Recurrence

This study investigates the geological background of the August 7-8, 2010 Zhouqu debris flows in the northwestern Chinese province of Gansu, and possible future occurrence of such hazards in the peri-Tibetan Plateau (TP) regions. Debris flows are a more predictable type of landslide because of its strong correlation with extreme precipitation. However, two factors affecting the frequency and magnitude of debris flows: very fine scale precipitation and degree of fracture of bedrock, both defy direct observations. Annual mean Net Primary production (NPP) is used as a surrogate for regional precipitation with patchiness filtered out, and gravity satellite measured regional mass changes as an indication of bedrock cracking, through the groundwater as the nexus. The GRACE measurements indicate a region (to the north east of TP) of persistent mass gain (started well before the 2008 Wenchuan earthquake), likely due to increased groundwater percolation. While in the neighboring agricultural region further to the north east, there are signal of decreased fossil water reservoir. The imposed stress fields by large scale increase/decrease groundwater may contribute to future geological instability of this region. Zhouqu locates right on the saddle of the gravity field anomaly. The region surrounding the Bay of Bangle (to the southeast Corresponding authors.


Introduction
On August 8, 2010 in the northwestern Chinese province of Gansu, more than 1000 people died when a debris flow devastated the small county of Zhouqu. The future occurrence frequency of this type of landslides is an immediate practical concern for mountainous regions. Landslides are partly an upscale process, with localized disturbances adding instability, by adding fluids or reducing root reinforcement [1], to pre-existing weathered regolith and granular soil particles spreading on slopes in a region predisposed to geo-hazards. In this study, the background of the Zhouqu landslides (Figure 1) is addressed. First, the geological and hydrological circumstances of this event are examined. Remote sensing observations over the past decade provide an opportunity to examine the bedrock cracking conditions over the peripheral Himalayan region, which includes both Zhouqu County and the 2008 Wenchuan earthquake affected areas.
Debris flow landslides [2] are more predictable because of their strong correlation with extreme precipitation. The 2010 extreme precipitation was caused by a mid-latitude frontal system, manifested as a shear zone on the 700 hPa flow fields, due to the high average elevation. The question "Why was the August 2010 Zhouqu landslide so powerful?" is addressed by using the high resolution precipitation output from Weather Research and Forecast (WRF, an advanced numerical weather prediction model [3]) to drive an advanced landslide modeling system (SEGMENT-Landslide [4]- [6]). The model simulates the entire life cycle of individual landslide and allows the spatial distribution of sliding material to be estimated.
Future extreme precipitation events can be examined using WRF and a recently developed technique [7] for perturbing the thermal background, to mimic a warming climate environment. The "pseudo climate change" approach is adopted to estimate possible changes in precipitation in a 2060-80 time-frame. The characteristics of future "Zhouqu" types of debris flows then can be determined from the landslide model.

Model and Methods
How much sliding material can be mobilized on an unstable slope is co-determined by slope hydrology and precipitation morphological characteristics (i.e., total amount as well as intensity distribution). The distribution of infiltrated/drained into the ground is a key parameter for estimating the landslide magnitude. When precipitation rate exceeds infiltration rate, there will be surface runoff. The surface runoff, when meets crevasses, percolates into deeper depth, creating a shear zone much deeper than uniform wetting can reach. The landslide magnitude thus is much larger for the intense storm case (cf. Figure 2(b) and Figure 2(c)). Also, it is possible that the root reinforcement effects do not work simply because the shear zone becomes far deeper than the root zone and there is a decoupling of the two. Figure 2 illustrates how rainfall distribution within the slope affects the landslide magnitude. The weight of rainfall is negligible in aiding the driving stress. However, if it can be effectively channeled into the interface between non-fracture bedrock and the overlain granular lump, the effect on reducing the maximum resistive stress is significant. In addition, the concentrated water column inside the crevasses will exert a hydrostatic horizontal direction differential pore pressure, aiding the driving stress in overcoming the resistive stress. Once the driving stress and the differential pore pressure together overtake the reduced (by drainage lubrication) maximum available resistive stress, the entire layer is set into motion; even through most soil inside this layer still is dry. In contrast, if the precipitation is retained entirely in the upper shallow layer and did not discharge into deeper layer, only a very shallow layer is mobilized and this layer, by entraining drier material, quickly loses its fluidity and may stop mid-slope. Thus, any future changes in precipitation morphology and   structure characteristics may have direct effects on storm-triggered landslides.

Method for Bedrock Crack Estimate
At present, there are no direct measurements of bedrock crevasses. However, the bedrock crevasses are "container" for ground water. Groundwater variations are measurable by gravity satellite (e.g., Paper [8]) such as the Gravity Recovery and Climate Experiment (GRACE [9]). In the global water cycle, precipitation over land generally exceeds evapotranspiration, the excess water usually as river runoff getting back to oceans to close the cycle.
As bedrock cracks open wider, more runoff will be diverted into drainage percolation and saved in groundwater reservoir (Equation (1) in Paper [10]). This is accompanied by local mass variations detectable by GRACE. In addition to net precipitation and drainage, other factors also may cause perturbations on regional gravity, such as irrigation projects and large-scale population and animal migrations. The region of interest is a remote area in China that has not been affected by these factors. Precipitation also is well known for its patchiness, especially over complex terrains [11]. In situ measurements typically are unable to "catch" the "hot spots". To close this dominant term in regional mass balance, the annual mean Net Primary Production (NPP), because of its high sensitivity to annual precipitation, is used as an indicator for precipitation trends [12] in our analysis of the geological background of landslides.

Method for Future Extreme Precipitation Estimate
Debris flows are more predictable than other form of landslides because of their strong correlation with extreme precipitation [13]. To resolve the mesoscale evolution of the August 8, 2010 event, the WRF uses a larger outer domain with 20 km grid spacing, along with a 4 km nested domain, and a 1 km inner-most domain (Figure 1(a)), within which SEGMENT-Landslide is actually run. Initial and lateral boundary conditions were based on Global Forecast System (GFS) final analyses on a 1 degree latitude-longitude grid; simulations were initialized at 00 UTC 04 August 2010 run for 6 days. For the nested domains, lateral boundary conditions are provided by the parent domain and one way nesting is used. Detailed model setting can be found from Paper [14]. Nested run of the weather model WRF, for 6 continuous days, starting from 00 UTC 04 August 2010, with the innermost domain set at 1 km resolution clearly indicates that the strong precipitation is caused by a frontal system that manifest as a 700 hPa wind shear zone of significant cyclonic vorticity, starting at about 06 UTC August 7. There is a series of terrain forced meso-γ scale cyclones developing around Zhouqu. These features are readily identifiable from the MODIS cloud top imageries at around 0635 UTC August 7 aboard the NASA EOS Aqua satellite (figures not shown).
In a future warmer climate and an enhanced hydrological cycle [15]- [19], more extreme weather conditions are expected. Lackmann [7] presented a method for estimating specific effects from climate warming on a certain pattern of extreme precipitation events. For simulating future (2060-2080) landslide events, the WRF is used in the same pseudo-global warming approach as proposed by Lackmann [7]. There is no simple way to anticipate how changes in future synoptic patterns might influence the frequency and severity of extreme precipitations that might have landslide consequences. The experimental design for simulating future extreme precipitation events addresses this issue by replicating the synoptic pattern that caused the August 7-8 precipitation, but with likely future large-scale thermodynamic changes imposed. While it is acknowledged that an identical synoptic pattern would not occur in the future, it is reasonable to assume that a similar pattern could occur. Because the model simulation is allowed to evolve dynamically for the duration of the synoptic event, the resulting changes, brought about by imposing thermodynamic changes, can include dynamical changes as well. Quantification of projected thermodynamic change due to increased anthropogenic greenhouse gases is accomplished using an ensemble of Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4) GCM simulations of monthly atmospheric parameters, under the A1B scenario. The August monthly fields were averaged among a subset of five GCMs (same as those used in [14]). A two-decade average temperature change was then computed from the August monthly averages from the 5-member GCM ensemble for the periods  and (2060-80). All available variables on isobaric levels were averaged. These spatially varying average fields were interpolated to the same 1˚ GFS grid as used for the initial and lateral boundary conditions in the control simulation. Differences between the 2060s and 2000s decades were computed for each grid cell, and these changes were subsequently added to the original (2010) GFS analyses. Given that the landslide event in question took place in August 2010 and that a change is added that was calculated from a 60-yr difference, the future simulation can be interpolated as approximating the thermodynamic conditions around the year 2060 that would arise from anthropogenic greenhouse gas forcing alone, under the A1B scenario. The same synoptic weather pattern from August 2010 also characterizes the future simulation, but this fictitious date does not correspond to a specific date in future.

Results
Nested WRF run with reanalysis lateral boundary conditions can successfully reproduce the precipitation distribution and intensity. SEGMENT-Landslide driven by the WRF provided meteorological parameters successfully reproduced the entire life cycle of the August 7-8, 2010 mudslide [13]. Here, we address a wider issue: the geological background of this event and future recurrence. Slopes are modified by slow and steady diffusive (e.g., natural weathering processes that generate sliding material) and fast advective (e.g., debris flows that transport sliding material) processes. The drastically different time scales determined that at any moment, landslides can occur only on a small fraction of the natural slopes (limited by sliding material availability and triggering mechanisms). In other words, most slopes are accustomed to environmental conditions (e.g., vegetation, chuting system and precipitation climatology) and in a tricky balance state. Drastic changes to any link may result in sliding. For regional area, changing of underground chutting system set the background and further changes in precipitation realizes sliding.
The GRACE measurements (2003-present) over the peri-Tibetan Plateau (TP) show systematic and (on ~100 km spatial scale) consistent signals of mass loss (blue color in Figure 3) and mass gain (red color shades in Figure 3). First, the significant mass loss from the Himalayan region is from the increased ice mass loss in recent decades. This melt water flows primarily into the rivers and discharged into the Indian and Pacific Oceans, with a turnover time generally of less than a week. For soil mantled slopes, loss of soil moisture is not a concern Figure 3. GRACE measured mass changes over the 2003-2013 period (in water depth equivalent, cm). For convenience, it is subdivided into four nested domains (d01 -04 or, as defined in the text, Regions 1 -4). Regions 2, 3 and 4 all show significant mass gain signal over the past decade, for different physical reasons (discussed in the text). Zhouqu (red star) located at the saddle of the gravity anomaly field (composed by increased ground water storage at Regions 2 and 3 and decrease in water storage over Himalayans and North China Plateau over the entire GRACE measurements period). Tengchong (Yunnan province, China, "Tc" in Region 4) has a very similar situation as Zhouqu. In a larger context, the entire region surrounding the Bay of Bengal is under horizontal shear placed by groundwater changes. The black dots are observed scarps. Blue crosses are SEGMENT estimated potential locals of storm-triggered landslides in the upcoming 30 years (stable at present but will slide in the future). Bold arrows schematically indicate the superimposed stress through groundwater temporal variations. for sliding. Only for those regions above 5000 m elevation (un-inhibited regions), removal of the inter-boulder ice may cause unstable configuration (granular theory can be found in [13] [20]). The three surrounding domains all have significant mass gains, much larger than explained by geological plate movements, which global wise only of 2 GT magnitude (~1 km 3 in volume), which is three orders of magnitudes smaller. There also have been no known large-scale population and animal migrations in the past decade. The only possible cause is in the water/fluids mass changes. Region 3, which is the primary agricultural region of China, with crops comprising the major vegetation over the region, may exhibit strong seasonal fluctuations in its biomass but annual variations still are a good representation of the precipitation amount. Region 2 is an arid region, which is clearly shown in the area total biomass (NPP in Figure 4). Region 2 has a similar total area as Region 3 but the total NPP is two orders of magnitude smaller. Although the total mass in Region 2 significantly increased during the 2003-2012 period, the NPP, as a good indicator of precipitation, did not. New data from MODIS indicates a small decreasing trend. In Region 4, the precipitation showed a significant decrease during the past decade (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), in strong discord with the mass increase trend. The precipitation did not increase for Regions 2 and 4 (and it reduced in Region 4), the increased mass can only be explained by an increased proportion into groundwater reservation rather than into runoff. Region 2 is less populated and is not a major agricultural region in China. Irrigation would only cause mass loss, which is the opposite of the observed mass increase and any attempts to recharging aquifers are not planned. This trend possibly can be a result of increased ground disintegration (i.e., cracks in bedrocks). The reason is not clear but may also have a natural component after the 2008 Wenchuan earthquake, considering that the signal of a mass increase starts at least from 2003, five years earlier. Whatever the causal factors for the disintegration of the bedrock over the two domains, there is reason to remain on alert for possible subsequent natural hazards that is impending. Landslides partially are an up-scaling process, with local sloping and geological features enhance the occurrence in a larger area pre-disposed for sliding. The growing bedrock crevasses and the accumulated fluids as groundwater work together making future landslides more severe and more frequent. Looking further ahead, the pattern of mass changes creates a saddle structure at the vicinity of Regions 2 and 3 and Zhouqu sits right on the saddle (Figure 3). The exerted stress may cause strain build-up and cause further earthquakes (by reducing the natural occurrence period). This composes a positive feedback to secondary geo-hazards such as flash floods and debris flows. Note that the mass gain trend for Region 2 and 4 began well before the 2008 Wenchuan earthquake, even before 2003, the starting time of the GRACE measurements. This points to possible geological causes, because there are no large hydraulic projects that can affect persistently the region of interest. To prove this hypothesis, WRF was run over entire Domain 1 at 5 km resolution continuously from 2000 until 2012. SEGMENT-Landslide, driven by WRF provided atmospheric parameters, simulated the unstable locations (red, cross markers in Figure 3; for clarity, only those involving sliding material larger than 10 4 m 3 are shown). Dots are observed landslide events. The ones in mass loss region over Himalayan are due to melt of permanent frozen soils. In SEGMENT-Landslide, the bedrock's degree of fracture is parameterized as a function of mass increase rate. The rainfall-triggered ones all lie within regions of apparent mass increase. It also is interesting to notice that the Zhouqu mudslide is only one among many cases occurred within the same saddle region. We notice that the region over Hengduan Mountains (Region 4) also is a hot spot for storm-triggered landslides (red crosses are simulated cases in the past decade). Due to lack of geological conditions, SEGMENT-Landslide skipped the regions within India and Burma. From the pattern of the gravity fields, we believe the entire region surround the Bay of Bangle is at present active in storm-triggered landslides.
The same precipitation total may mobilize a very different amount of sliding material, as indicated in Figure  2. The morphological characteristics of precipitation are critical for future landslide magnitude. To estimate possible changes in extreme precipitation events, WRF was run with thermodynamic changes computed from a 5-member IPCC GCM ensemble output, under the A1B emission scenario applied to the initial and lateral boundary conditions for this event. The resulting simulation is a replication of a highly similar synoptic pattern to that of the control simulation, but within a projected future thermodynamic environment. Under future warmer climate, the total precipitation over the region of interest (Domain 3 in Figure 1(a)) slightly increase (~1.6%). More importantly, the spatial and temporal distributions become even more uneven (the maximum rainfall intensity, not collocated with control run, can be 30% in difference; the maximum total precipitation, again occurs at different model grids, can differ as much as 12%). This favors more severe storm-triggered landslides.
Details of the projected future August 8, 2010 like storms over Zhouqu region are discussed in Paper [14]. In summary, the increased atmospheric water vapor as a straightforward thermodynamic consequence of warming, or a robust vapor increase consistent with the Clausius-Clapeyron relationship, holds [17] [21]. In addition, precipitation in Zhouqu area also is strongly affected by orographic effect. For example, the northern branch of the circumventing Qinghai-Xizhang Plateau airstream adds cyclonic vorticity to the frontal system. Due to the complicated topography and >3000 m surface elevation, there are no closed surface lows and there is only a shear line in the flow fields. The region in question is located within a frequent cyclogenesis region to the east of the TP; there are two such hot spots, and Zhouqu is within the northern one. In this sense, extreme precipitation at Zhouqu is produced by locally enhanced convective storms. Using WRF down-scaled ensemble climate model provided atmospheric conditions for the upcoming 50 years, SEGMENT-Landslide indicates that Region 2 and 3, which is to the north and northeast of TP and in the vicinity of Zhouqu is still under threat of future landslides. Region 4, primarily the Hengduan Mountain Ranges, especially the South-North oriented valleys of Nu Jiang, Jinsha, Lanchang and Yalong rivers (Mainly in Yunnan, Guizhou and Sichuan Provinces of China), also is under potential threat. In a larger context, the Bay of Bengal rim is under horizontal shear placed by groundwater changes and is prone to storm-triggered landslides in the upcoming years. The green markers are locations that are stable at present but will experience storm-triggered landslides before 2060.

Summary
Surrounding the TP, there are a train of regions suffering steady mass gain in the past decade. It is deduced that the mass increase in at least two regional areas is likely due to increased percolation of surface runoff, or re-charging of the fossil water reservoir, rather than an increase in precipitation amount. The reason for increased partitioning of precipitation into groundwater is not entirely clear but an increase in bedrock crevasses is one viable explanation. These regions of mass gain contrast with the mass loss over TP and North China Plateau and superimpose a horizontal stress field that favors geological hazards. Zhouqu and Wenchuan both sit on the saddle point of this pattern, likely not entirely by coincidence. As the gravitational signal starts well before the 2008 Earthquake, we do not think that the earthquake causes the crevasses enlargement, rather the changes in gravitational fields may have set the stage for both hazards. To make the situation worse, these two regions of ever-increasing fracturing are co-located with the atmospheric batroclinic cyclogenesis basins corresponding to the circumventing TP airstreams. As the climate warms, analyses indicate that extreme precipitation events are likely to become more common and future recurrence frequencies of such events therefore likely increase.