Estimating Return Intervals for Extreme Climate Conditions Related to Winter Disasters and Livestock Mortality in Mongolia

. Mass livestock mortality events during severe winters, a phenomenon that Mongolians call dzud, cause the 15 country significant socioeconomic problems. Dzud is an example of a compound event, meaning that multiple climatic and social drivers contribute to the risk of occurrence. Existing studies argue that the frequency and intensity of dzud are rising due to the combined effects of climate change and variability, most notably summer drought and severe winter conditions, on top of socioeconomic dynamics such as overgrazing. Summer droughts are a precondition for dzud because scarce grasses cause malnutrition, which in turn makes livestock more vulnerable to harsh winter conditions. However, these studies 20 typically look at a short time frame (i.e., after 1940), and few have investigated either the risk or the recurrence of dzud over a century-scale climate record. This study aims to fill the gaps in technical knowledge about the recurrence probability of dzud by estimating the return periods of relevant climatic variables: summer drought conditions and winter minimum temperature. We divide the country into three regions (Northwest, Southwest, and East Mongolia) based on the mortality index at the soum (county) level. For droughts, our study uses as a proxy the tree-ring reconstructed Palmer Drought 25 Severity Index (PDSI) for three regions between 1700-2013. For winter severity, our study uses observational data of winter minimum temperature after 1901 while inferring winter minimum temperature in Mongolia from instrumental data in Siberia that extends to the early 19th century. Using a Generalized Extreme Value distribution with time-varying parameters we find that the return periods of drought conditions are changing over time, with variability increasing for all the regions. Winter severity, however, is not changing with time. The medians of the 100-year return periods of the winter minimum temperature 30 in Mongolia over the past 300 years are estimated as -26.08°C for the Southwest, -27.99°C for the Northwest, and -25.31°C for the East. The co-occurrence of summer drought and winter severity increases in all the regions in the early 21 st century.

