Probabilistic Risk Assessment of Livestock Snow Disasters in the Qinghai-Tibetan Plateau

Understanding risk using a quantitative risk assessment offers critical information for risk-informed reduction actions, investing in building resilience, and planning for adaptation. This study develops an event-based probabilistic risk assessment model for livestock snow disasters in the Qinghai-Tibetan Plateau (QTP) region and derives risk assessment results based on historical climate conditions (1980–2015) and present-day prevention capacity. In the model a hazard module was developed to identify/simulate individual snow disaster events based on boosted regression trees. Together with 15 a fitted quantitative vulnerability function, and exposure derived from vegetation type and grassland carrying capacity, risk metrics based on livestock mortality and mortality rate were estimated. In our results, high risk regions include the Nyainqêntanglha Range, Tanggula Range, Bayankhar Mountains and the region between the Kailas Range and neighbouring Himalayas. In these regions, annual livestock mortality rates were estimated as > 2% and mortality was estimated as >2 sheep unit/km at a return period of 1/20 a. Prefectures identified with extremely high risk included Yushu in Qinghai 20 Province and Naqu, Shigatse, Linzhi, and Nagri in the Tibet Autonomous Region. In these prefectures, a snow disaster event with return period of 1/20 a or higher can easily claim a total loss of more than 200,000 sheep units. Our event-based PRA results provide a quantitative reference for preparedness and insurance solutions in reducing mortality risk. The methodology developed here can be further adapted to future climate change risk analyses and provide important information for planning climate change adaption in the QTP region. 25


Introduction
Livestock snow disasters are serious winter extreme weather events that widely occur in central-to-east Asian temperate steppe and alpine steppes (Li et al., 2018;Tachiiri et al., 2008).In the pastoral areas of these regions, heavy snow fall leads to thick and long-lasting snow cover, making forage unavailable or inaccessible (Fernández-Giménez et al., 2015).Together with extremely low temperature and strong wind, it severely inhibits natural grazing, claims considerable livestock mortality, and brings devastating impacts to the livelihoods of local herders, even threatening their survival (Wang et al., 2013a).In not report many details describing the final risk metrics.Ye et al. (2017) further extended the PRA framework to support insurance design and pricing.They focused on the risk of economic stress to local herders due to increased feeding expenditures induced by long-lasting thick snow cover in eastern Inner Mongolia but did not analyze livestock mortality issues.
In this study, we developed an event-based PRA method for present and future livestock snow disaster risk assessments for the QTP region.There are three major aims of this study: 1) Develop a hazard module that can identify/capture snow disaster event based on daily weather data.It is the basis for any event-based modelling attempts, and is particularly important for regions where historical records are absent and for future risk assessment where observations and records are not yet available and variabilities from future climate change will exist.2) Set up an event-based PRA framework for livestock snow disaster risk assessment by integrating snow disaster event (hazards), livestock vulnerability, and exposure together to derive a probabilistic quantification of risk.3) Derive the risk metrics for livestock mortality risk in the QTP and offer riskinformed reduction implications.
Worldwide, the QTP suffers from some of the highest livestock snow disasters due to its large area of snow cover area, longlasting snow cover days, and nomadic grazing.This region is also a hot spot in climate change.Quantitative risk assessments for the present day will likely be a significant source of information for disaster risk reduction.In addition, the framework can be adapted for livestock mortality in snow disasters in the context of future climate change analysis, and therefore support climate adaptation planning for local government and herding communities.

Materials and methods
We followed the PRA approach proposed by Carleton and Hsiang (2016), and applied the concept of event-based modelling in catastrophic risk models (Michel-Kerjan et al., 2013).The event-based modelling approach was framed using state-of-theart three-element risk modelling, hazard, exposure, and vulnerability (Kinoshita et al., 2018;Muis et al., 2015) to model losses claimed by individual events.Then PRA was achieved through repetition of individual event modelling, in which a large number of events were drawn from the full distribution of hazards, given the predicted losses/consequences from individual events, from which a full distribution of disaster loss can be obtained (Fig. 1).

Vulnerability function
Vulnerability is the critical function that links dose (hazard inputs) and response (loss estimates) (Carleton and Hsiang, 2016).
For livestock snow disasters, a set of vulnerability functions have been estimated linking livestock mortality (rate) to snow disaster duration, during disaster environmental stress, summer season vegetation productivity, and disaster prevention capacity (Fang et al., 2016;Wang et al., 2016).To fulfil the goal of event-based modelling, the vulnerability relationship must be built on event basis.Therefore, the results from (Li et al., 2018;Ye et al., 2018)

