2.1. Meteorological data

Geostatistical modelling allows data to be distributed from a set of weather stations to a grid. Detrended kriging and inverse distance weighting are two methods widely used to interpolate standard parameters from meteorological networks (see Reference Gonzales and CrampeGyalistras, 2003; Reference GyalistrasHaberlandt, 2007). While air-temperature distributions often feature an approximately constant lapse rate, this is certainly not the case for snow-depth gradients. Therefore, a special detrending algorithm was developed to cope with the typical nature of vertical snow distributions (Reference Foppa, Stoffel and MeisterFoppa and others, 2007).

Geostatistical detrended interpolation methods can only produce accurate results if the spatial density of input data is high enough to reflect (1) elevation gradients and (2) regional trends of the respective parameter. We thus allocated data from four monitoring networks, from 150 stations (Fig. 1), for six meteorological parameters and for 13 years: Two snow-monitoring networks run by the Swiss Federal Institute for Snow and Avalanche Research (SLF) supplied data on snow depth (HS) and daily new snow sum (HN). Additionally, the Swiss Federal Office of Meteorology and Climatology (MeteoSwiss) provided data on air temperature (TA), precipitation (PREC), wind speed (WS) and relative sunshine duration (SUN) from two further measuring networks. Most of the 150 stations featured data for only two or three out of the six parameters considered here. A data record (i.e. a 13 year series of daily data for one parameter at one station) was only considered if less than 2 years of data were missing/faulty. To ensure homogeneous data quality, all input data were thoroughly checked for plausibility, cropped where necessary and corrected where possible. Altogether, 286 data records were available for this study (Table 1).

Table 1.

Availability of snow and weather data: number of stations that provided a 13 year data record for a specific meteorological parameter. Records were not considered if more than 2 out of 13 years of data were missing/faulty

Fig. 1.

Canton Grisons, Switzerland: location of snow and weather stations (open circles) and territories of the 19 chamois populations (bordered by black lines).

HS: Snow depth; HN: new snow sum; TA: air temperature; PREC: precipitation sum; WS: wind speed; SUN: relative sunshine duration.

1. ¹Measuring network run by SLF.

2. ²Measuring network run by MeteoSwiss.