. A better understanding of the climate drivers of dzud and extreme events is also critical for preventive and responsive measures, such as weather index insurance. Weather index insurance recently became widely available, and its indemnities are paid based on realizations of a weather index such as rainfall and temperature that are expected to be highly correlated with actual losses, rather than on actual losses experienced by the policyholder (Barnett 65 . The advantage of index insurance is that the pre-determined index cannot be manipulated by the third parties. Payment is faster than loss-based insurance because payment will be made once the predetermined index exceeds the threshold. In contrast, for loss-based insurance, an insurance company must assess losses before making payment, which requires labor and time. The lower transaction costs of index insurance can also make it a more affordable product for the purchaser and thus a more viable offering for the insurance provider. The index-based livestock insurance program (IBLIP) 70 was institutionalized in 2014 to respond to the extreme climate disasters by the Government of Mongolia with help from the World Bank (Skees and Enkh-Amgalan, 2002;Mahul et al., 2015;Mahul and Skees, 2007). Hessl et al. (2018) analyzed variabilities of drought in Mongolia using long-term climate data. Such studies are still limited by the fact that there are few long-term instrumental records of climate in the region, and the records that do exist are often not continuous and contain missing data. Though historical documents record the occurrence of dzud from the 19 th century, 75 changes in climate in Mongolia have been observed in instrumental records only since 1940 (Batima et al., 2005). Some studies reported that the frequency of dzud has increased, such as after 1950 (Fernandez-Gimenez et al., 2012;Middleton et al., 2015) or after 2000 (Munkhjargal et al., 2020). Furthermore, another study concluded that it is expected to increase with future climatic changes (Bayasgalan et al., 2009) using Coupled Model Intercomparison Project's climate models.
Natsagdorj (2001) shows that the trends of drought and the dzud index, estimated by normalized monthly temperature and 80 precipitation, are increasing. However, these studies are based on observational data of dzud, which are available only from about 1940. It is critical to extend the time horizon in order to improve the reliability of the return period estimation of catastrophic dzud, and of the assertions of secular changes in the causal climate variables. Long-term climate proxies, such as tree rings, have the potential to do so by deriving recurrence periods of dzud and climate extremes, especially to improve index insurance products (Bell et al., 2013). Yet, one of the challenges of improving the reliability of recurrence estimations 85 is the lack of scientific understanding of the historical trends of past climate events due to the short meteorological record (Mahul and Stutley, 2010;Mcsharry, 2014;Rao et al., 2015).
To improve risk analysis of dzud, the investigation of extreme distributions of climate extremes is critical. D'arrigo et al.
(2001) inferred using millennial length tree-ring data that temperatures in Mongolia in the late 1990s and early 2000s were extraordinary. Based on well-calibrated and verified millennial-length tree-ring reconstruction of summer temperatures, Davi 90 et al. (2015) and Davi et al. (2021) show that the recent warming trend since the 1990s is anomalous in the long-term context in Mongolia. In addition, Davi et al. (2010) conducted spectral analysis to discover the periodicity of droughts in Mongolia by using tree-ring based reconstructed Palmer Drought Severity Index (PDSI). However, these studies do not estimate probability distributions of extreme climatic events or improve the reliability of the estimation of return periods of dzud for risk analysis. Here, we use the term "risk analysis" to refer to the analysis of the probability of an extreme event whose 95 consequences could be substantial (Rootzén and Katz, 2013), but not the analysis where risk refers to the combination of the probability of an event and its associated expected losses.

Objectives of the study
The objective of this study is to conduct risk analysis for the climatic variables that cause dzud, namely summer drought followed by extreme cold temperature and snowfall, in Mongolia while attempting to improve the reliability of the return 100 period estimation of dzud utilizing tree-ring proxies and historical data on climatic variables. The study also explores the implications of the risk analysis and return period estimation for index insurance using tree-ring data. To address these objectives, we posed the following research question: • How can the reliability of the return period estimation of climate extremes be improved?

105
There are two important climatic variables to predict dzud: summer drought conditions and winter temperatures (Lall et al., 2016;Rao et al., 2015). Notably, this study estimates return periods of extreme drought conditions, by using tree-ring based reconstructed PDSI from the Monsoon Asia Drought Atlas or MADA (Cook et al., 2010). It also estimates return periods of extreme cold temperatures in Mongolia. Since temperature data in Mongolia is only available from the early-to-mid 20th century, we simulate them from meteorological data in neighboring Siberia, which has records that extend back to the late 110 1800s, through a statistical model presented here. Tree-ring based temperature reconstructions in the region are typically limited to the summer growing season and do not capture winter temperatures. In Mongolia, the term "dzud" refers to high livestock mortality (Fernandez-Gimenez et al., 2012;Morinaga et al., 2003), however, we use climate variables to determine risk rather than mortality because mortality rate assumes that the size of the population does not matter. In fact, changes in livestock populations also matter since they can also be related to changes in 115 socio-economic factors, such as shortage of food supply, which can be related to non-climate factors. Other socio-economic factors also determine livestock herding loss, including the total number of animals and the density per square kilometer.
These numbers drastically increased after a transition to private ownership in 1990s (Douglas A. Johnson, 2006;Rao et al., 2015;Reading et al., 2006). The increased livestock population results in overgrazing and degradation of the grassland, which resulted in a decrease in the grassland carrying capacity and a high mortality rate (Bat-Oyun et al., 2016;Berger et al., 120 2013;Hilker et al., 2014;Liu et al., 2013).
In order to estimate a return period of an extreme climate event, extreme value theory (EVT) is useful (Cheng et al., 2014;Katz et al., 2002). There are ongoing debates about which methods are most suitable for estimating extremes, such as the return period and expected number of exceedance (Read and Vogel, 2015;Rootzén and Katz, 2013;Salas and Obeysekera, 2014). EVT informs us how to extrapolate a rare event which has not been experienced for a long time from existing 125 observational data with a short record. EVT is a widely used method for estimating the probability of extreme hydroclimatic events (Katz, 2010;Leonard et al., 2014;Slater et al., 2021), such as floods (Prosdocimi et al., 2015;Willner et al., 2018) (Willner et al., 2018), precipitation (Gao et al., 2018;Minářová et al., 2017), and compound events (Leonard et al., 2014). EVT helps us to formulate a risk management strategy by deriving a distribution of extreme climate events and estimating a possible extreme value for the future's preparedness. There are two main approaches in EVT: The block maximum approach 130 and the threshold approach, which will be described in Data and Methodology. The objectives of this study were to; • 1. Estimate return periods of extreme drought conditions by using reconstructed PDSI .
• 2. Estimate return periods of extreme cold temperatures in Mongolia by using long instrumental data from Siberia.
Conventionally, a stationarity process is assumed to estimate return periods,. Here, we consider the extension of the record 135 by explicit dependence on climate proxies. We examine how the return periods may change over time due to slowly and systematically changing climate conditions, persistence in the PDSI, or other climate records. Exploring the nonstationary approach to return period and risk opens "many opportunities" (Salas & Obeysekera, 2014). This has the advantage of reducing the bias in the near-term projection, assessment of the return period, and recurrence interval associated with the event. 140 We also explore the utility of using long-term climate proxies in the context of index insurance. In general, the index used for index insurance must be scientifically objective and easily measurable. The Index-Based Livestock Insurance Program (IBLIP) in Mongolia uses mortality rate as the index. Our study provides insights into the long term variations in the mortality rate due to climate with the goal of reducing the bias and the variance of the estimates of the probability of the index used, by identifying the trend, and using a longer record, respectively. 145 2 Data and Methodology

Data
PDSI is a standardized index that ranges from -10 (dry) and +10 (wet) based on a water balance model, accounting for 150 precipitation, evaporation, and soil moisture storage (Cook et al., 2010;Dai et al., 2004;Palmer, 1965). In this study, treering reconstructed PDSI values from 1700 to 2013 are taken from Monsoon Asia Drought Atlas (MADA) (Cook et al., 2010). MADA is a seasonally resolved gridded spatial reconstruction of drought and pluvials in monsoon Asia over the last 700 years, derived from a network of tree-ring chronologies (Cook et al., 2010). The MADA can also reveal the occurrence and severity of previously unknown monsoon droughts (Cook et al., 2010). We consider three regions (Northwest,155 Southwest, and East Mongolia) in Mongolia, as in Figure 1, based on clusters proposed by the previous studies Lall et al., 2016). These spatial clusters are based on the mortality data at the soum (county) level from 1972 to 2010, using hierarchical clustering (Johnson, 1967), which were adjusted with the spatial patterns of the Mongolian topography, climate zones, and mean precipitation in growing seasons. It is reasonable to use these clusters because the objective of the study is to improve risk analysis of Dzud and mortality of livestock in Mongolia. 160

