The recent increase in reports of mass removal processes (MRPs) worldwide has been attributed to causes such as global warming and climate variability, deforestation, the increase of population and the extension of human settlements into risky areas. However, it is not clear whether the increase in reports reflects an increase in episodes or settlement of population in unstable areas with a risk of disaster.
MRPs are triggered by a variety of both natural and anthropogenic occurrences, such as earthquakes, heavy rains, volcanic eruptions, exploding mines, collapse phenomena and others. Precipitation is usually the most common cause, because water infiltrates into the pores of the soil and leads to the instability of a slope  . Consequences can be devastating in areas with extreme weather events such as tropical cyclones or severe convective systems.
In Mexico, owing to its geographical location, complex topography and the high levels of social vulnerability, there are areas with a substantial record of disasters caused by MRPs, usually triggered by heavy rains  . The Sierra Norte de Puebla has characteristics that enhance the risk of disaster: a high poverty rate; precipitation ranging from 600 mm to >4000 mm per year in the highlands  and the rugged and mountainous terrain reflecting its location among three physiographic provinces, the Sierra Madre Oriental, the Plain Gulf Coast and Mexican-Trans Volcanic Belt  . Disaster prevention and adequate warning systems require an understanding of landslide mechanisms. The relationship between landslides and precipitation has been studied in order to obtain a threshold value  . For Gill and Malamud  , a threshold is the minimum amount of energy for a major hazard to trigger a secondary risk, and it can be defined in terms of factors such as soil, precipitation and hydrological conditions. The thresholds most used are those calculated probabilistically. Many countries have thresholds mapped for each region considering variables such as soil type, climate and geological characteristics  ; in Mexico, however, this information is scarce.
The aim of this paper is to assess a probabilistic precipitation threshold for landslide processes in the Sierra Norte de Puebla, Mexico. We hypothesize that the occurrence of landslides in the study area is mainly related to the characteristics of rainfall events (duration and intensity), while antecedent rainfall seems less influential. We used a Bayesian approach  . An important result was a regional amendment to the warning system proposed for Italian conditions by Brunetti  ; this should assist in warning local populations in the Sierra Norte de Puebla to evacuate the area and avoid disaster under certain storm characteristics.
2. Methods and Data
2.1. The Bayesian Method
Of the various statistical methods to obtain thresholds, a probabilistic model is more reliable than a deterministic one, since failure of slope stability is determined not only by rainfall but by a combination of numerous other physical parameters (e.g., soil thickness, groundwater conditions, shear resistance). The deterministic approach implies that there is no randomness in the development of future states of the system, and that there is one value (above or below the threshold) output for any given value input. These models require calibration over a well-specified type of failure, and detailed knowledge of all input parameters is generally difficult to obtain  . An approach more suited to the physical conditions of the Sierra Norte de Puebla is a regional threshold, since it is a direct application of conditional probabilities and allows stake holders to make decisions regarding risk  .
2.2. Sierra Norte De Puebla Landslides
The Sierra Norte de Puebla (19.47˚N - 20.83˚N, 97.12˚W - 98.32˚W) is a region formed by 63 municipalities in the north of the state of Puebla, bordering the states of Veracruz, Hidalgo and Tlaxcala (Figure 1). For planning purposes it is divided by the National Institute for Federalism and Municipal Development into two socioeconomic regions: the Sierra Norte and the Sierra Nororiental.
The present study reviewed landslides during 1989-2012 related to precipitation; sources consulted were the Disaster Effects Inventory System and also local and regional newspapers. Landslides included some of great magnitude that
Figure 1. Spatial distribution and number of mass removal processes events in the Sierra Norte de Puebla from 1989 to 2012, and of the weather stations used for the analysis. Digital elevation model 30 m and degree of slope 0˚ to 83˚.
caused considerable damage: the intense precipitation associated with tropical depression 11 of 1999 caused hundreds of MRPs and other effects in more than 27 municipalities, leaving dozens of communities cut off and more than 400 fatalities  ; heavy rains in September 2006 were followed by the collapse of Necaxatépetl hill in the municipality of Juan Galindo onto two microbuses, three trucks and one van, leaving 4 dead, more than 15 persons injured, and 2 missing; and Hurricane Dean (category 5) caused gravitational processes in 8 municipalities on 22 August 2007 when 74 municipalities of Puebla were declared a disaster area.
In total, 163 events occurred in the region during 1989 to 2012. Most MRPs occurred in areas with steeper slopes at the intersection of the Trans-Mexican Volcanic Belt and the Sierra Madre Oriental. Of the 44 municipalities with MRPs, the municipality of Huauchinango had 17 major movements triggered by rains, followed by Zacapoaxtla and Teziutlán in Sierra Nororiental, and Zacatlán in the Sierra Norte (Figure 1).
To calculate regional probability, those municipalities that had two or more movements on the same day but in different locations were recorded as having had a single movement triggered by the same event of rain. Approximate coordinates of each movement were deduced from descriptions in the local media. Details of the precipitation were obtained from weather stations of the Mexican National Weather Service and of the National Institute of Agricultural and Forestry Research; only stations with more thirty years of operation and more than 85% precipitation records were considered. Digital elevation model variables and degrees of slope were calculated for an overview of the terrain which served to map the area of influence of each weather station and because slopes more the 28˚ are vulnerable to abrupt and numerous landslides associated with heavy rains  . A maximum distance of 10 km from a MRP to the nearest weather station was defined in order to report rainy conditions that triggered the movement.
In total, 15 stations were selected on the basis of data quality, proximity to the MRP, relief and homogeneous distribution in the region (Figure 1). Hourly rainfall was analyzed for each of the 163 events, to determine the duration, and the cumulative rainfall in the days preceding the event. Such accumulation was defined visually in a graphic starting from three days of <0.5 mm of rain. The rain thresholds were divided into three types following Berti  : 1) well defined (Figure 2); 2) uncertain (Figure 3); and 3) indefinite (Figure 4). From the intensity and duration of rain the probability for each conditional scenario was calculated for each of 15 selected stations, together with a total conditional scenario considering all landslides that occurred in the study region and taking the highest precipitation value reported by station. The procedure was applied for one and two dimensions. Finally, with thresholds for both dimensions expressed in millimeters of rain and duration in days, an alert system was proposed as in Brunetti  .
Figure 2. Example of a well-defined trigger by rain. The red arrow indicates the day of the MRP occurrence. The threshold is appreciably more than 160 mm of accumulated rainfall.
Figure 3. Example of an uncertain triggering of rain. The red arrows indicate on 9 and 10 July 2008 MRP occurrences. The threshold is observable with 120 mm of accumulated rainfall.
3.1. Characteristics of Precipitation
Occurrence of MRPs was highest is July, followed by October, with minor incidents in May and November, the beginning and end of the rainy season. The rainiest months are July and September. The year 2008 was selected because it presented various hydro-meteorological phenomena associated with intense rains such as 8 tropical cyclones, 52 cold fronts, 37 tropical waves and five winter storms. Two representative surface stations were used to obtain the monthly accumulated precipitation. Although there is monthly variability, precipitation
Figure 4. Example of an indefinite triggering by rain. The red arrow indicates the MRP occurrence. The dotted line indicates the sharp increase in rain without movement.
usually occurs between 15:00 and 22:00 hours local time; whereas in the Sierra Nororiental region (Teziutlán station) between 4 and 5 pm and in the Sierra Norte (Huauchinango station) between 7 and 8 pm.
3.1.1. Tropical Cyclones
Summer rains are associated with the Inter Tropical Convergence Zone (ITCZ), the Mexican monsoon, the easterly waves, and hurricanes from the Pacific and Atlantic oceans  . Hurricanes are the strongest influence on the abrupt and recurrent increase in the intensity of precipitation. The MRP events triggered by a tropical cyclone according to historical information provided by the Mexican Meteorological Service were classified in terms of the El Niño/Southern Oscillation (ENSO) years (Figure 5). Of the 164 landslides, 107 (65.24%) were directly related to tropical cyclones, 41 to the rains of tropical depressions, 26 to hurricanes, 22 to tropical storms, 10 to tropical waves and 8 to cold fronts. This indicates the need for preventive measures when these phenomena approach the Mexican coast.
The ENSO is the most important rainfall modulator in the south-central region of Mexico and is linked to the MRP events to some extent. Accumulated rainfall in La Niña years is significantly higher than normal, and in El Niño years it is lower; both are positively correlated with MRP occurrence; of the 163 MRP events, 113 occurred during La Niña, 34 in normal years and 16 in El Niño years (Figure 6).
3.2. Probability Calculation
3.2.1. One Dimension
With precipitation data, the approximate threshold that triggered slope movement was determined; 67 MRP events were classified as well defined, 46 as
Figure 5. Frequency of mass removal process (MRP) events in the Sierra Norte de Puebla during ENSO and neutral years. (TD, Tropical depression; H, hurricane; TS, tropical storm; TW, tropical wave; CF, cold front).
Figure 6. Tri-monthly accumulated rainfall. MJJ (May-June-July), JJA (June-July-Au- gust), JAS (July-August-September), ASO (August-September-October).
uncertain and 35 as indefinite.
Following Berti  the one-dimension Bayesian analysis (Equation (1)) used three classes of precipitation: 0 to 30 mm, 30 to 60 mm and >60 mm.
P(A) = NA/NR is the prior probability of a MRP, represented by the relationship of the number of movements in the region to the total number of rain events during 1989-2012.
P(B) = NB/NR indicates the marginal probability represented by the number of rain events registered in a time period in each of the three ranges defined among total rainfall events.
P(B|A) = N(B/A)/NA is the conditional probability of B given A, represented by the MRPs during each of the three ranges defined by the total number of MRPs in the region during the period.
P(A|B) = N(B/A)/NB is the conditional or a posteriori probability, i.e. the probability of A given B, and is represented by the number of MRPs during each of the three ranges defined by the number of rain events in each of the ranges.
The total conditional probability (Figure 7) for rain rates of 0 to 30 mm, i.e. the probability that a movement will occur when it rains with this intensity, is
Figure 7. The conditional probability at three rainfall intensity ranges for each station and total (considering all landslide events and the higher precipitation value reported in the study region).
<0.01. The number is low, because most rain events in the region are <30 mm and very little movement arose during those days. For daily intensities of >30 to 60 mm, the probability is higher but still <0.1. In contrast, when the intensity is >60 mm, the probability is almost 0.5.
However, this one-dimension calculation omits the duration of the precipitation, a limitation overcome by a two-dimensional calculation.
3.2.2. Two Dimension
Bayes’ equation allows incorporation of the two variables of interest within the precipitation for the joint probability of triggering a MRP  . Other variables can be used in place of as previous rain; however, this will depend on prior interpretation to show a better explanatory power. Previous visual analysis of the data concerned duration and intensity of rain (Figure 8). Conventional thresholds are generally obtained in this way, establishing a lower limit of precipitation events that trigger movement, or by statistical formulas that result in a linear trend   . The threshold lower envelope proposed by Jibson as 
where I is the intensity (mm/day) and D is the duration (days) defines a line below which the values do not trigger any MRP (Figure 9). However, a global threshold may have more uncertainty. A global threshold proposed by Caine as 
Figure 8. Analysis of two-dimensional Log-I, Log-D. Puebla, proposed for the Sierra Norte de Puebla visually designed threshold, lower envelope (Equation (2)) and world global threshold (Equation (3)).
Figure 9. Probability calculation for intensity and duration of rain. The vertical axis indicates the probability calculated for each of the scenarios. The dotted line indicates the threshold that increases the likelihood of triggering movement in the Sierra Norte de Puebla. The white boxes indicate the lack of data for calculating probability and the red boxes the null probability of landslide occurrence.
clearly defines a zero probability of movement below the threshold but has the limitation of generating very general intensity values corresponding to <100 mm for one day. Therefore, it is proposed to raise the threshold for the Sierra Norte de Puebla, in order to approach the likelihood of triggering.
The equation (Equation (4)) was calculated for each of the stations with reference back to the joint base. There were nine ranks, 3 × 3: durations of 0 - 10 days, >10 to 20 days and >20 days; and intensity of 0 - 30 mm, >30 to 60 mm and >60 mm, resulting in nine combinations for conditional calculation (Table 1).
P(A) is the a priori probability; the number of MRPs for total rainfall events.
P(B,C) is the marginal probability; the number of rain events in the first intensity range (0 to 30 mm/day, corresponding to the first quadrant) and duration of 1 to 10 days by total rainfall events during that period.
P(B,C|A) is the conditional probability in 1 to 10 days by the total number of movements recorded in that period.
P(A|B,C) is the posterior or conditional probability; the number of registered movements when a rain event had an intensity of 0 - 30 mm and duration of 1 to 10 days, the total rainfall events of such intensity and duration that occurred in the period (Table 1).
The Bayes probability is calculated for each conditional scenario (Figure 9),
Table 1. Two dimensional probability matrixes for the Sierra Norte de Puebla. B, C indicate the joint probability of having a certain value (or range of values) of the two variables; the rainfall magnitude and the rainfall duration. The P(B,C|A) is the conditional probability of B given A and C, in this case the probability of having been a rainfall event of magnitude B and duration C when a landslide occurs; P(A) is the prior probability of A (or simply prior), in this case the probability that a landslide occurs regardless of whether a rainfall event of magnitude B and duration C occurs; P(B,C) is the marginal probability of B, in this case the probability of having been a rainfall event of magnitude B regardless of whether a landslide occurs; and P(A|B,C) is the conditional probability of A given B and C, in this case the probability of a landslide when a rainfall event of magnitude B and duration C occurs.
with higher probability of MRP occurrence when the rainfall exceeds 60 mm and from 1 to 10 days; in turn, the probability decreases with precipitation of 0 to 30 mm and from 1 to 10 days. Also, probabilities are nearly 0 for rainfall >30 mm at 20 to 30 days. This result could help decision makers when there is a high probability of precipitation of over >60 mm for >1 day. Also, a proposed new threshold is determined for the Sierra Norte de Puebla that best defines the probability of landslide according to the intensity and duration of rainfall.
The relative ease of the method allows incorporation of variables other than precipitation, such as soil type.
4. Warning System
The most important step in disaster risk management is prevention, so the present results should be incorporated in a warning scheme based on that proposed by Brunetti  . Five probability ranges are defined as guidance to the decision makers in the region (Figure 10).
First, in those regions of the sierra that increase in risk after a warning of a tropical cyclone, or in the rainiest months (particularly in La Niña years), or just with a forecast (numerical weather prediction, satellite images) of steady rain, the rainfall for the next 24, 48, 72 and 96 hours is calculated for each of the most significant 15 stations selected for their homogeneous distribution and reliability. According to the intensity thereby deduced for each time, the conditional probability is determined and classified into one of 5 categories: 1) well below the threshold, 0 < 0.029; 2) under the threshold, 0.0291 < 0.25; 3) on the threshold, 0.251 < 0.39; 4) above the threshold, 0.391 < 0.5; and 5) well above the threshold, 0.51 < 1.0. Each site is assigned a class number corresponding to a greater or lesser MRP risk. In those sites with values above or well above the
Figure 10. Example of landslide alert bulletin for the Sierra Norte de Puebla. Top left, location of 15 selected stations. Low left corner, classes of probability of a mass removal process. Top right, prediction of rain for the next few hours for each station selected; ↑ (increase), ↓ (decrease) and the correspondent assigned class of probability.
threshold, a higher risk of disaster is expected and an alert bulletin is issued to the appropriate authorities. This process can be performed automatically to work in real time with the information provided by the National Weather Service.
5. Discussion and Conclusions
The advantages of the Bayes method are as follows: application to almost any natural phenomenon in which a risk is involved; versatility to incorporate into a study the required variables; ability to develop appropriate thresholds at the scale required; and the relative ease of the method to integrate the uncertainty. Although the method has been applied in several countries with high reliability, in Mexico there have been no studies using this approach.
In this study, the two most important variables (intensity and duration of rain) in rainfall-related triggers of mass removal were considered in order to obtain a priori, marginal and largely conditional probabilities, since the uncertainty is crucial in the development of more suitable disaster prevention measures.
In general, the main limiting threshold factor is the obtaining of information and choice of variables. Inclusion of other relevant physical variables (slope, soil type, soil humidity, etc.) in probability calculation could enhance the study. However, the probability calculation is more complex because the additional information is often not available. Also, the scale can be seen as a limiting factor; although it would more accurate to determine a local threshold for each municipality, or regional one for the Sierra Norte and Sierra Nororiental, however there are no enough weather stations with reliable data especially in Puebla, despite their long history of disasters caused by rains. A better option may be to use precipitation data from remote sensing sources (satellite or radar).
Future studies could include more ranges of intensity and duration of precipitation and information about slopes, areas of population, vegetation, and soil mechanics inter alia, performing a multi-criteria evaluation.
The location of the Sierra Norte de Puebla among three physiographic provinces as well as its proximity to the Atlantic Ocean provided the conditions that increase the risk of MRPs and disaster for a highly vulnerable population. Extreme hydro-meteorological events during La Niña years intensify and prolong the precipitation; this is strongly correlated with MRP triggering so obtaining thresholds of these variables is of great importance. The Bayesian method allows conditional probability thresholds to predict future landslides according to the precipitation forecast, and these can be applied in warning systems in the region. The threshold classification can be modified to improve accuracy and thereby to prevent material and human losses.
This research was performed in 2012-2015 while the first author was preparing the undergraduate thesis at the Institute of Geography, UNAM, and was partially sponsored by the María Teresa Gutiérrez de MacGregor Fellowship. Dr Matteo Berti is gratefully acknowledged for making software available for checking our results and plotting Figure 9. Special thanks are owed to M. Sc. Luís Galvan for his help with GIS in producing some plots and slope calculation and M. Sc. Armenia Franco for precipitation data.