Data

This study is based on daily recordings of maximum and minimum temperature and precipitation from 70 meteorological stations in Bhutan covering the period 1990 to 2011 (Figure 1; Supplemental Material, Table S2;  http://dx.doi.org/10.1659/MRD-JOURNAL-D-15-00020.S1). At the stations, daily minimum and maximum air temperature were recorded manually from a thermometer in a standard screen at a 2 m height, while precipitation was measured by a rain gauge with a funnel and a container. The data were supplied by the Meteorology Division, Department of Hydro-Met Services, Ministry of Economic Affairs, Royal Government of Bhutan. Until 1995, about 10 stations were available; beginning in 1996, about 50 stations were available, and another 20 stations became available after 1996. However, the number of years of data from each station used in this study varies, because some years had to be excluded based on our screening of the data, as discussed in the next section.

Other independent variables used in our analysis were elevation, latitude, longitude, slope, and aspect. Station position and elevation were provided by the Department of Hydro-Met Services; station slope and aspect were extracted from a digital elevation model (Aster global digital elevation model [GDEM], version 3, resolution 30 m) using the surface and slope functions in the program ArcGIS from Esri (ArcGIS Desktop, release 10.2.2). These characteristics are shown in Supplemental Material, Table S2 ( http://dx.doi.org/10.1659/MRD-JOURNAL-D-15-00020.S1).

Data screening and validation

Years with fewer than 360 days of observational data were omitted from the analysis. Of the 70 stations, 24 stations had less than 10 years of data. Stations with less than 10 years of data provided 147 station-years of data, compared with 694 station-years for stations with 10 or more years of data. The average annual rainfall of stations with less than 10 years of data was 1597 mm, compared with 1898 mm for stations with 10 or more years of data. This difference was found only during the summer monsoon period, where there are large differences between stations. Average mean temperature was 15.9°C for stations with less than 10 years of data and 17.7°C for stations with data from 10 or more years. This difference was similar in all seasons and reflects that stations with less than 10 years of data were generally located at higher elevations and latitudes than those of stations with longer data series.

There were consistent temperature trends depending on altitude. Sites at high elevations (>3500 m above sea level [masl]) all recorded temperatures of less than −10°C in winter; those at low altitudes (<100 masl) recorded summer temperatures of more than 35°C. Most interior parts of Bhutan receive no rain in winter, whereas the low-elevation southern region receives relatively heavy rainfall in summer. Because of their consistent nature, the data are believed to be of good quality and sufficient for the present study despite the short data series for some stations. To ensure adequate coverage, our statistical analysis included data from all stations, irrespective of data series length, but it was designed to account for the length of the data series (this is discussed in more detail later). Because of these limitations (short and nonuniform data series), we did not attempt to analyze potential trends in climatic conditions across the study area.

Computation of data averages

The daily mean temperature was calculated over the observation period as the average of maximum and minimum temperature, which was subsequently averaged for each station for the 4 seasons: winter (December–February), spring (March–May), summer (June–August), and autumn (September–November). Details are presented in Supplemental Material, Table S1 ( http://dx.doi.org/10.1659/MRD-JOURNAL-D-15-00020.S1).

Test for stationarity of the data series

The autoregressive properties of seasonal temperature and precipitation were tested separately for the 5 stations with more than 20 years of data (Khomachu, Punakha, Samtse National Institute of Education [NIE], Sunkosh, and Wangdue) using the AUTOREG procedure of SAS (SAS Institute, Cary, NC, USA) with the Durbin-Watson statistic. Only stations with more than 20 years of data were included in this analysis, because testing of autoregressive properties requires data series of that length.

The trends for quarterly mean temperature and cumulative precipitation were also tested in this analysis. For the mean temperature, significant first-order positive autoregressive properties were found for spring (1 station), summer (3 stations), and annual (2 stations) temperatures. One of the stations (Wangdue) showed a significant increase (P < 0.05) in mean annual temperature over time for all seasons, whereas another station (Sunkosh) showed a significant decline in summer mean temperature. The seasonal and annual precipitation data were square root transformed to obtain variance homogeneity. For precipitation, significant positive autoregressive properties were found for spring (1 station) and summer (2 stations). One of the stations (Punakha) showed a significant decline over time in precipitation for winter and spring. A significant decline in spring precipitation was found for Khomachu, and a significant increase in summer precipitation was found for Samtse NIE.

Temperature lapse rate

The temperature lapse rate (decrease with altitude) was calculated using linear regression of monthly mean temperature against station elevations:

where T_i and E_i are the mean temperature (in degrees Celsius) and elevation (in meters) of station i, respectively; α is the intercept; β is the slope; and ε_i is the error term. The temperature lapse rate per 100 m was estimated by multiplying β by 100.

Regression models of geographical variation in temperature and precipitation

Because the autoregressive properties for temperature and precipitation were not marked and consistent, we chose to apply a mixed linear model for the analysis of seasonal and annual temperature and precipitation data, using the MIXED procedure of SAS for all 70 stations. (The square root of precipitation was used to obtain variance homogeneity.) Seasonal and annual temperature and precipitation were related to elevation, latitude, latitude squared (to account for nonlinear relationships between seasonal precipitation and latitude), longitude, slope, and aspect. A stepwise procedure with backward elimination was used to exclude variables not contributing significantly to explaining variation in temperature and precipitation. Backward elimination was chosen because the small number of potential variables made it possible to manually test models for significant influence of individual variables and to select models based on the Akaike information criterion (AIC). Models were thus selected that had the lowest AIC and significant (P > 0.05) contribution of independent variables. Station number and year were used as random variables to capture the random structure of the data and to allow for decadal patterns. This was also done to increase the robustness of the parameter estimates and to allow for different periods and durations of the different station series.

The following equation was used for temperature:

where T_ij is the mean temperature for station i in year j (in degrees Celsius); E_i is the elevation (in meters); La_i is the latitude (in degrees); β₀, …, β₂ are coefficients; and A_i and B_j are the random effects of station and year, respectively. The following equation was used for precipitation:

where is the cumulated precipitation for station i in year j; La_i is the latitude (in degrees) at station i; Lo_i is the longitude (in degrees); S_i is the slope (in degrees); β₀, …, β₃ are coefficients; and A_i and B_j are the random effects of station and year, respectively.

The linear equations developed with the linear mixed models (Tables 1, 2) were used to estimate annual and seasonal mean temperature and precipitation across Bhutan. The topographical variables used to compute these climatic variables were derived from the Aster GDEM dataset, and the climate variables were thus mapped at a 30-m resolution.

TABLE 1 

Relationship of seasonal and annual mean temperature (in degrees Celsius) to station elevation and latitude. There were 841 observations for the parameter estimation.

TABLE 2 

Relationship of square root of seasonal and annual mean precipitation (in millimeters) to station latitude, squared latitude, longitude, and slope. There were 841 observations for the parameter estimation.

Logistic regression analysis

Average temperature and precipitation across Bhutan, estimated using these equations, were related to land cover information using the Land Use Planning Project dataset (Figure 2; MOAF 2014). A number of points were targeted randomly in each land cover type in a regression analysis; 1 point per square kilometer was added sequentially, for a total of about 38,000 randomly selected points. For each point, the land cover class (ie vegetation type) and its area were extracted, together with average climate values, from the relevant raster.

FIGURE 2 

Major land cover types in Bhutan. (Map by Ugyen Dorji; data source: MOAF 2014)

The resulting data on land cover and use were analyzed by logistic regression to determine their relation to topographical and climatic variables using JMP software, version 11 (SAS Institute, Cary, NC, USA). The following variables were tested in the logistic regression: elevation, seasonal and annual mean temperature, seasonal and annual rainfall, and topographical position index (TPI, as described later in Equations 4 and 5). The models were evaluated with respect to R² and AIC, which provides a measure to compare performance of different models and to evaluate where the preferred model has the lowest AIC. For each land cover class (MOAF 2014), the mean value and the 10th and 90th percentiles were computed for the following variables: elevation, TPI, annual precipitation, and mean annual temperature. TPI was estimated based on the Aster GDEM using a spatial resolution of approximately 30 m. The computational method applied was that of named deviation from mean elevation (Wilson and Gallant 2000):

where z₀ is the elevation in a raster cell, is the mean altitude within a distance R from this cell, and SD is the local standard deviation of altitudinal values within the same area:

Weiss (2001) and Jenness (2006) applied this method for classification of topographical position and renamed the ratio TPI. De Reu et al (2013) provided a comprehensive overview of the method and its applications in mountainous landscapes. Small TPI values indicate a valley, large values indicate position on a ridge, and intermediate values indicate intermediate slopes or flat areas. In combination with a map of slope gradients, it is possible to classify the landscape into these categories. For this study, the TPI was used directly as a continuous variable.

Temperature

The mean temperature varies considerably both between months or seasons and among geographical areas. This is illustrated in Figure 3 for 4 characteristic stations in Bhutan, where both the mean annual temperature and the temperature range over the year differ greatly among stations. The 4 weather stations recorded mean monthly temperatures ranging from <5 to >15°C in winter and from <15 to >25°C in summer.

FIGURE 3 

Mean monthly temperature (dots) and precipitation (bars) for 4 weather stations in Bhutan representing different temperature regimes. The error bars show the standard deviation.

On average, a 100-m increase in elevation will decrease temperature by 0.42–0.58°C (R² values of 0.71–0.92). The temperature lapse rate has a bimodal pattern (Figure 4). Higher lapse rates are exhibited in February, March, and November, and lower lapse rates occur in January, July, and December.

FIGURE 4 

Monthly estimated lapse rate by regression of average monthly mean, minimum, and maximum temperature (averaged over observation period) against station elevation. The graph to the right shows the R² of the linear regressions.

Using a mixed linear regression model, we tested the effects of various factors (elevation, latitude, longitude, and slope) on mean seasonal and annual temperature for the 70 stations included in our data set. This showed a significant relationship of temperature to elevation and latitude (Table 1). Elevation had a stronger influence on temperature than did latitude. Winter (December–January) mean temperature had a significant relationship with elevation; spring (March–May), summer (June–August), autumn (September–November), and annual mean temperatures had significant relationships with elevation and latitude. The root mean square error (RMSE) ranged from 0.84–1.44°C. Elevation was significantly correlated with latitude (r  =  0.66), which indicates a risk of multicollinearity in the parameter estimation. However, the AIC consistently showed lower values for the models that included both elevation and latitude, indicating that these models had better predictive power than the model that only included elevation. The greatest improvement in AIC was obtained for summer temperatures.

Across Bhutan, the estimated mean temperature based on the temperature models (Table 1) ranged from −20 to 20, 0 to 30, and 0 to 20°C for winter, summer, and annual calculations, respectively (Figure 5). The large spatial variation in temperature and its link to elevation can be clearly seen, with highest temperatures in the lowlands bordering India, as well as in the valleys. The warmest areas were located in the southern foothills, with gradually decreasing temperatures as elevation increases northward across Bhutan. For one of the warmest stations, Phuntsholing (Figure 3), the highest temperatures occurred from May to September, whereas for the high-elevation station at Yotongla, relatively warm conditions occurred only from June to August. Relatively cool conditions at Phuntsholing occurred only in December and January but at Yotongla occurred for a longer period from November to March.

FIGURE 5 

Mean winter (Dec–Jan), summer (Jun–Aug), and annual temperatures (in degrees Celsius) and mean annual precipitation (in millimeters) in Bhutan. (Figure by Peder K. Bøcher; data sources: MOAF 2014 and Meteorology Division, Department of Hydro-Met Services, Ministry of Economic Affairs, Royal Government of Bhutan)

Precipitation

Precipitation varied considerably across Bhutan. Some areas, such as Khomachu, received a mean annual precipitation of less than 800 mm, whereas other areas, such as Phuntsholing, received more than 4000 mm, with 800 mm in the month of July alone (Figure 3). The highest precipitation was recorded in the southern and northern regions of Bhutan (Figure 6). The overall geographical pattern of precipitation distribution was the same for winter, summer, and the full year.

FIGURE 6 

Relationship between seasonal precipitation and latitude for climatological stations in Bhutan.

Our analysis of the relationships between precipitation and geospatial factors attributed variations in precipitation primarily to latitude, with the same quadratic relationship between mean accumulated precipitation and latitude in all seasons (Figure 6). We chose to represent this relationship by a quadratic function (latitude and latitude squared). December–February and March–May mean accumulated precipitation showed significant relationships with station latitude and latitude squared (Table 2); June–August, September–November, and annual mean accumulated precipitation also showed significant relationships with longitude and slope. RMSE ranged from 6.0–31.8 (in millimeters, square root transformed).

Land cover

The effect of geographical and climatic factors on vegetation and land use was explored by logistic regression, which relates a categorical dependent variable (in this case, land cover or use, as shown in Figure 2) to continuous predictor variables (Table 3). Of the individual predictor variables, elevation and annual mean temperature had the largest prediction capacity. With 2 predictor variables, elevation and precipitation or temperature and precipitation provided the best prediction. With 3 variables, the best prediction was obtained using elevation, temperature, and precipitation, and with 4 variables, TPI was included.

TABLE 3 

Land cover and land use predictor variables computed by logistic regression based on elevation, TPI, precipitation, and temperature at annual and seasonal scales.

The areas where vegetation cover and land use are most influenced by human activity (settlements, agriculture, horticulture, plantations, and improved pasture) have mean annual temperatures above 9°C and below 24°C (Table 4). However, annual precipitation varies greatly within each of these land use classes. Improved pasture is found in areas with temperatures lower than those in the other land use classes—at higher elevations but, according to the TPI, still in valleys. Settlements are also confined to valleys, whereas agriculture and plantations are also found on the sloped areas. Agriculture and horticulture occur at elevations lower than those of other land uses.

TABLE 4 

Elevation, TPI, and annual precipitation and temperature values for the different land cover and use classes (%, percentile). (Table extended on next page.)

TABLE 4 

Extended.

Other natural land cover classes showed greater variation than those influenced by humans. For example, forested areas varied greatly in both temperature and precipitation. Nevertheless, the different forest types separated well in terms of temperature and precipitation (Table 4). Broadleaf forest was found in regions with temperatures equivalent to those of agricultural areas and with more rainfall than occurs in conifer areas. The different conifers were largely found to be separated in terms of temperature, with chir pine in warmer areas, blue pine in cooler areas, and fir in cold areas. Mixed conifer forest covered a wider range of cold to cool temperatures, and the lowest temperatures were found for scrub forest, which also extended to the highest elevations. Scrub forests were also mostly found in areas with steeper slopes (greater TPI values) than in other forest types.

Meadows and marshland areas were found in areas with mostly sloping terrain and in cool or cold temperatures, with marshland having generally colder conditions than meadows. The marshland extends to elevations also covered by snow and glaciers.

The remaining land covers (rock outcrops, landslips, bodies of water, and glaciers or snow) were generally found in areas with rainfall higher than that in most other land cover classes (Table 4). Bodies of water (lakes and rivers) were found in valley bottoms not associated with any specific temperature range. Glaciers were found in areas with average temperatures of about 0°C or lower, rock outcrops were in regions with steep slopes, and landslips were on lower slopes.

Temperature

The spatial and temporal temperature pattern in mountainous regions is complex because of the influence of geographical position, elevation, aspect (slope direction), and slope angle (Dobrowski et al 2009). The highly mountainous and broken terrain in the Himalayas causes an extreme variation in climate, making it difficult to adequately measure and analyze the various climate parameters (Lambert and Chitrakar 1989). This is true for temperature variability in Bhutan (Figures 3, 5).

The strong influence of elevation on temperature, as seen in our study for Bhutan (Table 1), has also been reported for Nepal (Kattel et al 2013). The decrease in temperature with altitude is not only a fundamental physical property (Blandford et al 2008; Dobrowski et al 2009) but also a major factor in an ecosystem’s suitability for different species. This is seen in particular in the separation of different forest types based on temperature (Table 4).

The existence of a positive, albeit somewhat less significant, relationship of temperature with latitude may be because of the decrease of precipitation from south to north (Figure 6), where the north-central regions could have less cloud cover, in terms of extent and duration, and thus more sunshine than the south because of heavy precipitation in the south, especially during the summer monsoon season. The effect of latitude on mean temperature is therefore also stronger during the summer season.

The lapse rate of 0.42 to 0.58°C for each 100 m of change in elevation in Bhutan is slightly smaller than the widely used lapse rate of 6.5°C km⁻¹ (Blandford et al 2008; Dobrowski et al 2009). The high correlation between temperature and elevation (R²  =  0.71–0.92) (Figure 4), with a standard error of 1.30 to 1.62°C for predicting temperature from elevation alone, suggests that temperature in our study area is mostly controlled by elevation. In other regions of the world, this may not always be the case. Rolland (2003) reported an annual lapse rate of 0.54 to 0.58°C per 100 m, similar to ours, in alpine regions. However, Dobrowski et al (2009) reported substantially different lapse rates from the United States; Blandford et al (2008) reported monthly lapse rates for Idaho that are markedly lower than ours for minimum and average temperatures and higher for maximum temperature.

In our study, higher lapse rates were consistently observed for minimum temperatures (0.44–0.58°C) and lower lapse rates were seen for maximum temperatures (0.42–0.54°C) (Figure 4), as previously described by Diaz and Bradley (1997). The R² values for the mean and minimum temperatures are also consistently higher than those for the maximum temperature, indicating a stronger correlation between elevation and minimum temperature, which may be related to differential effects of radiation and cloudiness on minimum and maximum temperatures. This may also be related to the effects of air humidity on lapse rate, which probably has greater significance for the maximum temperature, where the relative humidity may vary more during daytime than at night.

The temperature lapse rate has a bimodal annual pattern, as previously reported by Kattel et al (2013). Higher lapse rates were exhibited before and after the summer monsoon season, and lower rates were seen during the summer monsoon, as well as during winter. In contrast, Blandford et al (2008) reported that the lapse rate for maximum temperature is greater in summer than in winter, while minimum and average temperature lapse rates are higher in spring and lower in midsummer. Rolland (2003) also reported higher temperature lapse rates in summer. This indicates that a relatively complex pattern of temperature lapse exists in different regions of the world, with temperature lapse rates and patterns governed by other factors, such as local conditions and climate; the quality and extent of the data and the analysis method used may also play a role (Rolland 2003).

Precipitation

The north–south gradient in precipitation across Bhutan is pronounced (Figure 5), as shown clearly through the strong quadratic relationship with latitude (Figure 6). Latitude has the strongest effect on precipitation; other variables, such as longitude and slope, have a weaker effect (Table 2). However, the lower density of weather stations in northern Bhutan increases uncertainty about the mapped precipitation in this region. It is likely that the quadratic function used to describe the relationship between latitude and precipitation led to an overestimation at higher latitudes in northern Bhutan. The high rainfall rates in southern Bhutan are consistent with measurements in south Asian monsoonal regions such as India, where precipitation of 770 to 1020 mm in a single day has been recorded at many weather stations (De et al 2005).

The strong correlation with latitude may be partly explained by the monsoon. During summer, when rainfall is strongest, the monsoon brings moist air from the Indian Ocean northward (Webster et al 1998), and the humidity of the air declines with greater distance from the ocean. The mountains are believed to be responsible for the northward extension of the monsoon climate on the Asian continent (Hahn and Manabe 1975), and the mesoscale (up to 1 km in height) mountain ranges have been identified as the main anchors for monsoon rainfall in South Asia (Xie et al 2006). Thus, the nonlinear relationship of precipitation with latitude in Bhutan may be linked to the mesoscale mountain ranges anchoring monsoon rainfall and resulting in rapid reduction of precipitation once the monsoon advances over mountains northward in summer. This in principle holds true for the Himalayan winter monsoon that has maximum precipitation in northern hilly regions, with decreasing values toward the south. Precipitation is within the range of 100–150 mm in December–March, and this is important for the winter crops and for replenishing the Himalayan snow, which provides water in other seasons (Yadav et al 2012).

Precipitation shows a negative relationship with longitude: monsoon rain increases westward. This may be linked to the summer monsoon wind blowing from the southwest to the northeast (An 2000). Bhutan seems to receive less precipitation in the east than in the west, although no general trend has been identified for the Himalayan range (Anders et al 2006). Moreover, the orographical effect in the Himalayas is complex and does not present a uniform obstacle to air flow (Bhatt and Nakamura 2005). Therefore, Bhutan could have its own precipitation trend because of the complex terrain and high variability in climatic conditions. Such detailed trends may only be detectable when studying a smaller geographical area, as in the present study, and may not be evident when compared to the whole Himalayan range.

The negative correlation of precipitation with slope angle may be explained by a sheltering effect of stations located near steep slopes. However, the effect of slope is smaller than that of latitude.

We did not find significant relationships between elevation and precipitation. The anticipated orographical effects on precipitation may not be related to the elevation of the specific site, because there may be a spatial displacement of the effect of elevation on precipitation. This displacement could be particularly large in Bhutan, where there is a great variation in elevation over short distances. However, there was a significant tendency toward reduced precipitation on steeper slopes in summer and autumn. This may support the assumption of a spatial displacement between the effect of elevation and the location of rainfall.

Effects of climate and topography on land cover

Land cover is affected by both topographical and climatic conditions, as shown by the predictor variable of the logistic regression in Table 3. Human land use is mostly limited to valleys in the warmer southern and central parts of the country. The average annual temperatures of 17–18°C for the agricultural and horticultural areas indicate that these are mostly used for warm-season crops such as rice and maize (Moore and Lobell 2014; Sánchez et al 2014), whereas cool-season agricultural areas are largely limited to pastures. This may be linked with the developments of agriculture in Bhutan, which likely followed similar patterns to those of other regions in Asia, such as southwest China (Guedes and Butler 2014), where temperatures largely determined the spread of rice cultivation.

Likely effects of climate change

Because temperature correlates strongly to elevation, with an average temperature lapse rate of roughly 0.5°C per 100 m of change in elevation, land cover types may respond differently to climatic warming, unless the land cover is determined by other geographical characteristics. Many of the land cover classes seem to be highly influenced by temperature, and it is therefore likely that they will shift with climatic warming.

Changes in land cover linked to changes in precipitation are more difficult to estimate, because there is only a weak and uncertain correlation between precipitation and vegetation type. Nevertheless, ecosystems such as conifer forest and improved pasture are located in regions receiving relatively low amounts of precipitation and have relatively narrow ranges between the 10th and the 90th percentiles, indicating higher sensitivity to precipitation. Conifers, glaciers or snow, marshland, improved pasture, horticulture, and agriculture have narrow temperature ranges between the 10th and the 90th percentiles, indicating that they are sensitive to temperature change (Table 4).

The present relatively low elevation range of agriculture and horticulture indicates that there may be limitations to the expansion of these land uses besides climatic ones. This could be linked to land accessibility and infrastructure.

The improved pasture system has a temperature range of less than 6°C (Table 4) and could be out of its optimum temperature range within less than 100 years, if a temperature increase of 0.06°C/y (Shrestha et al 2012) continues into the future. In contrast, water bodies have a temperature range of more than 20°C and will be out of temperature range in about 400 years. Improved pasture has a less than 500 mm annual precipitation range (Table 4) and could also be out of its precipitation range in less than 100 years, if precipitation continues to increase by 7 mm y (Shrestha et al 2012). Land cover may change dramatically in location, size, and characteristics with changes in precipitation and temperature, because land cover is highly influenced by topographical and climatic conditions (Table 4).