Preliminary Analysis
The correlation in PDSI values from 1700 to 2013 between three clusters is shown in Table 1. The Mann-Kendall trend test is used to examine the trends of the PDSI data (Kendall, 1948;Mann, 1945). The Mann-Kendall test shows that there are no monotonic trends in the PDSI data for all clusters (Table 1). Yet, times series of tree-ring reconstructed PDSI by clusters show that there is significant centennial-scale variability, which is important to consider since they suggest that there are 170 persistent regimes that can last for decades to centennial time scales (Figure 2 (a)-(c)). Though these may occur randomly or reflect systematic cyclical behavior, their consideration in a risk management strategy is critical.  The autocorrelation function (ACF) and Partial ACF of all the regions show that there are significant autocorrelations in the 180 PDSI data in all clusters (in Figure S1 and S2). The development of a time series simulation model that uses these long lead correlations would help inform the risk analysis associated with the persistent regimes identified earlier. Autoregressive Integrated Moving Average (ARIMA) models with different orders were evaluated based on the Bayesian Information Criterion (BIC), which can account for fitting errors for the Bayesian conditional mechanism of models. Please note that the BIC is standard information-theoretic criteria whose relative magnitudes allow one to choose one model over another 185 (Akaike, 1979;Burnham and Anderson, 2004). The order of the best ARIMA models in each cluster is (3,0,0) for the Southwest, (1,0,2) for the Northwest, and (1,0,0) for the East. These ARIMA models will be used later to forecast the effective return periods of droughts.