The snow and meteorological data were interpolated to grid using WinMet (http://www.hydrologie.unibe.ch/PREVAH; Reference Visscher, Merrill, Fortin and FrairViviroli and others, 2007), where for each parameter the methodology was optimized specifically. This procedure resulted in six grids per day (i.e. 13×365×6 = 28 470 grids), featuring the spatial distribution of HS, HN, TA, PREC, WS and SUN with 500m gridcell spacing. Note that geostatistical detrended interpolation mainly reflects elevational gradients and spatial pattern resolved by the station network. But the method does not correct for aspect and slope.

At this stage, we added an avalanche risk index as a seventh potential environmental driver. The spatial distribution of this parameter required a different approach. The daily Swiss avalanche danger bulletin (http://www.slf.ch) provided the basis for calculating a site-specific avalanche risk index. Based on topographic adjustment procedures (see Reference Johnson, Parker, Heard and GilinghamLeuthold and others, 1997), the regional danger index was modified by means of site-specific weighting factors. These weighting factors quantitatively took into account how well aspect, slope and altitude matched the dangerous areas as specified in the bulletin.

2.2. Mortality data

For management purposes, the Fish and Game Department of canton Grisons (AJF) has systematically collected data on chamois mortality since 1992. All carcasses found, whether by the professional game wardens on patrol or by hunters or hikers, are to be reported to the AJF. The carcasses (including incomplete leftovers) are evaluated by trained staff according to a fixed protocol. Besides place and date of discovery, for each carcass found sex, age, cause of death and date of death are determined or estimated wherever possible.

Determination of the date of death was especially relevant for this study. Sometimes carcasses were only found after months, especially after melt-out and again in later summer at the beginning of periods of high activity of hikers/ hunters (preparing for the autumnal hunt). Roughly 45% of the carcasses were older then 2 weeks. Based on finding date and carcass age (in a set of categories), we recalculated the date of death, albeit with decreasing accuracy for older carcasses. For very old carcasses the assessment distinguished between ‘died in summer’ and ‘died last winter’, thus allowing each fatality to be attributed to at least a specific summer/winter season.

Causes of death were recorded in 34 categories, which we aggregated to the six main categories listed in Table 2. In the following we only consider natural mortality and thus disregard the categories ‘Shot’ and ‘Road kills’. This procedure yielded a total of 6500 data records of chamois that died of natural causes over the course of the study period.

Table 2.

Incidence of different causes of death

The chamois inhabited 19 administrative units, which have been chosen to reflect self-contained populations (Fig. 1). The boundaries between these units typically comprise major rivers, highways and major roads, while part of the outer rim was defined by rather arbitrary borders. These boundaries may not be completely insuperable for chamois, but they are useful for considering populations as being geographically separated, at least to a certain degree.

To compare mortality rates between populations (see section 3.3), the mortality counts had to be normalized by the respective population sizes (hereafter referred to as relative winter mortality). There were only rough data available on annual variations in population sizes. Instead we used one accurate population-size estimate per population for the whole period and accounted for interannual effects by a random factor in our mixed-effects model (see section 2.4). Population-size estimates were reconstructed based on hunting and carcass data over 13 years. Given the hunting level vs natural mortality rates, and the fact that all hunted chamois were analyzed for age and sex, these reconstructions were presumably quite accurate.

¹Only these categories were counted as natural mortality.

²This category contains the sub-category ‘avalanche accidents’ (see Fig. 4).

³Outside the regular hunt.

Note that the populations are managed, whereas hunters may take down about 4000 chamois annually in autumn. The sex-, age- and population-specific hunting quotas were adapted each year according to actual density and age structures. Thus, population structures were rather stabilized by management, and density effects should not greatly confound our analysis of aggregated winter mortality data. Since we see no reason to suppose a different average presence/activity of game wardens and hunters between years, we assume that the probability of detecting a carcass is approximately constant on a year-to-year basis. However, the seasonal statistics may be biased in this respect. Large predators (wolf/lynx/bear) were largely absent in the study area, with only occasional occurrences of wolf or lynx. Given that predation was only an insignificant contribution to the overall mortality (Table 2), variability in predator abundance should not confound our results.

2.3. Data aggregation

For a first overview of the mortality data, the records were summarized in monthly statistics (see section 3.1). The temporal classification was based on the date of death, i.e. not on the date when the carcass was found. In this context, we could only include mortality records that allowed the date of death to be determined within an accuracy of at least a few weeks (83% of data).

For the further analysis, we concentrated on the variability of winter mortality, which enabled season-specific effects of environmental drivers to be assessed. We did not consider environmental drivers with a time lag of >1 year, since the hunting regime (see section 2.2) was expected to equalize long-term effects. The analysis required merging the mortality data (single case records) with the snow and weather data (daily meteorological grids). For this purpose, we first pooled the mortality records in 19×13 data fields containing per population per winter season the cumulated number of mortalities with a date of death between 15 December and 15 June. Then each of the 19×13 data fields was complemented with specific snow and weather data: per winter season and meteorological parameter we produced a set of seasonal grids by aggregating the daily grids over five concurrent or preceding periods (Fig. 2).

Fig. 2.

Scheme for temporal aggregation of data. The number of mortalities with a 15 December–15 June date of death is referred to as winter mortality. Associated snow and meteorological data are temporally aggregated in concurrent or preceding periods: melt-out season (15 April–14 June), snow-free season (15 June–14 October), snow-up season (15 October–14 December) and snow-covered season (15 December–14 April).

Note: The following acronyms for aggregation methods are used in the data analysis: mean: mean value; sum: seasonal sum; days5†: number of days with HN > 5 cm, WS > 5ms^–1; days4.5^†: number of days with ARI > 4.5; days–10^†: number of days with TA < –10˚C; days25^†: number of days with HS/HN > 25 cm.

¹ARI is a quasi-continuous variable with 66 possible levels; ARI > 4.5 can be interpreted as high risk.

²Averaged/summed over snow-covered and melt-out season.

^*Variable x was transformed using log(x). ^†Variable x was transformed using sqrt(x). ^‡Variable x was transformed using asin(sqrt(x)).

For the spatial aggregation of these meteorological grids, the AJF provided data on population-specific seasonal habitats (known traditional winter and summer habitats). These allowed the snow and weather data to be evaluated for each population specifically within their expected seasonal residence areas. We used the summer habitats as aggregation domain for melt-out and snow-free seasons, and the winter habitats as aggregation domain for snow-up and snow-covered seasons respectively.

As aggregation methods we used averaging, summing and conditional counts (Table 3). Ultimately, the procedures described above resulted in 19×13 data fields containing per population and winter season the number of mortalities, and 58 associated population, habitat and season-specific environmental indices (see respective entries in Table 3).

Table 3.

Variables used in mortality statistics: units, range, methods used in seasonal spatio-temporal aggregation, and transformations

Finally, geomorphological habitat characteristics were included in the dataset as additional explanatory variables for analyzing the variability of mortality between populations (see section 3.3). We therefore evaluated elevation, slope and aspect of the population-specific seasonal habitats using a digital elevation grid of 100m horizontal resolution.

2.4. Statistical methods

As a matter of principle, the mortality of chamois displayed variability between years and between populations. The changes in mortality counts between years may reflect the interannual variability of environmental factors. On the other hand, steady differences of relative mortality rates between populations may, among other factors, relate to persistent regional climate gradients. That is why we analyzed between-years effects and between-populations effects of snow and weather separately.

We used linear mixed-effects models with population as random factor in order to investigate between-years effects. With this method the influence of intrinsic factors on mortality was accounted for by the random effect (population). The extrinsic fixed effects (snow and weather) were then tested against the residuals. Between-population effects were analyzed the same way, except that in this case we used year as the random factor. Response variables and explanatory variables were transformed where appropriate using standard filters (Table 3) to meet the assumption of normality and homoscedasticity.

The selection of a best-fit model followed stepwise backward–forward selection methods in S-PLUS (Insightful Corporation, Seattle, WA, USA) based on Akaike’s information criterion (AIC; Reference AkaikeAkaike, 1981). However, where problematic we suppressed the presence of strongly correlated explanatory variables in the best-fit models by manually removing one of two correlated variables and rerunning the model selection procedure.

3.1. Mortality and season

Averaged with respect to year and populations, the mortality data showed a clear seasonal pattern (Fig. 3). Mortality peaked around February to March, when typically snow depths were highest (Fig. 4). The intra-annual variability was dominated by fatalities related to poor general conditions and diseases. This category represented 49% of all cases and 63% of the cases with known cause of death respectively. Accident-related deaths (fall, avalanches) tended to occur towards the end of the winter in March to May. Ten per cent of all fatalities were attributed to avalanches. Predation rates were marginal.

Fig. 3.

Natural mortality of chamois by month, sex and cause of death. Carcasses that could not be attributed to a sex were counted as 0.5 male and 0.5 female. The temporal classification is based on date of death.

Fig. 4.

Snow depth, new snow sum, avalanche risk index, and fatalities due to avalanches (see Table 2) in terms of months. All data scaled to arbitrary units.

Note that the data shown in Figures 3 and 4 only represent 83% of all records, since an approximate date of death could not be determined for 17% of the carcasses. As most of the undatable carcasses were found around melt-out in May, many of these cases may represent fatalities of the snow-up period. This period might thus be slightly underrepresented in the mortality statistics (Fig. 3). Note also that the mortality counts in August and September are presumably biased due to the autumnal hunt. The activity of hunters is likely to increase the carcass detection probability and might even induce additional (natural) fatalities due to stress. However, elevated mortality counts in August and September also reflect occurrences of Mycoplasma conjunctivae, which cause additional fatalities primarily during these 2 months (Reference GiacomettiGiacometti and others, 2002).

¹Calculated according to Reference WiersemaXu (2003), R ² _fix+random represents the variance of the response variable explained by fixed and random effects, while R ² _fix represents the variance explained by fixed effects alone.

3.2. Winter mortality and weather

Between-years effects consider the influence of environmental drivers on the interannual variability of mortality, i.e. we determine which snow and weather conditions entail years with increased numbers of fatalities. In the following, we concentrate on the variability of winter mortality, which enables timing issues between environmental drivers and mortality as response to be assessed.

The variance analysis revealed that winter mortality counts of chamois were predominantly associated with mean snow depths during the snow-covered period (Table 4), as also reported in previous studies (e.g. Gonzalez and Crampe, 2001). The avalanche risk index during the snow-up season was identified as a secondary predictor for mortality. It is remarkable that neither snow depth during melt-out nor any measure of the winter length appeared in the best-fit model. Hence, spring weather conditions seem less critical. Moreover, we did not find any association between winter mortality and our explanatory variables from the preceding snow-free season.

Table 4.

Linear mixed-effects model for the interannual variability of the mortality of Alpine chamois. 71% of the variance is accounted for by random and fixed effects together, and 36% by the fixed effects alone

Alternative models for the variability of winter mortality between years were similar and also listed snow depth and avalanche risk index as the only two explanatory variables, but aggregated over crosswise exchanged periods and/or with different methods. The second-best model (relative increase of AIC by 2%) consisted of fixed effect 1: HS, snow-up season; fixed effect 2: ARI, snow-covered season. Both factors, aggregated over the same periods, were considerably correlated and to a certain degree interchangeable in the models.

3.3. Winter mortality and habitat

Analyzing between-populations effects represents the question whether variations in the mortality rates of different populations can be linked to differences in the climatic suitability of their respective habitats. The analysis is performed on relative winter mortality, i.e. the absolute number of carcasses per population and year divided by the respective population size. This normalization was necessary to compare mortality between populations with different sizes.

¹Southerly aspect means lower mortality rate.

The best-fit model (Table 5) involved a strong association between autumn weather and relative winter mortality: habitats that often see sunny weather conditions during snow-up are inhabited by populations with consistently lower mortality rates. The same applies to habitats with below-average snow depths in winter (see Reference Viviroli, Gurtz and ZappaWiersema, 1989). Neither habitat slope nor altitude related to mortality. On the other hand, south-facing habitats were inhabited by populations with significantly lower mortality rates, although this effect was minor. These findings exactly match the habitat preferences of chamois in winter (e.g. Reference Loison, Jullien, Menaut, Hofgaard, Ball, Danell and CallaghanMichallet and others, 1999). The only alternative model with similar performance (relative increase of AIC by 1.5%) was almost identical but listed accumulated new snow instead of mean snow depth, both aggregated over the snow-covered period.

Table 5.

Linear mixed-effects model for the variability of the relative mortality rate of Alpine chamois between populations. 46% of the variance is accounted for by random and fixed effects together, and 27% by fixed effects alone