Hazard
Hazard module is critical in our event-based PRA method, and its needs to identify individual snow disaster events, and provide provide the exact timing (starting and ending dates) of each event, based on which Duration and Wind can be derived.Such information, nevertheless, are not so straightforward to obtain.For the historical period, there are no ready-touse snow disaster event datasets at the grid level.The number of meteorological stations capable of observing snowfall in QTP is limited and are primarily located to the eastern and southern part of the region.For future risk assessment, no projections of snow disaster events are provided in climate scenario datasets, although models have been developed to simulate daily snow depth (Yuan et al., 2016).Therefore, a snow disaster event identifier/simulator was developed here to identifying/simulating snow disasters.
A snow disaster is a weather process with snow fall, low temperature, and snow cover, with certain length of durations, according to the Chinese national standard for Snow Disaster Grades in Grazing Regions of China (GB/T20482-2017) and China Meteorological Administration (CMA) standard for Meteorological Grades of Urban Snow Hazards (QX/T 178-2013).A snow disaster event designation largely depends on the snow weather process and observer's decision (manual record).To mimic a meteorological observer's decision to designate a snow disaster event, our snow disaster event identifier/simulator has considered two major questions.First, whether a specific day would be regarded as a snow-disasterday (SDD) given weather information of the day and previous days.The key is the modelling the binary response variables (Yes/No), which can be conducted with either regression or classification methods.Second, whether two SDDs, exactly neighbouring or a couple of days away from each other, should be regarded as one snow disaster event.The key is to assemble many single SDDs into snow disaster events, which can be accomplished using smoothing and filtering.In response, three major steps were considered (Fig. 2 For this step, boosted regression tree (BRT) modelling was used to establish a multi-variate and non-linear relationship 5 between SDD and various weather information.BRT modelling methodology was chosen due to its promising power for both explanatory and predictor purpose in many ecological and environmental modelling scenarios (Elith et al., 2008;Hastie et al., 2009).Other machine learning methods such as random forest can also be used here but less likely to outperform BRT according to the literature (Oppel et al., 2012;Youssef et al., 2016).To fit a BRT model, historical snow disasters were first turned into SDD flags, if a date was included in a historical snow disaster it was flagged with "1", and "0" otherwise.10 Variables used to explain and predict days that would be considered SDD was inspired by the two standards, GB/T20482-2017 and QX/T 178-2013.We included daily snow depth (SD, cm), daily maximum (maxWind), mean (meanWind) and minimum wind (minWind) speed (m/s), daily maximum (maxT), mean (meanT), and minimum (minT) temperature (℃), daily precipitation (Pre, mm), and average daily precipitation since the last snow fall day (precipitation > 0.1 mm) (aveP,

Historical snow disaster records
Step 1: BRT modeling of Snow-disaster-day probability

Climate data (station)
Step 2: Smoothing (window size = 5~25 by 2 ) Step 3: Filtering BRT model fitting was conducted using the package dismo (Hijmans et al., 2011) in R 3.3.3.Given the type of response variable, the Bernoulli distribution family was used.To identify the best combination of model parameters, we compared the combinations of lr = (0.01, 0.005, 0.001) and tc = (1, 2, 3, 5), as recommended by Anderson et al (2016).The maximum number of trees was set to 10,000, which proved sufficient for convergence.For each combination of parameters, we applied the predictor selection process using the gbm.simplify function to best obtain a balance between prediction power and number of predictors requested.The result with the least cross-validation deviance was retained.To achieve the most promising goodness-of-fit, historical snow disaster records obtained from CMSDS (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) were used for fitting.This part of the data exactly matched the CMSDS record and most precisely followed CMA's definition of a snow disaster.
Records for 1980-2007 were used later for validation and calibration purposes.
(2) Steps 2 and 3: assembling single SDDs to events by smoothing and filtering A fitted BRT model can help predict the probability of a single day being judged as a SDD.To predict/rebuild snow disaster events, these single day probabilities must be deemed snow disaster events, an ensemble of multiple SDDs.Because the explicit output from the BRT suffered from prediction errors, simply using a threshold to turn probabilities into "0/1" values would yield a set of "busy" snow disaster events, e.g., high frequency but small duration (Fig. A1, "Step 1").Therefore, a smoothing treatment is needed to filter out isolated single SDDs and fill the small gaps between two neighbouring events.
There are two parameters essential to changing the frequency and duration of identified snow disaster events: the smoothing window size and filtering threshold.In general, using larger window size for smoothing can filter out noises and reduce the frequency of events, while using lower threshold can increase the duration of single events.In order to best match the annual occurrence and the duration of single events, the two parameters were tuned through calibration using the full dataset of historical records between 1980 to 2015.We considered moving averages with window sizes from 5 d (minimum duration of a single disaster as defined by CMA) to 31 d (one month) in steps of 2 days, in combination with thresholds of 0.10 -0.5 in steps of 0.02.The timing and duration of events derived from our model for any given pairs of window size and threshold were compared with historical records, including the frequency distribution of annual occurrence of single events, the frequency distribution of single event duration, and the timing of each single event.Through tuning, the combination of parameters that yielded the best matches were recorded.China, and the QTP has been used as a focus region for validation (Chen et al., 2011;Yang et al., 2010).The 3-h dataset was aggregated to daily for input to the BRT model to rebuild gridded snow disaster events.Based on the identified events, the variables Duration and Wind were computed as inputs to the vulnerability function.From the 35 winters' events identified, we calculated the annual frequency and mean (single) event duration of snow disasters, as well as their return period values (Fig. A3).

Exposure
Exposure measures the distribution of assets or population exposed to hazards (Kinoshita et al., 2018).In our framework, it should provide the spatial distribution of herd size, and help turn the outputs from event loss modelling and livestock mortality rate (the response variable in the modelled vulnerability function) into mortality (death toll).
A full gridded distribution map of herd size in the QTP remains unavailable.The derivation of its spatial distribution was supported by the rule-of-thumb for "forage-livestock balance," as written in the Forage-livestock Balance Management Approach.The Approach was issued by the Ministry of Agriculture of China in 2006 to mitigate severe over-grazing in the pastoral areas of China (Shang et al., 2012).Using the Approach, herd size at the county level was strictly controlled under carrying capacity computed according to a thorough investigation of local grassland resources.Therefore, a gridded carrying capacity map can be a good approximation of the actual herd size distribution.
We estimated the spatial distribution of herd size by turning grassland distribution data with the look-up  A2; Table A1).

Loss modelling
Snow disaster event losses measured with livestock mortality rate (death toll/ herd size) were modelled by taking requested inputs into the vulnerability function.Duration and Wind were outputs from the hazard module.Growing season aggregate precipitation P was computed from the climate forcing data.For GDP, we used the fixed value for 2015 obtained from the Event mortality rates were then aggregated into annual mortality rates, considering the possibility of multiple events per location annually, although unlikely.In aggregation, we assumed that the second snow disaster event can only have an impact on livestock surviving from the first event, and so on.Therefore, the annual aggregate loss rate in a given grid is  is the modelled loss rate of the ith event in a year, and N is the total number of events.The annual aggregate mortality rate can finally be turned into a death toll by multiplying exposure, the herd size in a given grid.
Event/annual mortality (death toll) can then be derived by multiplying event/annual loss rate in any given location by its herd size.For each grid, 35 annual loss records were modelled (there are 35 winters in 36 years), including both mortality and mortality rate figures.The number of event loss records differ by location, depending on the identified number of events for each grid.

Risk metrics
In the risk metrics, modelled losses of discrete event/annual losses were turned into probability distribution of losses.We followed standard risk metrics by deriving the average and return period values (Michel-Kerjan et al., 2013;Shi and Kasperson, 2015) of annual mortality rate and death toll for each grid.Due to our limited time span of repetition, return periods of 1/10 a (once in ten years, the 90 th percentile of the distribution), 1/20 a, and 1/50 a were considered while 1/100 a usually used in flood/ earthquake studies (Kinoshita et al., 2018) was not considered.The kernel density method was employed to fit non-parametric distributions to derive those return period values by grid.We used the Gaussian kernel function and its corresponding optimal window width in the fitting process according to the "rule-of-thumb" for optimality (Deng et al., 2007;Silverman, 1986).In addition, aggregate mortality rate and death tolls at municipal level were also derived using zonal statistics, so as to better validate the result with historical losses, and provide policy implications.

The trained BRT model and tuned parameters in rebuilding snow disaster events
The trained BRT model retained six variables but excluded SD, minWind, and Pre.In the final model, we used lr = 0.001 and tc = 5.It had a training data Area-Under-the-Curve (AUC) score of 0.858, and a cross-validation AUC of 0.829, indicating good predicting performance (Youssef et al., 2016).For the six variables entered in the final model (Fig. 3), maxT has the highest relative contribution (32.77%), while meanWind has the lowest (5.78%).After tuning the window size with the moving average and threshold, we found the best results with a window size of 21 and a threshold of 0.18.The derived results well captured the timing of occurrence of historical events (Fig. A1) and matched the empirical cumulative density functions (ECDF) for historical durations (Fig. 4), for both event and annual aggregate durations.In historical records, two or more events in a single year at a single location are rare.Therefore, ECDFs for historical single event duration and annual aggreation duration were quite close to each other in Fig. 3. Two-sample Kolmogorov-Smirnov tests were also conducted to verify the degree of agreement between ECDFs.For single event duration (observed vs. predicted) the test statistics was 0.138, and its corresponding p-value was 0.118.The annual aggregate duration (observed vs. predicted) test statistics was 0.131, and its corresponding p-value was 0.189.Therefore, the prediction model well captured statistical features of historical snow disaster duration and the predicted results can be used for event loss modeling.

Model-derived snow disaster events, 1980-2015
With the tuned model, the timing of snow disaster events were identified in the historical period 1980-2015.
Correspondingly, the annual occurrence frequency and duration of snow disaster events were derived (Fig. 5).In the figure, 5 non-grassland areas, including permanent snow areas, were masked using the vegetation map.Across the entire plateau, the annual average frequency was below 0.2 in most regions, i.e., on average, snow disasters occur every 5 years in these regions.
Higher frequency regions were primarily located in major mountains, including the Tanggula Range and Nyainqêntanglha Range in the central part of the plateau, and the Kailas Range and neighboring Himalayas.These regions are higher elevation and spatially close to permanent snow-covered areas.For major pastoral production regions, i.e., the Naqu prefecture in the 10 central QTP, the annual average frequency was 0.2 to 1, echoing the local proverb, "small disaster once in 3 years, and a major disaster once in 5 years" (Ye et al., 2017b).The distribution of mean event duration of snow disasters was consistent with annual frequency, indicating strong controls from elevation and topology.For most regions, mean event duration was below 7 d.For typical pastoral regions, i.e., Naqu Prefecture, a snow disaster can last for more than 12 d on average.Mean event duration can last for more than 14 d in high elevation mountainous areas, including the Himalayas in the southwest and alpine meadows to the east end of Bayankhar mountains, which is nearly 10% of the total grids with valid values.

Risk in terms of livestock mortality rate
The assessed livestock snow disaster risk measured by annual mortality rate is presented in Fig. 6.As presenting the full probability distribution of livestock mortality rate by grid is not viable, these figures include the annual average and three return-period mortality rate maps (1/10a, 1/20a, and 1/50a), upon which the non-pasture areas were masked.Spatial distributions of mortality rate at different return-periods are highly consistent (Fig. 6).The pattern is very similar to the pattern of annual aggregate snow disaster duration (Fig. 5), confirming the dominant influence of snow disaster duration.
High-mortality rate regions are primarily located in the major mountainous areas, including the Tanggula Range and Nyainqêntanglha Range in the central QTP, the Kailas Range and neighboring Himalayas in the southwest QTP, Bayankhar mountains in the east QTP, and southern part of the Kalakoram Range and the west-end of the Kunlun Mountains in the northwest corner of the QTP.Classified by administrative districts, high mortality rate regions include the Yushu and Guoluo prefectures in Qinghai Province and Naqu, south Ngari, and Linzhi Prefectures in the Tibet Autonomous Region.In these regions, the annual average mortality rate ranges from 0.5% to 7.6%.The distribution of annual average mortality rate is extremely positively skewed.The cumulative percentage of grids with annual average mortality rates <0.5%, 1%, 2%, and 5% are 50%, 68%, 80%, and 98%, respectively.At return period of 1/20 a, estimated mortality rates range from 0.09% to 23%.The cumulative percentage of grids with 1/20 a mortality rates <0.5%, 1%, 2%, and 5% are 12%, 24%, 45%, and 89%, respectively.

Risk in terms of livestock mortality
Risk metrics in terms of livestock mortality were then derived by multiplying the mortality rate by exposure (Fig. 7).Again, annual average mortality, and the mortality at 1/10 a, 1/20 a, and 1/50 a were all reported.

Spatial patterns of livestock snow disaster risk in the QTP
Our results illustrate the spatial distribution and offer quantitative metrics of risk in terms of livestock mortality and mortality rate due to snow disasters in the QTP.The spatial pattern of risk agrees with earlier studies covering this region quite well.From an empirical perspective, the literature frequently mentions Easter Inner Mongolia, the Northern Tianshan Mountains in Xinjiang, and Northeastern QTP as centers of snow disaster around China (Gao, 2016;Hao et al., 2002).
Within the QTP, high frequency snow disaster regions that are mentioned repeatedly in the literature include Yushu, Guoluo, Naqu, Shigates, and Nagri (Bai et al., 2011), which have all been identified in our study.As for risk assessment, our results also agree well with earlier studies.For instance, regions between the Kailas Range and neighboring Himalayas, southern Qinghai Province (mainly Yushu and Guoluo), and the northwestern corner of the QTP are all considered as higher risk regions in both qualitative (Liu et al., 2014b) and quantitative (Shi, 2011;p.106-107)risk assessment results.In northern and western Naqu Prefecture, and the central-to-western end of the Nyainqêntanglha Range, our results are consistent with the national snow disaster risk map (Shi, 2011; hereafter termed risk maps), which are of higher-to-the-highest risk.
Nevertheless, these regions are considered the lowest of lower risk in the results presented by Fenggui Liu et al. (2014).
For the magnitude of annual average mortality rate, our results were smaller than those in the risk maps of China (Shi, 2011;p.104-107); in general, our results had values about half those previously reported.For the high-risk regions, annual average mortality rates were generally ≥ 2% in our results, but ≥ 4% in the risk map results.Our result had a vast low risk region with annual average mortality rates <0.5%, but the threshold was 1~3% in the risk maps.In terms of mortality, our results matched historical records better.For instance, the most severe snow disaster in southern Qinghai Province since 1960 was in 1996, mostly in Guoluo and Yushu (Bai et al., 2011), the deadliest disaster in nearly 65 years.It claimed a loss of 1.08 million livestock (equivalent to 1.5~2 10 6 sheep units, assuming the local herd structure for 2016).According to our risk metrics, an annual loss of 2  10 6 (summing Guoluo and Yushu) would have a return period over 100 years (Table 2).
Another example is the 1997-1998 snow disaster in Naqu, the most severe snow disaster since 1960, leading to a loss of 0.82  10 6 livestock (equivalent to have 1.2  10 6 sheep units, assuming the local herd structure for 2016) (Wen, 2008).Again, an annual loss of 1.2  10 6 sheep units has a return period of 50 to 70 years according to our metrics.As we have assumed stronger prevention capacity using the GDP values for 2015, it is natural to have a higher return period than empirical values.
If the mortality rate estimated in risk maps were used instead, then the corresponding return-period could be underestimated by approximately 1/20 a, and consequently, snow disaster risk would be exacerbated.

Advantages of the event-based PRA
Our study differ from the existing literature largely in its event-based PRA framework.Such a framework derives unique information that can hardly be obtained by earlier methods that are based on annual basis, which are important for preparedness decisions and insurance solutions.For instance, the event-based framework enables the estimation of frequency distribution of the occurrence and duration of single disaster events.Overall, our analysis indicates that snow disasters are frequent in terms of annual occurrence but more than one snow disaster a year is unlikely (Fig. 5).Given this finding, prior counter-measures should be implemented to build prevention capacity to handle the one event in every year (Mechler et al., 2010).
Our results for single event duration provide important quantitative references for hay and fodder storage, which were not achieved by earlier annual basis analyses.For the majority of higher-risk regions, once a snow disaster occurred, on average it lasted for 12 d (Fig. 5).At return periods of 1/10 a and 1/20 a, the durations of single events were up to 21 days and 28 days, respectively (Fig. A2).At return periods of 1/50 a, single events could even last for more than 40-50 days.The regional average duration of a 1/20 a event in Naqu, Yushu, Guoluo, and south Ngari, was estimated to be 24, 22, 26, and 26 days, respectively.From a preparedness perspective, the amount of hay and fodder storage needed, from combined herder households and local government reserves, can be readily estimated from our results once their goal of preparedness capacity is set, i.e., capable of managing a 1/10 a event.Alternatively, our results can also help local regions measure their preparedness capacity given their amount of hay and fodder storage.
Our event-based PRA results also provide solid technical support for insurance solutions.Earlier studies that assess risk on annual basis using annual aggregate snow-cover days, or snow depth variables are not capable of doing so as insurance indemnities must be clearly triggered by specific events.The capability of finding the frequency distribution of event occurrence and event duration provides necessary information to help the design of insurance trigger schemes.Meanwhile, our results can readily support the calculation of actuarially fair premium rates and cat-risk loadings by applying deductible conditions to turn event loss records into event-based insurance losses, and then calculating annual average and return period insurance loss-cost ratios (Wang and Zhang, 2003).

Risk-informed implications
Our results imply that the present level of preparedness in local regions are far from sufficient.According to our local survey data, a sheep unit, even under the minimum level of supplementary feeding, consumes 1 kg of hay and 0.5 kg of fodder (Ye et al., 2017b).Nyainrong County in central Naqu Prefecture prepared a 4.4-million yuan fund to prepare winter storage of hay and fodder in 2016.The total amount of hay bought (460 ton) can only support supplementary feeding of county-wide livestock (1167.6 thousand sheep units) for at most 3~5 days.Such a level of preparedness can only endure a snow disaster with return period less than 1/5 a.This discrepancy also explains the frequent losses from snow disasters in these regions.
It is not straightforward to improve such preparation capacity to manage higher return-period events.As the herd size has mostly reached the defined carrying capacity, hay harvests for winter seasons from local grassland are much less feasible (Shang et al., 2012).Agro-pastoral regions, i.e., the northeast, southeast, and southern parts of the QTP, mostly the agropastoral regions in Qinghai, Gansu, Sichuan, and Shigatse in Tibet, can obtain some support from crop straws.Pure pastoral regions, i.e., Naqu and Nagri, are problematic, due to the lack of local resources, and the high transportation cost to import from other regions.To enhance preparedness capacity, inter-prefecture overall planning will be needed to arrange support from agro-pastoral regions, e.g.Shigatse and Lhasa, to pastoral regions.In addition, forage harvesting bases within small parts of pastoral regions with relatively large precipitation and grassland productivity will also be favoured.However, it remains less likely that local regions will be able to provide sufficient hay and fodder to endure long-lasting events, i.e. > 14 d or an equivalent event with a return period larger than 10 a.
Due to the difficulty in improving prevention capacity, insurance schemes are needed to provide relief (Mechler et al., 2010).
The insurance product could be conventional (indemnity-based), where the post-disaster loss-adjustment is conducted on herder household bases.A more favourable choice is to adopt an index-based structure using livestock mortality rate as predicted by our model, in which disaster duration and growing season aggregate precipitation are two critical indices.The index-based structure is much more efficient in vast and sparsely populated regions by expediting the loss-adjustment process, and consequently, should save considerable associated costs (Ye et al., 2017b).Furthermore, using a modelpredicted mortality rate can encourage local communities to invest in prevention and preparedness, and reduce moral hazards (Miranda and Farrin, 2012).

Limitations
Several limitations in our risk assessment model must be mentioned.First, our hazard module to rebuild/predict snow disaster still suffers from uncertainty.We obtained a good AUC score from the BRT model for identifying snow disaster days, and also good agreement in timing of occurrence and distribution of duration for longer-duration events.However, the performance in capturing small disasters of short-duration, i.e., < 5 d, still needs improvement.
Due to the lack of exact spatial distribution of sheep units, the exposure data was derived according to the computed carrying capacity by grassland types.The total herd sizes computed differ within 20% to those reported in the statistical yearbooks at the prefecture levels for pastoral regions.The difference for agro-pastoral prefectures in the eastern QTP was much larger, as agriculture provides extra foodstuff for raising livestock and supports livestock farms in addition to free grazing.Fortunately, as snow disaster threats have little impact on livestock kept on farms, our result is less likely to severe underestimate mortality in these locations.For prefectures in the central and western part of the region, our result is of practical significance for use as a reference for livestock mortality under the long-run forage-livestock balance, irrespective of the present true herder size.
Finally, our risk metrics were derived from events rebuilt from historical climate data, but not from stochastic simulations.
Consequently, we have a limited number of events and annual loss records.We are only confident for risk metrics less than 1/35 a. Metrics for any higher return periods were derived from extrapolation and must be used with caution.This limitation can be resolved by inputting stochastic climate datasets using a stochastic weather simulator.
Fig. 1 Probablistic Risk Assessment Framework for Livestock Snow Disasters in the Qinghai-Tibetan Plateau (QTP) region Fig. 2 Technical flow of the snow disaster event identifier/simulator Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-182Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 6 September 2018 c Author(s) 2018.CC BY 4.0 License.mm/d).Snow cover obtained from satellite imagery data was not considered because it is unavailable for future periods and predicting disasters, although, in this study, we used only historical data.Historical snow disaster event data with the time of each event for each meteorological station were used to train the BRT model.These data were obtained from two sources.Records for 1980-2007 were obtained from W. Wang et al. (2013) while records from 2008-2015 were obtained from the China Meteorological Science Data Sharing Service System (CMSDS, http://data.cma.gov.cn).Data for the predictors were also obtained from CMSDS, including 106 national reference stations in the region.The dataset contains daily observations of maximum, mean and minimum temperature, maximum and mean wind speed, and precipitation.
Nat. Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-182Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 6 September 2018 c Author(s) 2018.CC BY 4.0 License.Finally, the fitted BRT model together with the tuned parameters of smoothing and filtering was applied to generate all snow disaster events during 1980-2015 by grid.The China meteorological forcing dataset (He and Yang, 2011) obtained from the Scientific Data Centre of Cold and Arid Regions http://westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21 was used.It offers variables, including precipitation, air temperature, wind speed, and sunshine duration at spatial resolution of 0.1° × 0.1° and temporal resolution of 3 h.We used this dataset because it focuses on the cold and arid regions in western Nat. Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-182Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 6 September 2018 c Author(s) 2018.CC BY 4.0 License.statisticalyear books of Qinghai, Tibet, Sichuan, Gansu, and Xinjiang.County level GDP values were assigned to each grid within its boundary.We used constant GDP values for 2015 for two reasons.First, the results can be directly treated as a stationary time series for estimating the probability distribution, as the influence of prevention capacity improvement has been removed.Second, it meets the goal of risk assessment, to estimate the likelihood of potential loss in the near future given present-day prevention capacity (prevention capacity equivalent to year 2015).

Fig. 3
Fig. 3 Relative influence of variables predicting a snow-disaster-day.Blue bars are relative importance of each factor, and the sum of relative importance is 100%.

Fig. 4
Fig. 4 Empirical cumulative density functions for historical and model-predicted snow disaster duration Fig. 5 Gridded annual frequency/ annual average occurrence (a) and mean event duration (b) of snow disasters from model predictions Fig. 6 Gridded livestock snow disaster risk in terms of mortality rate (%) in annual average values (a), and return-period values of 1/10a (b), 1/20a (c) and 1/50a (d).The grid size is 0.1° × 0.1°.
Fig. A2 Spatial distribution of livestock exposure estimated from vegetation distribution

annual mortality Risk metrics: expected/ return-period mortality (rate) Carrying capacity Herd size distribution Vegetation distribution Hazard Exposure Vulnerability Event duration Climate data (gridded) Snow disaster event identifier/simulator Timing (dates) Growing season precipitation Vulnerability function (Li et al. 2018; Ye et al. 2018) Event loss modeling (mortality rate)
were considered.Following their suggestion from the multi-method comparison, we chose the predictive version of the generalized additive model, Given such a relationship, the vulnerability is a truly dose-response function between livestock mortality rate (mortality/herd size) and snow hazard intensity together with other environmental stressors and prevention capacity, as proposed by Event/ Editorial Committee of Vegetation Map of China, 2007), which offers detailed information about the spatial distribution of 11 vegetation type groups, 55 vegetation types, 960 plant formations, and more than 2000 dominant species in vector data.To match the look-up table and map information, we merged some vegetation types and used only the major grassland types (percentage area >0.5%) according to the FAO survey (2005) (Fig.
(Zhang et al., 2014)ype to carrying-capacity relationship.For the look-up table, we adapted the plan of Xin et al. (2011) for Qinghai and criteria outlined by the Land Management Administration of Tibet Autonomous Region(1994)for Tibet after reviewing various criteria(Zhang et al., 2014).For grassland distribution, we used the Vegetation Map of the People's Republic of China (1:1 million) (