Climate variables
Models that use climate variables as covariates are explored for developing a nonstationary risk model. These data are 190 summarized in Table 2. We use high-resolution gridded datasets at Climate Research Unit (CRU) at University of East Anglia for monthly temperature, and summer (May-August) and winter (November-February) precipitation for the three clusters (Harris et al., 2014). All the gridded points within each cluster are averaged. We also used average monthly temperature data from instrumental records in Siberia, including Irkutsk (1882), Minusinsk (1886), and Ulan Ude (1895-1989. We also use the Arctic Oscillation (AO) index, which comes from two sources: the Joint Institute for the 195 Study of the Atmosphere and Ocean (JISAO) and the National Oceanic and Atmospheric Administration (NOAA). The two records were merged into one record (e.g. Kaheil and Lall, 2011)). The AO index is closely associated with summer and winter climates in East Asia (He et al., 2017). In particular, the negative phase of AO is associated with more frequent cold air outbreaks in East Asia, including Mongolia (Cohen et al., 2010;He et al., 2017;Yu et al., 2015). Finally, please note that though dry conditions of PDSI is negative, all the analyzed PDSI values below are presented in reversed values because the 200 used R package, extReme (Gilleland and Katz, 2016), will capture the maximum values.

Methodology 205
Extreme Value Analysis (EVA) is utilized in this study. In EVA, the distribution of many variables can be stabilized so that their extreme values asymptotically follow specific distribution functions (Coles et al., 2001). There are two primary ways to analyze extreme data. The first approach, the so-called block maxima approach, reduces the data by taking maxima of long blocks data, such as annual maxima (Coles et al., 2001). The Generalized Extreme Value (GEV) distribution function is fit to maxima of block data, as given by 210 (1) 2 where, + = max{ , 0} , > 0, − ∞ < , < ∞.
Equation (1) covers three types of distribution functions depending on the sign of the shape parameter  The Fréchet distribution function is for  while the upper bounded Weibull distribution function is for  (Gilleland and Katz, 2016).
The Gumbel type is obtained in the limit as →0, which results in The second approach, the so-called threshold excess approach, is to analyze excesses over a high threshold (Coles et al., 2001). The Generalized Pareto Distribution (GPD) has a theoretical justification for fitting to the threshold excess approach (Gilleland and Katz, 2016), as given by where μ is a high threshold, x>μ, scale parameter >0 and shape parameter −∞ < < ∞ . The shape parameter determines three types of distribution functions: heavy-tailed Pareto when >0, upper bounded Beta when <0, and the 220 exponential is obtained by taking the limit as → 0, which gives The extreme value models can be applied in the presence of temporal dependence (Coles et al., 2001), as given below: By examining the times series of the PDSI values and winter minimum temperature, we can enhance the understanding of how return periods of droughts, and extreme cold weather have changed over time. The best GEV and GPD models are selected based on Maximum Likelihood Estimation (MLE) and BIC (Katz, 2013). Also, diagnostic plots were used to assess whether the best GEV and GPD models are reasonably fits to the distributions . Stationarity is assessed through the comparison of the BIC applied to a set of candidate models formualted using equation 5, with terms that include time or not.

Main Results and Discussion
The procedure is implemented as follows: 235 1. Fit GEV distributions to the tree-ring reconstructed PDSI values, allowing for non-stationarity by making µ, σ, and/ or ε a function of time.
2. Fit GEV distributions to the tree-ring reconstructed PDSI values using climate variables (AO index, summer precipitation, snow, and minimum temperatures).
3. Evaluate models based on BIC. 240 4. Using the best GEV model, return periods are estimated. 5. The above procedure is repeated for GPDs fit to the tree-ring reconstructed PDSI values.

Fitting GEV to the Tree-Ring Reconstructed PDSI for Return Period Estimation
We construct two types of models: (1) stationary and nonstationary extreme value models, and (2) nonstationary models using climatic variables as covariates. First, we consider polynomial models in time of the order of 0 to 2 for both the 245 location and scale parameters of the GEV distribution, resulting in seven models to be tested, including the stationary model, for each region. In addition, autoregressive (AR) models are examined. The models are evaluated based on the BIC (Table   3). The best GEV models and its maximum likelihood estimates (MLE) with 95% confidence intervals are as follows (Table   3, Figure 3):   (3) AR (1) Note: L stands for the location parameters, S stands for the scale parameters. 0 means a constant in the parameter, 1 is temporally linear, and 2 is temporally quadratic for each parameter. These results suggest that in the long run, a stationary model for PDSI in Mongolia may be appropriate. Only the Southwest has nonstationarity in the scale parameter. This nonstationarity in the scale parameter for the Southwest, with a mean 280 coefficient of 0.002 relative to the constant value of 0.95, means that over 100 years the variability could increase from 0.9 to 1.05. If we take 0.9 to be a mid-period estimate, this would be rather a modest change. This could be a real feature or an artifact of the non-constant reconstruction variance from the tree ring reconstruction algorithm.
Next, we estimate parameters of the GEV distribution functions fit to the PDSI values by including other climate variables 285 such as AO index, summer precipitation, snow, and minimum temperatures as covariates from 1903 to 2010. Summer precipitation is a mean of May to August of a previous year, while snow is mean of values from November of a previous year to February of the year. Also, AO index data starts in 1903. Thus, we use data starting 1903 though the data itself exists since 1901. The minimum temperature is a minimum value from November of a previous year to October of the year. The GEV models with the lowest BIC for each cluster and MLEs with the 95% confidence intervals are as follows (Table 4 and 290 In the GEV models, climate variables (precipitation and snow) are important covariates for the extreme values of the PDSI values and improve the model performance (Table 3). These climate variables have no inter-year dependence that is 300 significant based on ARIMA, and hence there is no memory in these variables and the best model is stationary model.
Consequently, no near-term forecast is feasible. The time series of effective return periods of 100-year events for the GEV distribution functions fitted to the PDSI using the 310 climate variables are shown in the Southwest, Northwest, and East from 1903 to 2010 ( Figure 5). This shows that variabilities of return periods of 100-year events of the PDSI values become larger over time in all the regions. Before 1940, the variabilities are small possibly because the instrumental data records began in 1940's. Even after 1940's, it also shows that the magnitude of 100-year events has increased in the last half of the data series. A PDSI value of 3 used to be a 100 year event around 1920. Yet, around the beginning of the 21 st century, it has increased to be between 4 and 5. However, 315 considerable inter-annual and decadal variability is evident.

Variabilities of return periods of 100-year events of the PDSI values become larger over time in all the regions. The blue horizontal line is the mean of the effective return periods while the red one is its median. Please note that the vertical axis is shown by the reversed values of PDSI values, meaning that a positive value is a drought condition.
The relationship between significant climate covariates and reversed reconstructed PDSI values based on the best GEV 325 models for each return period of 10, 50, and 100 years events are shown in Figure 6 and Figure 7. This shows that less precipitation leads to higher reversed reconstructed PDSI values, meaning more likelihood of droughts. Consequently, with this model, future projections of precipitation could be helpful to predict drought severity and frequency.  is precipitation, the y-axis is snow, and the z-axis is reversed reconstructed PDSI values. The left? cube is for 10-year events, the central is for 50-year events, and the right is for 100-year events. Since the best GEV model contains 340 precipitation and snow as covariates, the model for the Northwest is cubic. This shows that less precipitation leads to higher reversed reconstructed PDSI values, meaning more likelihood of droughts.

Fitting GPD to the Tree-Ring Reconstructed PDSI for the Return Periods Estimation
To fit a GPD, a threshold needs to be selected. We selected a threshold of 1.0 (please see the Supplement S1 for the detailed explanation of how we chose the threshold). GPDs are fit to the tree-ring reconstructed PDSI values from 1700 C.E. as both 345 stationary and non-stationary models ( Table 5). The model of stationarity is best in terms of BIC for all clusters. Being similar to the GEV cases, we analyze the other climate variables after 1903. Table 6 shows that the best model of GPD 355 is the one with a constant in the scale parameters in terms of BIC for all clusters. MLEs estimated by the best GPD models are shown in Figure 8. The table shows that for catastrophic droughts, climate variables are not a significant covariate, although the differences in BIC values in the Southwest and Northwest between the ones with constants and with AO index are small. The estimated effective return periods based on these best GPD models are listed in Table 7. In Table 7, the difference between the values for 10, 50, and 100 years is slight because the shape parameters estimated from the GEV for 360 each case are negative. This means that the data is negatively skewed and this leads to an implicit upper bound for the process. As a result, each of the quantiles is restricted by that upper bound and ends up quite close to each other.

Results Based on GEV and GPD Models
In this section, we fitted the GEV and GPD distribution functions to the PDSI values. Results are the following: • The results show that the PDSI values will follow the distributions with ε<0, namely the Weibull distribution for the GEV models and the upper-bounded Beta distribution for the GPD models. This information can be used to estimate return periods of extreme drought. 375 • For the Southwest, the non-stationary models performed better if we look at GEV without a threshold. However, with a threshold of 1 for the GPDs, the stationary models perform better than the non-stationary models, which indicate that all trends in reconstructed PDSI values are influenced by small events, not by extreme events; i.e. extreme events are stationary. For both the Northwest and East, stationary models performed better for both the GEV and GPD models. 380 • Compared to the models with constants in the parameters, the GEV model with the climate variables are better in terms of the BIC value. Therefore, establishing a relationship between drought conditions and climate variables, particularly precipitation and snow, is useful in understanding the dynamics that determine dry conditions. However, compared to the models with constants in the scale parameters, the GPD models with the climate variables don't lead to the improvement of the model performance in terms of the BIC criteria. 385 • In terms of BIC, the models of a GPD fitted to tree-ring reconstructed PDSI values show better performance than the GEV models.
• Because of the third point, the effective return periods based on the GEV models change with the climate variables.
In contrast, the effective return periods based on the GPD models are constant: for example, a 100-year event is the PDSI value of -4.17 for the Southwest, -4.76 for the Northwest, and -3.87 for the East. This suggests that while the 390 magnitude of the annual maxima seems to change with the base climate conditions, the frequency of the extreme events beyond a threshold is not affected that much.

Simulating Annual Minimum Temperature in Mongolia Using Siberia Data
Instrumental winter temperature data in Mongolia is limited before 1950. Also the gridded climate database that cover Mongolia starts after 1901. Thus, we attempt to estimate the Mongolia data from longer records from Siberia, which can get 395 back to 1820. Existing studies suggest the winter temperatures between Mongolia and Siberia are correlated spatially, driven by polar jet dynamics (He et al., 2017;Iijima and Hori, 2018;Munkhjargal et al., 2020). First, winter temperatures in Mongolia will be simulated by using instrumental temperature data from Siberia (in Section 3.2). By using the simulated winter temperature in Mongolia, return periods of extreme cold temperature during winters will be estimated in Section 3.3.
The procedure was implemented as follows: 400 1. Conduct correlation analysis between Siberia and Mongolia data to select which station data are informative for temperature in Mongolia.
2. Impute missing data of instrumental data in Siberia 3. Fit a GEV and GPD to the winter minimum temperature in Mongolia with the Siberia data 4. Simulate winter minimum temperature of Mongolia from Siberia data based on the best GEV model. 405 5. Calculate effective return periods of 10, 50, and 100 years from the simulated winter minimum temperature of Mongolia. First, correlation analysis is conducted to see which station data in Siberia is useful for Mongolia data. Temperature data in both Mongolia and Siberia is monthly data. We use minimum temperature and average temperature during the winter time (October to April) to remove the seasonality. We remove the seasonality because if both series data have seasonality (or 410 periodicity), the correlation between them will be high just for that reason. Removing seasonality helps us identify if the anomalies from the periodic behavior are correlated, or namely if they share similar dynamics in effects induced by atmospheric circulation beyond the seasonal cycle. Data are taken for the common periods when all the points have data (i.e. between 1901 -1990). Irkutsk data alone is used since it alone shows significant correlations between the temperature data in Mongolia and Siberia (Results of Pearson and Spearman correlation coefficients and scatter plots are shown respectively 415 in Table S Next, we checked the structures of missing data from Irkutsk. Some years are missing all monthly records. We impute 420 Irkutsk's data with pattern matching methods, which is equivalent to k-nearest neighbors, by Gibbs sampling using predictive mean matching method (Van Buuren and Groothuis-Oudshoorn, 2011). Using winter minimum temperature from the Irkutsk data in Siberia (TminIrkutsk) as a covariate, we fit the Mongolia winter minimum temperature (Tminmongolia) based on the GEV and GPD models.

Fitting GEV to the Winter Minimum Temperature in Mongolia 425
The results for GEV models based on BIC are shown in Table 8. Models with Siberia data both in the location and scale parameter are the lowest BIC for the Northwest. For the Southwest and East, the one with Siberia data in the location parameter and constant in the scale parameter shows the lowest BIC (Table 8). The best models for each region are shown in Figure 9 and in the following:

Fitting GPD to the Winter Minimum Temperature in Mongolia
For GPD, we select 23 (-23 degrees in reality) as a threshold (Please see Section S.1, Figure S6 and Figure S.7 in 445 Supplemental regarding how the thresholds were selected). In this case, the one with the Irkutsk's data in the scale parameter has the lowest BICs for all clusters as Table 9 shows.

Results based on GEV and GPD models
In this section, we fitted the GEV and GPD distribution functions to the winter minimum temperature in Mongolia. The results are as follows: • All the results show that the winter minimum temperature will follow the distributions with ε<0, namely the Weibull distribution for GEV and the upper-bounded Beta distribution. 460 • Based on BIC, GPD models show better performance in both Southwest and East regions, while the GED models show better performance in Northwest.

Return Periods of the Winter Minimum Temperature in Mongolia Simulated from Siberia Data
Next, we simulate the Mongolia winter minimum temperature based on data from Irkutsk Siberia for 197 years using the parameters estimated by the best GEV model. Then, using this simulated Mongolia winter minimum temperatures, we 465 estimate the 90% confidence intervals of return periods of 10, 50 and 100 year events for each cluster (Figure 11). The median of 100 year return periods are -26. 08, -27.99, and -25.31 Celsius degrees for the Southwest, Northwest, and East.
The variations in these density plots come from both the statistical properties of the Mongolia data itself and variations attributed to Siberia data.

Binary Index
Based on the thresholds used in the GPD approaches, we also explored if the frequency of the co-occurrence of summer drought and cold winter temperatures has changed over time. First, we counted cases as a binary value of 1 when both summer drought and cold winter temperatures in Mongolia are below thresholds (-1 for PDSI values and -23 degree for winter temperatures), otherwise zero. Then, we fitted the local binomial and Poisson regressions while computing 480 generalized cross-validation statistic to determine the smoothing parameter (Loader, 2006) (please see Section S.2 in the supplement for used alpha values for each cluster). Figure 12 shows Northwest and East have similar long-term trends: decreasing trends in the early 20th century and the slightly increasing trends since 1990s. The Southwest shows an increase in the mid-20th century and has started happening again in the early 21st century.

Conclusions
Meteorological data in Mongolia is limited in length with many missing values. Therefore, we utilize longer records from paleoclimate proxy data and meteorological data from neighboring Siberia. The motivation was to improve risk estimation for dzud in Mongolia. Based on extreme value theory, this study derives fitted distributions for drought and winter extreme 495 cold conditions. The study also improves the estimation of return periods of extreme drought conditions and winter temperature, using tree-ring reconstructed self-calibrated PDSI, and station records from Mongolia and Siberia.
GEV models without a threshold show that there is a trend in tree-ring reconstructed PDSI data in the Southwest, while there is a no trend in PDSI in both the Northwest and East. However, the threshold approach indicates that extreme events in reconstructed PDSI values are stationary, indicating that catastrophic drought conditions were stationary for the last 300 500 years. si The study estimated the extreme distributions of drought and winter minimum temperatures in Mongolia. The PDSI values follow the distributions with ε<0, namely the Weibull distribution for the GEV models and the upper-bounded Beta distribution for the GPD models. Also, the results of the study show that the winter minimum temperature follow the distributions with ε<0, namely the Weibull distribution for GEV and the upper-bounded Beta distribution. These estimated 505 distributions can be used to improve the risk calculations for livestock index insurance in Mongolia.
Based on the results of our GEV fitted to the PDSI values, we show that climate variables, such as precipitation and snow, are important covariates for the extreme values of the reconstructed PDSI values. However, for catastrophic drought events, climate variables are not significant covariates based on the results of the GPD model fitted to the PDSI values.
The GEV model also shows that the return periods of drought conditions are changing over time and variability is increasing 510 for all the regions. Yet, based on GPD, the return periods of drought conditions are constant: for example, the actual values of the PDSI for the 100-year events are: -4.17 for the Southwest, -4.76 for the Northwest, and -3.87 for the East. The median of 100-year return periods of the winter minimum temperature in Mongolia is -26.08 Celsius degrees for the Southwest, -27.99 Celsius degrees for the Northwest, and -25.31 Celsius degrees for the East.
This study improves the return period estimation of droughts and winter minimum temperature. Summer drought and winter 515 temperature are important predictors for livestock mortality since they explain 48.4% of the total variability in the mortality data, along with summer precipitation and summer potential evapotranspiration . Therefore, this long-term estimation of return periods of these significant predictors can be used to improve risk analysis of high livestock mortality in order to prepare for the winter catastrophes through early warning systems and index insurance.
A binary index for the co-occurrence of threshold exceedance of drought severity and temperature was developed and its 520 temporal variation assessed. The index shows that all the regions have increasing trends of these co-occurrence. Begzsuren et al. (2004) identify that mortality rates are highest in combined drought and dzuds years than those with dzuds or drought alone while examining the co-occurrence of these extreme events with 51 years of observational data. This implies that the increasing trends of the co-occurrence would pose severe socioeconomic impacts on the country's livestock industry. (Rao et al., 2015) 525 Our study estimates the return intervals and underlying probabilistic characteristics of the climate variables. Index insurance requires a proper threshold and the understanding of underlying distributions of risk events. Thus, the estimation of extreme value distributions and return periods has the potential to improve livestock index insurance, which is implemented in Mongolia by the Government of Mongolia with the help of the World Bank (Mahul et al., 2015). Insurance is priced by considering the uncertainty associated with the estimation of the probability of exceeding the threshold at which the pay-out 530 occurs. The estimation of the uncertainty is reduced as the length of record (in our case from the paleoclimate extension) increases. At present, no one in the industry is using paleoclimatic information to extend and reduce coverage costs, but there is interest in using it to understand the clustering of pay-outs. Furthermore, the results of this study increase understanding of how extreme climatic events in arid regions, which are sensitive to anthropogenic climate change, are changing. The urgent needs to improve resilience of the society to this winter disaster is even more unequivocal. 535

Data Availability Statement
All relevant and publicly available data will be shared via a public data repository if the paper is accepted for publication; data sources are clearly specified throughout the paper.

Author contributions 540
MH, ND, MW, and UL designed the research. MH conducted data analysis. MR and CL prepared the dataset. MH prepared the manuscript with contributions and feedback from all co-authors.