Robin Pouteau

Serge Rambal

Jean-Pierre Ratte

Fabien Gogé

Thierry Winkel
email: thierry.winkel@ird.fr

Downscaling MODIS-derived maps using GIS and boosted regression trees: the case of frost occurrence over the arid Andean highlands of Bolivia

Keywords: altiplano, Andes, Bolivia, boosted regression trees, DEM, downscaling, frost risk mapping, MODIS, physiography, temperature lapse rate, topoclimate model, satellite land surface temperature, seasonal variation, spatial variation

Frost risk assessment is of critical importance in tropical highlands like the Andes where human activities thrives at altitudes up to 4200 m, and night frost may occur all the year round. In these semiarid and cold regions with sparse meteorological networks, remote sensing and topographic modeling are of potential interest for understanding how physiography influences the local climate regime. After integrating night land surface temperature from the MODIS satellite, and physiographic descriptors derived from a digital elevation model, we explored how regional and landscape-scale features influence frost occurrence in the southern altiplano of Bolivia. Based on the high correlation between night land surface temperature and minimum air temperature, frost occurrence in early-, middle-and late-summer periods were calculated from satellite observations and mapped at a 1-km resolution over a 45000 km² area. Physiographic modeling of frost occurrence was then conducted comparing multiple regression (MR) and boosted regression trees (BRT). Physiographic predictors were latitude, elevation, distance from salt lakes, slope steepness, potential insolation, and topographic convergence. Insolation influence on night frost was tested assuming that ground surface warming in the daytime reduces frost occurrence in the next night. Depending on the time period and the calibration domain, BRT models explained 74% to 90% of frost occurrence variation, outperforming the MR method. Inverted BRT models allowed the downscaling of frost occurrence maps at 100-m resolution, illustrating local processes like cold air drainage. Minimum temperature lapse rates showed seasonal variation and mean values higher than those reported for temperate mountains. When applied at regional and subregional scales successively, BRT models revealed prominent effects of elevation, latitude and distance to salt lakes at large scales, whereas slope, topographic convergence and insolation gained influence at local scales. Our results highlight the role of daytime insolation on night frost occurrence at local scale, particularly in the early-and mid-summer periods when solar astronomic forcing is maximum. Seasonal variations and interactions in physiographic effects are also shown. Nested effects of physiographic factors across scales are discussed, as well as potential applications of physiographic modeling to downscale ecological processes in complex terrains.

Introduction

Low air temperature is one of the most important factors controlling vegetation zonation and key processes such as evapotranspiration, carbon fixation and decomposition, plant productivity and mortality in natural and cultivated mountain ecosystems [START_REF] Chen | Microclimate in forest ecosystem and landscape ecology[END_REF]Nagy et al., 2003).

Depending on vegetation structure, landscape position or soil properties, frost can damage plant tissues thus affecting forest, pasture and crop productivity [START_REF] Blennow | Models of low temperature and high irradiance and their application to explaining the risk of seedling mortality[END_REF]. These damages have consequences for human populations, particularly in the tropics where highlands often remain densely populated [START_REF] Grötzbach | Mountain peoples and cultures[END_REF]. In the Andes of Argentina, Bolivia, Chile, Ecuador and Peru agriculture thrives at altitudes up to 4200 m [START_REF] Del Castillo | Quinoa in Bollivia: an ancestral crop changed to a cash crop with "organic fair-trade[END_REF] and treeline reaches its world's highest elevation up to 5100 m [START_REF] Hoch | Growth, demography and carbon relations of Polylepis trees at the world's highest treeline[END_REF] in spite of night frost occurring on more than 300 days practically spread all over the year [START_REF] Garcia | Agroclimatic constraints for rainfed agriculture in the Bolivian Altiplano[END_REF][START_REF] Gonzalez | Altitudinal and seasonal variation of protective and photosynthetic pigments in leaves of the world's highest elevation trees Polylepis tarapacana (Rosaceae)[END_REF][START_REF] Rada | Low temperature resistance in saplings and ramets of Polylepis sericea in the Venezuelan Andes[END_REF][START_REF] Troll | The Cordilleras of the tropical Americas: aspects of climatic, phytogeographical and agrarian ecology[END_REF]). In the southern Andes, sparsely vegetated areas juxtaposing extended flat plains around salt lakes and steep slopes on the cordilleras and volcanos, display semi-arid and desert landscapes largely dominated by terrain structure. Subjected to the night/day and sunlit/shaded slope contrasts characteristic of the mountain climate, this environment is well suited for examining the influence of regional and landscape-scale physiography on the local climate regime, and particularly frost occurrence.

Several studies on topoclimate in highlands showed that elevation and slope are the main explanatory variables in modeling local climate spatial variability [START_REF] Chuanyan | Methods for modelling of temporal and spatial distribution of air temperature at landscape scale in the southern Qilian mountains, China[END_REF]. By means of digital elevation models and astronomical equations, the potential insolation (incoming solar radiation) has been included as an additional independant variable in some of these models, substantially improving their capacity to predict free-air as well as soil-surface temperature distributions [START_REF] Benavides | Geostatistical modelling of air temperature in a mountainous region of Northern Spain[END_REF][START_REF] Blennow | Models of low temperature and high irradiance and their application to explaining the risk of seedling mortality[END_REF][START_REF] Fridley | Downscaling climate over complex terrain: high finescale (<1000 m) spatial variation of near-ground temperatures in a montane forested landscape (Great Smoky Mountains)[END_REF][START_REF] Fu | A geometric solar radiation model with applications in agriculture and forestry[END_REF]. Among the physiographic variables, elevation and slope steepness are known to influence cold air drainage at night and, hence, the distribution of frost risks at landscape scale [START_REF] Lundquist | Automated algorithm for mapping regions of cold-air pooling in complex terrain[END_REF]Pypker et al., 2007a). In the daytime, slope aspect and physiographic shading effects control the effective radiation load per unit of soil areas, resulting in very contrasted values of daily maximum soil temperature [START_REF] Fu | A geometric solar radiation model with applications in agriculture and forestry[END_REF].

Minimum night temperature might be sensitive to insolation during the previous day, since soil surface warming during that day could dampen soil radiative cooling in the next night. Though challenged by studies on minimum air temperature variations in moderately high mountains under temperate climate [START_REF] Blennow | Modelling minimum air temperature in partially and clear felled forests[END_REF][START_REF] Dobrowski | How much influence does landscape-scale physiography have on air temperature in a mountain environment[END_REF], this hypothesis should be tested in the central part of the Andes, where low latitude, high elevation (typically ranging between 3600 and 4200 m) and sparse vegetation result in much greater radiation load and thermal contrasts across shaded and sunlit areas.

Besides, using a downscaling approach, [START_REF] Fridley | Downscaling climate over complex terrain: high finescale (<1000 m) spatial variation of near-ground temperatures in a montane forested landscape (Great Smoky Mountains)[END_REF] noticed that the lack of relationship between daytime radiation and nighttime temperature is true at local scale lower than 1000 m but not at regional scale (Great Smoky Mountains, USA), where variations in radiation balance across locations do influence nighttime temperature distribution particularly in cooler situations. Considering the Andes, recent work by [START_REF] Bader | A topography-based model of forest cover at the alpine tree line in the tropical Andes[END_REF] and Bader et al. (2008) used topographic modeling and remote sensing data to examine the response of vegetation distribution to climate warming, but we should go back to [START_REF] Santibañez | Topoclimatic modeling for minimum temperature prediction at a regional scale in the Central Valley of Chile[END_REF] and [START_REF] François | Frost risk mapping derived from satellite and surface data over the Bolivian Altiplano[END_REF] to find studies on the links between frost climatology and physiography over this region. These early works were not continued, and the case of the Andean highlands remained poorly documented in spite of the potential interest of that region, densely populated and representative of the tropical mountains vulnerable to global warming [START_REF] Vuille | Climate Change and tropical Andean glaciers: past, present and future[END_REF].

Analyzing topographic effects on free-air or land surface temperature also led to reevaluate the simplifying assumption of a generic environmental lapse rate (the decrease in free-air temperature as elevation rises, typically assumed to be -0.6 °C per 100 m), commonly applied in hydrological and ecological studies to extrapolate air temperature in mountain areas. In fact, several studies show temperature lapse rate variations due to seasonality, height above the ground, or ground surface characteristics [START_REF] Blandford | Seasonal and synoptic variations in near-surface air temperature lapse rates in a mountainous basin[END_REF][START_REF] Dobrowski | How much influence does landscape-scale physiography have on air temperature in a mountain environment[END_REF][START_REF] Fridley | Downscaling climate over complex terrain: high finescale (<1000 m) spatial variation of near-ground temperatures in a montane forested landscape (Great Smoky Mountains)[END_REF][START_REF] Lookingbill | Spatial estimation of air temperature differences for landscapescale studies in montane environments[END_REF][START_REF] Marshall | Near-surface-temperature lapse rates on the Prince of Wales Icefield, Ellesmere Island, Canada: implications for regional downscaling of temperature[END_REF], though no detailed reports were published for the Central Andes. Most of the above mentionned studies used multiple linear regression for modeling the influence of physiography on free-air or soil-surface temperatures. The present study resorts to an advanced form of regression, the boosted regression trees (BRT). BRT use the boosting technique to combine large numbers of relatively simple tree models to optimize predictive performance. BRT have been used successfully in human biology [START_REF] Friedman | Multiple additive regression trees with application in epidemiology[END_REF], land cover mapping [START_REF] Lawrence | Classification of remotely sensed imagery using stochastic gradient boosting as a refinement of classification tree analysis[END_REF], biogeography [START_REF] Parisien | Environmental controls on the distribution of wildfire at multiple spatial scales[END_REF], species distribution [START_REF] Elith | A working guide to boosted regression trees[END_REF], and soil science [START_REF] Martin | Optimizing pedotransfer functions for estimating soil bulk density using boosted regression trees[END_REF]. They offer substantial advantages over classical regression models since they handle both qualitative and quantitative variables, can accomodate missing data and correlated predictive variables, are relatively insensitive to outliers and to the inclusion of irrelevant predictor variables, and are able to model complex interactions between predictors [START_REF] Elith | A working guide to boosted regression trees[END_REF][START_REF] Martin | Optimizing pedotransfer functions for estimating soil bulk density using boosted regression trees[END_REF]. Though direct graphic representation of the complete tree model is impossible with BRT, the model interpretation is made easy by identifying the variables most relevant for prediction, and then visualizing the partial effect of each predictor variable after accounting for the average effect of the other variables [START_REF] Friedman | Multiple additive regression trees with application in epidemiology[END_REF].

The aims of the present work were: i) to explore how regional and landscape-scale physiography influence frost occurrence in Andean highlands through integration of field and remote sensing data, digital terrain analysis, and GIS, ii) to downscale regional frost occurrence maps at a level relevant for farming and land management decisions using BRT models. This study was focused on the austral summer period, from November to April, when frost holds the greatest potential impact for local farming activities.

Material and methods

Study area and regional climate

The study area was located at the southwest of the Bolivian highlands, near the borders of Argentina and Chile,between 19°15' and 22°00' South and between 66°26' and 68°15' West. This region, boarded by the western Andes cordillera, is characterized by the presence in its centre of a ca.100 x 100 km dry salt expanse, the Salar of Uyuni, while another salt lake, the Salar of Coipasa, lies at the north of the study area. The landscapes show a mosaic of three types of land units: more or less extended flat shores surrounding the salt lakes (elevation ca. 3650 m) and an alternation of valleys and volcanic relieves (culminating at 6051 m) in the hinterland. The native vegetation of this tropical Andean ecosystem, also known as puna, consists of a mountain steppe of herbaceous and shrub species (e.g. Baccharis incarum, Parastrephia lepidophylla, and Stipa spp.) [START_REF] Navarro | Mapa de vegetación de Bolivie a escala 1:250.000[END_REF] traditionally used as pastures but progressively encroached by the recent and rapid expansion of quinoa crop (Chenopodium quinoa Willd.) [START_REF] Vassas | El auge de la quinua: cambios y perspectivas desde una visión social[END_REF].

Due to its low latitude and high elevation, the study area is characterized by a cold and arid tropical climate. Average precipitations vary between 100 and 350 mm year -1 from the South to the North of the region [START_REF] Geerts | Agro-climatic suitability mapping for crop production in the Bolivian Altiplano: a case study for quinoa[END_REF], presenting an unimodal distribution with a dry season from April to October.

The annual average temperature (close to 9 °C) hides daily thermal amplitudes higher than seasonal amplitudes, of up to 25 °C [START_REF] Frère | Estudio agroclimatológico de la zona andina[END_REF]. These particular thermal conditions lead to high frost risks throughout the year. Advections of air masses from the South Pole represent only 20% of the observed frosty nights [START_REF] Frère | Estudio agroclimatológico de la zona andina[END_REF] and are four times less frequent in the summer than during the winter, when the intertropical convergence zone goes northward [START_REF] Ronchail | Advections polaires en Bolivie : mise en évidence et caractérisation des effets climatiques[END_REF]. Therefore, the main climatic threat lies in radiative frost, occurring during clear and calm nights. As reported by local peasants, frost occurence shows a strong topographical and orographical dependence, as well as a marked seasonality. This seasonality lead us to split the active vegetation period into three time periods characterizing the mean regional climate dynamics in the summer rainy season: November-December when precipitation and minimum temperature rise progressively, January-February when precipitation and temperature are at their maximum, and March-April when both begin to decrease.

Data

Meteorological ground data

In the study area, daily air temperature records were available in three meteorological stations: one November 2006 to 31 December 2007 in Jirira. This dataset was temporally and spatially insufficient to interpolate frost risks at a regional scale, but it allowed to establish the relationship between air temperature and remotely sensed land surface temperature.

Remotely sensed data

The two sensors, Terra and Aqua, of the satellite system MODIS give daily images of the Earth radiative land surface temperature. Images from the fifth version of the MYD11A1 MODIS product were concatenated and projected in the UTM-19S (Universal Transverse Mercator 19 South) coordinate system using the MODIS projection tool. In this way, daily 1-km resolution images of the radiative land surface temperature (Ts) over the study area were obtained. Ts images recorded by the Aqua sensor around 2 a.m. were used as they were closer to the Tn data recorded at ground level and closer to the true physiological conditions experienced by the vegetation [START_REF] François | Frost risk mapping derived from satellite and surface data over the Bolivian Altiplano[END_REF]. Time series of nominal 1-km spatial resolution MODIS data were downloaded from NASA's EOS data gateway (https://wist.echo.nasa.gov/) from 20 July 2001 to 25 April 2006 and from 01 January 2007 to 31 December 2007. Due to the particular surface properties of the salt lakes of Coipasa and Uyuni in terms of surface moisture and radiative emissivity, parameter estimations were considered dubious there and Ts data for the salt lakes were discarded from the analysis. This database was managed and analyzed using the ENVI 4.2. software (ITT Visual Information Solutions, www.ittvis.com). The statistical correspondence between Tn data recorded in the meteorological stations of Salinas, Irpani and Jirira and Ts data of the pixels including these three localities was examined by linear regression and Pearson correlation.

Apart from Ts measurements during clear nights, MODIS images also bring information about the possible presence of clouds between the Earth surface and the satellite at the time of the record. The information of these "flagged" images is valuable for our purpose since radiative frost would not occur during cloudy nights. The frequency of cloudy pixels in the daily MODIS images was thus calculated and used in the frost occurrence calculation (see below).

Digital elevation model and physiographic predictors

The SRTM digital elevation model [START_REF] Farr | The shuttle radar topography mission[END_REF] with a 90 m horizontal resolution and a vertical accuracy better than 9 m was used after resampling to 100 m to make easier the correspondence between the digital elevation model and the MODIS images at 1-km scale. In a GIS environment using Idrisi Kilimanjaro, Envi 4.2. and ArcMap 9.2. softwares, eight physiographic variables were calculated at a 100-m resolution for each location to examine their potential role in the spatial determinism of frost and to downscale frost maps to levels closer to those of frost impacts on anthropic activities (Table 1). The compound topographic index (CTI) was used as an index of cold air drainage [START_REF] Gessler | Modeling soillandscape and ecosystem properties using terrain attributes[END_REF], with low CTI values representing convex positions like moutain crests and with high CTI values representing concave positions like coves or hillslope bases. Three insolation variables (DPI, MPI and API) were calculated by the ArcMap 9.2. solar analysis tool. They express the amount of radiative energy received across all wavelengths over the course of a typical seasonal day (DPI), or from sunrise to 12:00 (MPI), or from 12:00 to sunset (API) of such a day. These insolation variables account for site latitude and elevation, slope steepness and aspect, daily and seasonal sun angle, and shadows cast by surrounding heights. API was calculated to examine the specific influence of insolation in the afternoon just before the considered night, with the hypothesis that high soil surface insolation and warming would affect the soil energy balance, and thus reduce the risk of radiative frost in the following night (without regards to other potential factors such as soil albedo, soil water content, air humidity, etc. see [START_REF] Garcia | Dynamics of reference evapotranspiration in the Bolivian highlands (Altiplano)[END_REF]). Similarly, MPI was calculated as a surrogate to the early morning insolation, with the hypothesis that areas in the shade of surrounding heights in the early morning would experience cooler conditions for a longer time, thus being more vulnerable to frost than sunlit areas. In the calibration procedure, these 100-m resolution variables were upscaled at 1-km resolution by averaging 10 x 10 pixel clusters in the DEM 100-m images, thus fitting the 1-km resolution of the remotely sensed frost occurrence maps.

Physiographic modeling of frost occurrence over regional and subregional domains 2.3.1. MODIS-derived frost occurrence

Frost is detected by remote sensing when surface temperature appears negative on cloudfree images. Based on the standard meteorological threshold of 0 °C, frost occurrence (R) for a specific time period was therefore defined as follows: R = Prob (Ts < 0 °C) * F (1)

where: R = frost occurrence at the 0 °C threshold (relative probability ranging from 0 to 1), Prob (Ts < 0 °C) = probability of the surface radiative temperature being lower than 0°C, F = frequency of cloudless days in the considered period. Note that "frost occurrence" is used here instead of frost risk to differentiate our estimates based on 6-year daily Ts values from climatological estimates based on longer data series.

In order to calculate the probability Prob (Ts < 0 °C) , the distribution of the random variable Ts during successive time periods (namely: November-December, January-February, March-April) was studied, checking its normality through the Kolmogorov-Smirnov test. For each 1-km pixel and each time period, Ts mean and standard deviation, cloudless day frequency (F) and, finally, frost occurrence (R) were calculated from the available nighttime remotely sensed data series (n = 366 , 355, and 361 in the ND, JF and MA periods respectively). Maps of observed (remotely sensed) frost occurrence at 1-km resolution were then generated by applying Eq. ( 1) for each time period.

Frost occurrence models over regional and subregional domains

A subsample of 1-km pixels (n = 7500) was randomly selected for the calibration of the physiography-frost occurrence relationships over the entire study area (hereafter called "regional models"). Regional BRT were built for each seasonal period (November-December, January-February, March-April) using the gbm package version 1.6-3 developped under R software (R Development Core

Team 2006). A bag fraction of 0.5 was used which means that, at each step of the boosting procedure, 50% of the data in the training set were drawn at random without replacement. The loss function (LF),

defining the lack-of-fit, used a squared-error criterion. The learning rate or shrinkage parameter (LR), the tree size or tree complexity (TS), the number of trees (NT) and the minimal number of observations per terminal node (MO) were the main parameters for these fittings, and were set through a tuning procedure [START_REF] Martin | Optimizing pedotransfer functions for estimating soil bulk density using boosted regression trees[END_REF]. LR, determining the contribution of each tree to the model, was thus taken equal to 0.05. NT, the maximal number of trees for optimal prediction was set to 2000. For optimal prediction, TS, the maximal number of nodes in the individual trees, was set to a value of 9, and MO was set to 10 observations per terminal node. For sake of comparison, multiple linear regression models were calculated at the regional scale using the same predictors and the same calibration datasets as for the regional BRT. These " regional MR" were built using the Statistica package (StatSoft France 2005).

The regional BRT and MR models were validated comparing observed (remotely sensed) and predicted frost occurrence over the entire study area in the three time periods. This was made excluding the pixels used for calibration, which resulted in a 49353 pixels validation set. The predictive capacity of the models was analysed examining the observed versus predicted values plots, the bias (B), the root mean square error of prediction (RMSE), and the coefficient of determination of the regression between estimated and observed values (R 2 ). Once validated, the BRT were interpreted, looking first at the relative contribution of the physiographic variables to the predictive models, and then considering the partial dependence of the predictions on each variable after accounting for the average effect of the other variables.

In order to test the scale-dependence of the predictors, a similar BRT procedure was applied over a smaller spatial domain defined by a selected range of regionally varying factors, namely: latitude between 19°5 and 20° South and elevation lower than 4200 meters (total area = 7775 km²). This spatial domain corresponds to the Intersalar, the area of major agricultural activities in the region, where local populations cultivate quinoa and rear camelids up to an altitude of ca. 4200 meters. A new set of 7500 training pixels was randomly selected from this smaller domain to calibrate these "subregional BRT", using the same values of fitting parameters and a similar validation procedure as in the previous analysis. Excluding the training pixels, this validation was conducted on the remaining 275 pixels of this smaller domain. The relative contribution of the physiographic predictors to the subregional models was also examined. The interactions between predictive variables were considered by joint plots of their partial dependence in the subregional models.

Downscaling frost occurrence prediction at 100-m resolution

Once validated, the regional BRT were applied on each pixel of the DEM 100-m image in order to downscale frost occurrence from 1-km to 100-m resolution in the three considered time periods. An indirect validation of these 100-m frost occurrence maps was then conducted by aggregating 10 x 10 pixel clusters of 100-m frost predictions and comparing the resulting 1-km predictions to the observed (remotely sensed) frost occurrence at 1-km resolution. A qualitative validation was also conducted by examining the capacity of these 100-m maps to display well-known local patterns of frost and cold air distribution over complex terrains.

Estimation of the land surface temperature lapse rate

The regression of Ts recorded over sloping areas versus elevation of the corresponding pixels allowed to calculate average values of the land surface temperature lapse rate at night for successive dates in each considered time periods. Sloping areas were defined as terrains with slope steepness greater than 3° and elevation lower than 5000 meters. This elevation limit was chosen to discard highaltitude sites possibly covered with snow or ice which superficial thermal properties modify lapse rate estimations [START_REF] Marshall | Near-surface-temperature lapse rates on the Prince of Wales Icefield, Ellesmere Island, Canada: implications for regional downscaling of temperature[END_REF]. The resulting sampling area represented 22529 km², covering an elevation range of 1341 m (from 3659 to 5000 m).

Results

Climate information and remotely sensed frost occurrence evaluation

The frequency analysis of daily minimum air temperature (Tn) recorded at Salinas over 10 discontinuous years (1989, and the 1998-2006 period) was made using the standard climatological threshold of 0 °C (Fig. 1). the Ts data series was statistically significant in all the cases (P < 0.05). This allowed to apply a normal probability density function in equation ( 1) in order to generate the 1-km resolution maps in Fig. 2 showing the regional patterns of frost occurrence variations in three successive time periods as derived from satellite observations.

Physiographic modeling of frost occurrence

Model validation of regional and subregional BRT models

The results of the statistical comparison between observed and predicted frost occurrence at 1-km resolution in three time periods are presented in Fig. 3 and Table 2. The regional BRT clearly outperformed the MR models, the latter being affected by strong non linearities in the high frost occurrence range, showing its poor predictive capacities in the early and late summer periods. On the other hand, the regional BRT negligibly overestimated the satellite observations with practically no bias whatever the time period. The RMSE and R² values showed that BRT predictions were fairly good in January-February (RMSE = 0.057, R² = 0.90), and only slightly more dispersed by the beginning or the end of the summer season (RMSE between 0.07 and 0.08, R² between 0.78 and 0.83). With an error generally less than 8% on predicted frost occurrence values, the regional BRT thus appear suitable for predicting frost occurrence from physiographic variables alone. Table 2 shows similar performances of the regional and subregional BRT, with only higher bias for the subregional model.

Hierarchy of physiographic variables in regional and subregional BRT

Regional BRT showed large effects of elevation, distance from the salt lakes, latitude and CTI on frost occurrence, while slope and insolation variables had only marginal influence (Fig. 4a). Comparing the three time periods, the relative contributions of the predictive variables showed some variations, with the distance from the salt lakes dominating in the initial period (November-December), while elevation became more important from January to April, and particularly in the mid-summer period. In contrast, latitude, CTI and insolation variables kept a fairly constant effect with similar contributions at the beginning and at the end of the season. Calibrating BRT models over a limited spatial domain reveals slightly different patterns in the contributions of the predictors (Fig. 4b): distance from the salt lakes gained importance on elevation in the three time periods, and daily potential insolation showed noticeable contribution until mid-summer (though the influence of its morning and afternoon components remained marginal). The weights of CTI and latitude were intermediate whatever the time period, while the contribution of slope became important by the end of the season. The plots of partial dependence for frost occurrence in the regional BRT (Fig. 5 see Table 1 for variables abbreviations).

More details emerge from the partial dependence plots of interactions in subregional BRT (Fig. 6 illustrating the January-February period, with similar results in the other two periods). For instance, up to a distance of 10 km from the salt lakes (LDS ≈ 4.0) the effect of elevation on frost occurrence is low and nearly constant below 3900 m, and then rises gradually above that level. But farther than 10 km away from the salt lake borders, frost occurrence first decreases as elevation rises up to 3900 m and then increases sharply above. This suggests that thermal inversions at night are more frequent at distance from the salt lakes. Considering the interaction of elevation with CTI, while frost occurrence at low elevation increases gradually up to CTI values of 9 and more sharply thereafter (concave situations), at high elevation the effect of CTI is already important at values below 8 and then rose only marginally. The effect of landscape concavity thus appears prominent at low elevation, where cold air can accumulate in local depressions, while it becomes negligeable at high elevation where crests and peaks dominate. Fig. 6. Joint partial dependence plots of some interactions between topographic variables in the subregional BRT of the January-February period (see Table 1 for variables abbreviations).

Table 2. Validation statistics of frost occurrence multiple regression (MR) and boosted regression trees (BRT) models calibrated over the entire study area (regional models) or the Intersalar area (subregional models). B: bias; RMSE: root mean square error of prediction; R2: determination coefficient of the regression line between observed and predicted values. ND = November-December; JF = January-February; MA = March-April.

Fine resolution frost occurrence maps

Regional BRT were used in prediction to downscale frost occurrence maps from the 1-km to the 100-m scale. The statistical validation conducted on reaggregated 1-km pixel clusters shows good fit between predicted and observed (remotely sensed) frost occurrence values (Table 3), with RMSE of predicted values of the same order than in the BRT regional model directly applied at the 1-km resolution (Table 2). The 100-m scale maps display topoclimatic variations resulting in a detailed zonation of frost occurrence. As an example, Fig. 7c-d shows that flat areas surrounding Mount Tunupa are more prone to frost occurrence than the slopes of the volcano up to an elevation of approximately 4000 m, while sites located at higher altitudes are naturally colder. All over this area, east-facing slopes appear less exposed to frost than west-facing slopes. In some particular places at the west, cold air stagnation is also identifiable in the lower parts of local depressions. Such details were not visible on the 1-km resolution map (Fig. 7b).

Lapse rate estimation

Table 4 shows significant seasonal variations in the average lapse rate in land surface night temperature calculated over the study area, with statistically stronger values in the mid-summer period (-0.64 °C/100 m) in comparison with the beginning or the end of the season (ca. -0.60 °C/100 m). The linear relationship between elevation and temperature was also greater in mid-summer compared to the other two periods. The spatio-temporal variability in lapse rates was high since the coefficients of variation were between 24 and 35% in the successive time periods, with lower variation in the midsummer.

Discussion

The present study provides the first application of MODIS products to characterize frost occurrence over a 45000 km² area in the Andean highlands. Frost occurrence over the summer period was either calculated directly from land surface temperature remotely sensed at a 1-km scale, or estimated and downscaled at 100-m by means of physiographic modeling. To our knowledge this is also the first application of BRT in physiographic modeling. Both techniques are complementary in characterizing frost occurrence: remote sensing brings spatialized and repetitive information on land surface temperature and physiographic features, while BRT allow to explore the relative contribution of physiographic factors at various scales and, hence, to downscale satellite information to a level appropriate to farming and land management applications.

Application of remote sensing data for frost occurrence characterization

As pointed by [START_REF] François | Frost risk mapping derived from satellite and surface data over the Bolivian Altiplano[END_REF] at least three factors may potentially affect the relation between Ts and Tn records: a difference in time (ca. 6 a.m. for minimum night air temperature versus 2 a.m. for satellite radiative temperature), a difference in height (1.5 m above the soil surface for meteorological data versus land surface temperature for satellite data), and a difference in spatial resolution (ca. 100 m 2 footprint for local meteorological data versus 1 km 2 for satellite data). In spite of this, these authors observed only a stable shift of some degrees between Tn records in the Bolivian altiplano and Ts registered at 2 a.m. with a 1-km spatial resolution by the NOAA/AVHRR satellite. A similar result was found in the present study showing a linear and highly significant correlation of MODIS land surface temperature at night with minimum air temperature recorded in meteorological stations (R 2 = 0.81).

However, this validation may be biased since, as in most mountain areas in the world, the available meteorological records are likely not representative of the most elevated and isolated parts of the study area. Regarding the MODIS satellite, recent studies have improved and validated its calibration algorithm for land surface temperature in various situations encompassing Bolivian highlands, semiarid and arid regions, or nighttime/daytime overpasses [START_REF] Wan | New refinements and validation of the MODIS Land-Surface Temperature/Emissivity products[END_REF][START_REF] Wang | Validating MODIS land surface temperature products using long-term nighttime ground measurements[END_REF]. This ensures the reliability of MODIS temperature data in the study area in spite of its specific location in cold and arid tropical highlands. After verifying for the normality of the distribution of Ts data, the probability of frost occurence can be easily calculated from the available satellite data series (eq. ( 1)). It should be noted that the available series of 6-year daily records was long enough to statistically characterize frost occurrence at the standard meteorological threshold of 0 °C, but not at lower temperature levels due the scarcity of observations of severe frost events over the considered period. Frosts at -4 or -7 °C would, however, be more relevant for agroclimatic purposes since they correspond to the frost tolerance levels of major Andean crops such as potato and quinoa [START_REF] Bois | Response of some Andean cultivars of quinoa (Chenopodium quinoa Willd.) to temperature: effects on germination, phenology, growth and freezing[END_REF][START_REF] Garcia | Agroclimatic constraints for rainfed agriculture in the Bolivian Altiplano[END_REF][START_REF] Geerts | Agro-climatic suitability mapping for crop production in the Bolivian Altiplano: a case study for quinoa[END_REF][START_REF] Jacobsen | Plant responses of quinoa (Chenopodium quinoa Willd.) to frost at various phenological stages[END_REF]. This limitation should progressively disappear as the MODIS archives grow and allow for the statistical evaluation of less frequent (and more severe) frost events.

Spatial and temporal patterns of frost occurrence 4.2.1. Frost occurrence as affected by regional physiography and climate seasonality

The frost occurrence maps derived from MODIS data at 1-km resolution (Fig. 2) clearly show the influence of regional-scale physiography like the mountain distribution or the proximity of the salt lakes, as well as the seasonal variation of frost occurrence over a 6-month period. In their attempt to map agroclimatic suitability in the Bolivian altiplano, [START_REF] Geerts | Agro-climatic suitability mapping for crop production in the Bolivian Altiplano: a case study for quinoa[END_REF] notice that frost risk is difficult to interpolate spatially. Nevertheless, based on data from 41 ground climatic stations of the altiplano, they achieve a description of regional frost risk patterns that are globally confirmed by our satellite maps, with lower frost probabilities in the Intersalar region and higher probabilities at the south-west of the salt lake of Uyuni. [START_REF] Geerts | Agro-climatic suitability mapping for crop production in the Bolivian Altiplano: a case study for quinoa[END_REF] also mention that their kriging interpolation was improved by incorporating a WNW anisotropy due to the combined north-south influence of the Lake Titicaca and the west-east effect of zonal winds. These zonal winds affect the entire altiplano and largely control the synoptic weather types [START_REF] Garreaud | The climate of the Altiplano: observed current conditions and mechanisms of past changes[END_REF]. As demonstrated in other cold regions or mountain areas in the world [START_REF] Blandford | Seasonal and synoptic variations in near-surface air temperature lapse rates in a mountainous basin[END_REF][START_REF] Dobrowski | How much influence does landscape-scale physiography have on air temperature in a mountain environment[END_REF][START_REF] Marshall | Near-surface-temperature lapse rates on the Prince of Wales Icefield, Ellesmere Island, Canada: implications for regional downscaling of temperature[END_REF], the occurrence and spatial patterns of the zonal winds could be important drivers of the seasonal variation in land surface temperature and temperature lapse rate found in the study area.

Another driver of frost seasonality is cloud cover which, in a typical tropical unimodal rainy season, results in progressive overcasting at the beginning of the rainy season, maximum cloud cover in the mid-season, and then progressive decrease by the end of the season. This cloud cover pattern is recorded daily at the 1-km scale by the MODIS satellite and was integrated in the calculation and mapping of frost occurrence (equation 1, Fig. 2).

The seasonal change in sky cloudiness also influences temperature lapse rates. [START_REF] Blandford | Seasonal and synoptic variations in near-surface air temperature lapse rates in a mountainous basin[END_REF] thoroughly discuss the effect of seasonal and synoptic conditions on lapse rate calculated for average or daily extreme values of near-ground temperature in temperate mountains. They outlined that minimum temperature lapse rate are shallower when air masses are dry and cold, which is explained by increased frequency of cold air drainage and temperature inversions under clear-sky and dry air conditions at night. Minimum temperature lapse rate values are also more variable during seasonal transition between summer and winter, due to fluctuating weather regime at that time and higher frequency of temperature inversions. Our estimates of minimum temperature lapse rate in successive time periods (Table 4) are consistent with both assertions, showing steeper and less variable values in mid-summer (January-February) when the sky is more cloudy and air conditions are relatively humid, temperate, and stable. These estimates approximating -0.6 °C/100 m appear fairly high compared to minimum temperature lapse rate values, typically ranging from -0.15 to -0.35 °C/100 m, in mountains of mid-latitude regions [START_REF] Blandford | Seasonal and synoptic variations in near-surface air temperature lapse rates in a mountainous basin[END_REF][START_REF] Dobrowski | How much influence does landscape-scale physiography have on air temperature in a mountain environment[END_REF][START_REF] Harlow | Derivation of temperature lapse rates in semi-arid south-eastern Arizona[END_REF]. In subtropical mountains however, De Scally (1997 ) states that, due to their high thermal regime, the temperature lapse rate is generally higher than in mid-latitude mountains. The high lapse rate values thus quoted for the Himalaya (De Scally, 1997 ) or the Andes [START_REF] Frère | Estudio agroclimatológico de la zona andina[END_REF][START_REF] Snow | The climates of northern South America[END_REF]cited by Pielke & Mehring, 1977;[START_REF] Trombotto | Monitoring of mountain permafrost in the Central Andes, Cordon del Plata, Mendoza, Argentina[END_REF] refer unfortunately to mean daily or mean annual temperatures which cannot be compared directly to our estimates of minimum temperature lapse rates in specific time periods.

Astronomic forcing is another cause of seasonality in the topography-frost relationship, explaining why topographic controls, usually treated as stationary, show actually pronounced seasonal variations, as pointed by [START_REF] Deng | Multi-scale linkages between topographic attributes and vegetation indices in a mountainous landscape[END_REF] in the case of topography-vegetation relationships. In our study, seasonal changes in the effects of slope steepness and aspect on frost occurrence are illustrated by varying SLO and DPI contributions (Fig. 4b). They are explained by astronomical forcing resulting in insolation values in the early and mid-summer higher by 25% in average than in the late-summer period (see Table 1), thus giving higher influence of DPI over SLO in the former two periods.

Fig. 2 shows that the shores of the salt lakes are less prone to frost, while highlands at the west and south of the region are continuously exposed, even in mid-summer (January-February) when frost occurrence is generally low. This latter situation is obviously due to extreme elevation, while the "milding" effect of the salt lakes could be due to the specific thermal properties of these vast salted extenses. [START_REF] François | Frost risk mapping derived from satellite and surface data over the Bolivian Altiplano[END_REF] observed warmer night temperatures over the Coipasa and Uyuni salt lakes and suggest that water covering these lakes part of the summer, as well as the higher thermal conductivity and thermal capacity of the salted substratum, could explain that their borders remain warmer than the surrounding areas.

From regional to subregional physiographic influences on frost occurrence

Multiple regression methods were used in previous studies to evaluate the relative contribution of topographic factors to near-ground temperature and frost occurrence. These studies were generally conducted in mid-or high-latitude mountains [START_REF] Bennie | Characterising inter-annual variation in the spatial pattern of thermal microclimate in a UK upland using a combined empirical-physical model[END_REF][START_REF] Chuanyan | Methods for modelling of temporal and spatial distribution of air temperature at landscape scale in the southern Qilian mountains, China[END_REF][START_REF] Dobrowski | How much influence does landscape-scale physiography have on air temperature in a mountain environment[END_REF], and often in densely forested areas at mid-altitude [START_REF] Blennow | Modelling minimum air temperature in partially and clear felled forests[END_REF][START_REF] Lindkvist | A frost assessment method for mountainous areas[END_REF][START_REF] Pypker | Using nocturnal cold air drainage flow to monitor ecosystem processes in complex terrain[END_REF]. Our study explored an extended agricultural region at its extreme elevation limit in cold and arid tropical highlands. In this context, BRT models clearly outperformed multiple regression models, probably due to their capacity to include nonlinear effects and interactions between predictors [START_REF] Martin | Optimizing pedotransfer functions for estimating soil bulk density using boosted regression trees[END_REF]. At the regional scale, BRT analyses show that elevation, distance to the salt lakes and latitude were the physiographic features most contributing to frost occurrence variations, while features directly or indirectly related to slope or topographic convergence (SLO, CTI, DPI, API, and MPI) were less important. These regional BRT models based on physiographic features alone explain between 78% and 90% of the variation in frost occurrence observed in different time periods (Table 2).

In their study of the influence of physiography on the distribution of climate variables across the United States, [START_REF] Daly | Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States[END_REF] outline that the effects of elevation and proximity to large water bodies exceed those of other topographic factors at large scales, whereas the effects of slope and landcover features become prominent at relatively smaller scales. In complex terrains, local variations in slope aspect and steepness create a mosaic of hillslopes experiencing contrasting climatic regimes [START_REF] Daly | Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States[END_REF], while topographic depressions are another landscape feature commonly associated with cold air drainage and frost occurrence [START_REF] Lundquist | Automated algorithm for mapping regions of cold-air pooling in complex terrain[END_REF][START_REF] Pypker | Using nocturnal cold air drainage flow to monitor ecosystem processes in complex terrain[END_REF]. These local landform features emerge as forcing factors of frost occurence at the local scale, where the range of variation in elevation and latitude became limited while that in local landscape features remained large. When applied to the reduced spatial domain of the Intersalar, our BRT analyses indeed showed that elevation lost some importance at the benefit of daily potential insolation (DPI), slope steepness (SLO), or topographic convergence (CTI) (Fig. 4b). As a major characteristic of the physiography of southwestern Bolivia, the vast salt lakes of Coipasa and Uyuni remained influential at that local scale, as

shown by the high contribution of the distance to the salt lakes (LDS) in these BRT subregional models.

The good fit of frost occurrence values predicted by BRT models applied either at the 100-m or the 1km resolution (Tables 2 and3) validates the use of BRT regional models for local frost occurrence estimations since these models were able to seize both the influences of large-scale factors like latitude and elevation, and of local factors like slope steepness, insolation, and landscape position. The resulting local mosaic of cold depressions and warmer slopes at particular elevations and exposures is illustrated by the map in Fig. 7c. We hypothesized that small-scale variations in soil warming due to differential insolation in the day before (or the morning after) a given night could influence soil cooling and thus radiative frost at night in particular places. In fact, the contribution of DPI appeared significant in the BRT subregional model, at least from November to February when potential insolation is at its seasonal maximum (Table 1), thus leading to highest contrasts in soil energy balance between sunlit and shaded locations. However, the small contributions of the afternoon or morning components of insolation (API and MPI) (Fig. 4b) seem to belie the idea that potential insolation is directly involved in frost vulnerability at particular places. Actually, [START_REF] Blennow | Modelling minimum air temperature in partially and clear felled forests[END_REF] states that the larger amount of heat stored into the ground in sunlit places cannot compensate for soil cooling at night since this cooling occurs within a few hours after sunset. Nocturnal soil cooling should be still faster under clear-sky conditions at high altitude. This is, however, in contradiction with the common perception of lower frost occurrence in sunlit slopes, particularly in stony terrains and shallow soils supposed to benefit from the thermal stability provided by the rocks. Microclimate stability associated to rock outcrops has been documented by [START_REF] Rada | Low temperature resistance in saplings and ramets of Polylepis sericea in the Venezuelan Andes[END_REF] in the paramo ecosystem of Venezuelian Andes at lower elevation (3800 m) and under wetter conditions (969 mm of annual precipitation). It is likely that the much drier conditions of the puna ecosystem in southern Bolivia reduce the thermal inertia of the soils, thus leading to a very fast soil cooling at night. Apart from astronomical forcing discussed previously, the varying importance of CTI, DPI and SLO in the BRT models, as well as the interactions between them (Fig. 6) reflect complex spatio-temporal relations between insolation and landform factors producing multiplicative or mitigating effects on near-ground temperature. At the microscale level, unobserved soil and vegetation properties might also interfere with landform features. Soil moisture and vegetation cover, for example, are known to influence the radiative balance at the soil surface, and might contribute to buffer the near-ground temperature from cold extremes in particular places [START_REF] Fridley | Downscaling climate over complex terrain: high finescale (<1000 m) spatial variation of near-ground temperatures in a montane forested landscape (Great Smoky Mountains)[END_REF][START_REF] Geiger | The climate near the ground[END_REF].

Ecological implications

The relationships between ecological patterns and processes change across spatial and temporal scales, with singular complexity in mountain areas (e.g., [START_REF] Deng | Multi-scale linkages between topographic attributes and vegetation indices in a mountainous landscape[END_REF][START_REF] Saunders | Hierarchical relationships between landscape structure and temperature in a managed forest landscape[END_REF].

Regarding air or soil surface temperature in mountains, nested factors are interacting, from regional synoptic weather forcing to local topoclimatic situations and microscale variations in vegetation cover and soil moisture. All these factors in turn may dominate the distribution of temperatures, depending not only on the dynamics of the situation (turbulent or stable, nighttime or daytime conditions…) but also on the spatial and temporal scale of interest (from macroscale to microscale, from seasonal to instantaneous). In this way, macroscale conditions of clear sky and calm nights are required for radiative frost to occur, but the frequency and severity of these frost events are further increased by low site position (or conversely, extremely high location) and, at still smaller scales, by vegetation sparseness and soil surface dryness or roughness [START_REF] De Chantal | Frostheaving damage to one-year-old Picea abies seedlings increases with soil horizon depth and canopy gap size[END_REF][START_REF] Fridley | Downscaling climate over complex terrain: high finescale (<1000 m) spatial variation of near-ground temperatures in a montane forested landscape (Great Smoky Mountains)[END_REF][START_REF] Langvall | Effect of shelterwood density on nocturnal near-ground temperature, frost injury risk and budburst date of Norway spruce[END_REF][START_REF] Oke | The temperature profile near the ground on calm clear nights[END_REF]. In the Andean highlands, instantaneous near-ground minimum temperature may be 4°C lower in a sparsely vegetated area compared to a neighboring forest understory [START_REF] Rada | Low temperature resistance in saplings and ramets of Polylepis sericea in the Venezuelan Andes[END_REF]. Similar fine scale variations in minimum air temperature occur within cultivated canopies despite the low plant cover of most Andean crop species (see [START_REF] Winkel | Assessing the protective effect of vertically heterogeneous canopies against radiative frost: the case of quinoa on the Andean Altiplano[END_REF] for the quinoa crop). These local variations in minimum near-ground temperatures may be sufficient for some part of the vegetation to escape lethal freezing. Potential frost impacts on vegetation operating at regional and subregional scales may thus be over-shadowed by microscale variability in minimum

temperature. Yet, contrary to what occurs in dense forests where plant interactions within canopies are significant (Bader et al., 2008;[START_REF] Turnipseed | Airflows and turbulent flux measurements in mountainous terrain: Part 1. Canopy and local effects[END_REF], the sparse and low vegetation typical of the Andean highlands is likely to exert an influence limited to small spatial scales, with topography and coarse scale factors controlling most of the variation in minimum air temperature. In fact, [START_REF] Blennow | Modelling minimum air temperature in partially and clear felled forests[END_REF] outlines that topographic influences on minimum air temperature increase in parallel with decreasing vegetation cover.

Practical implications

The latter consideration implies that agroclimatic applications, such as crop zonation or suitability assessments, require a multi-scale approach, ideally complementing frost risk characterization at the topoclimatic scale by an evaluation of the local effects of crop practices on canopy structure and soil surface moisture and roughness. Though limited to topography-frost relationships, our attempt of downscaling frost occurrence at a 100-m scale usefully expands previous works on regional agroclimatic zoning in the Bolivian altiplano [START_REF] François | Frost risk mapping derived from satellite and surface data over the Bolivian Altiplano[END_REF][START_REF] Geerts | Agro-climatic suitability mapping for crop production in the Bolivian Altiplano: a case study for quinoa[END_REF]. To our knowledge, this is the first time that such a detailed zonation of topoclimate is reported for this region, providing fine-scale information helpful for land management and rural planning [START_REF] Theobald | Ecological support for rural land-use planning[END_REF].

Considering the scarcely available meteorological records in the study area, these 100-m scale maps bring new information about the spatio-temporal variation of frost occurrence, allowing now to localize exactly the seasonal pattern of frost typical of the Andean summer period [START_REF] Frère | Estudio agroclimatológico de la zona andina[END_REF][START_REF] Troll | The Cordilleras of the tropical Americas: aspects of climatic, phytogeographical and agrarian ecology[END_REF]. For instance, the virtual zero value of frost frequency in January and February derived from meteorological records at Salinas (Fig. 1), covers in reality a wide range of situations with still significant frost occurrence in mid-summer, as on the nearby border of the salar of Coipasa or the western and southern part of the study area (Fig. 2). In fact, the recent expansion of quinoa crop in the region was firstly and mostly located in flat areas near the Coipasa and Uyuni salt lakes [START_REF] Vassas | El auge de la quinua: cambios y perspectivas desde una visión social[END_REF], which exemplifies the complex trade-offs between agroclimatic risks and economic expectancies operating in farmers' decision making [START_REF] Luers | The surface of vulnerability: an analytical framework for examining environmental change[END_REF][START_REF] Sadras | Dynamic cropping strategies for risk management in dryland farming systems[END_REF].

Perspectives

Through remote sensing of land surface temperature and modeling of topographic features implemented within a boosted regression procedure, we were able to explicitly downscale frost occurrence at the landscape scale. The method developed here may be adapted to climatic or ecological processes other than frost. Rainfall distribution could be a candidate since its spatio-temporal patterns clearly depends on landscape characteristics over complex terrains. Current litterature outlines the importance of landform as a factor of rainfall variability in the Andes [START_REF] Giovannettone | Probing regional orographic controls of precipitation and cloudiness in the central Andes using satellite data[END_REF], though most studies were conducted at the coarse spatial resolution appropriate to continental scale climatology [START_REF] Garreaud | Interannual rainfall variability over the South American altiplano[END_REF][START_REF] Misra | Dynamic downscaling of seasonal simulations over South America[END_REF][START_REF] Vuille | 20th century climate change in the tropical Andes: observations and model results[END_REF]. Canopy energy budget, soil water balance, or ecosystem productivity are other ecological processes tractable for topographic modeling [START_REF] Bradford | The relative importance of light-use efficiency modifications from environmental conditions and cultivation for estimation of large-scale net primary productivity[END_REF][START_REF] Rana | Estimating energy fluxes from sloping crops using standard agrometeorological measurements and topography[END_REF][START_REF] Urban | Forest gradient response in Sierran landscapes: the physical template[END_REF]. Such applications depend firstly on the availability of remotely sensed proxies for the considered process. For example, the remotely sensed daily amplitude in surface temperature and vegetation indices can be used to derive daily evapotranspiration [START_REF] Wang | Estimation of evaporative fraction from a combination of day and night land surface temperatures and NDVI: a new method to determine the Priestley-Taylor parameter[END_REF]. Similarly, satellite estimates of absorbed photosynthetically active radiation may serve to evaluate net primary productivity [START_REF] Bradford | The relative importance of light-use efficiency modifications from environmental conditions and cultivation for estimation of large-scale net primary productivity[END_REF][START_REF] Turner | Assessing FPAR source and parameter optimization scheme in application of a diagnostic carbon flux model[END_REF]. An additional requisite for the calibration of these applications consists in local ground measurements for the variable of interest. This is a major issue in the case of the tropical highlands where, similarly to what occurs for meteorological data, reliable datasets on matter and energy fluxes at ground level are and will remain scarce [START_REF] Vergara | Visualizing Future Climate in Latin America: results from the application of the Earth Simulator[END_REF]. The methods and results presented here can contribute to a better understanding of the potential risks associated with climate and land use changes in complex terrains, so that decision-makers can develop efficient strategies to improve the ecological sustainability of natural and agricultural ecosystems in vulnerable mountain areas.

at

  Salinas de Garci Mendoza (19°38'S, 67°40'W) managed by the SENAMHI (Meteorology and Hydrology National Service, Bolivia) where daily minimum air temperature (Tn) was recorded in 1989 and from 1998 to 2006, and two others at Irpani (19°45'S, 67°41'W) and Jirira (19°51'S, 67°34'W) where meteorological stations set up by the IRD (Research Institute for Development, France) recorded semi-hourly air temperature from 23 November 2005 to 18 February 2006 in Irpani, and from 6

Fig. 1 .

 1 Fig. 1. Frequency analysis of daily minimum air temperatures lower than 0 °C registered at Salinas in 1989 and from 1998 to 2006.

Fig. 2 .

 2 Fig. 2. One-kilometer resolution maps of MODIS-derived frost occurrence in southern Bolivia in three successive time periods (the color scale at the right shows the frost probability) During the austral summer (November-April), two periods of frequent below zero temperatures surround a ca. 80-day time interval of low frost occurrence, from the beginning of January to the end of March.Daily Tn values recorded at screen height by the meteorological stations of Salinas, Irpani and Jirira were highly correlated to Ts data remotely sensed over these three localities at night by the MODIS satellite (Tn = 0.97*Ts + 0.93, R 2 = 0.81, n = 750). The percentage of cloudy pixels on the MODIS images gives a general information about the seasonal pattern of cloud cover in the study area: the January-February period was the most overcast with on average 46 ± 8.2% of the study area masked by clouds in each daily satellite image, while this percentage fell to 30 ± 7.4% and 29 ± 6.6% in the November-December and March-April periods respectively. The test of Kolmogorov-Smirnov applied to

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison of frost occurrence values observed in three time periods with frost occurrence predicted by the MR (multiple regression) and BRT regional models (ND = November-December; JF = January-February; MA = March-April, n = 49353 in each time period)

Fig. 5 .

 5 Fig. 5. Partial dependence plots of the six most influential physiographic variables in the regional BRT in three time periods (continuous line: November-December, dashed line: January-February, dotted line: March-April; ticks at the inside top of the plots show deciles of site distribution across the variable;

Fig. 7 .

 7 Fig. 7. Elevation map of the Mount Tunupa area (a), and frost occurrence in the March-April period mapped at 1-km resolution from MODIS observations (b), at 100-m resolution (c) and in 3-D view using regional BRT (d). Frost occurrence is scaled between 0 and 1 as the probability of daily occurrence of negative Ts values in the March-April period.

Table 1 .

 1 Ranges of physiographic variables observed over the study area.

	Variable	Minimum	Maximum	Unit
	LAT latitude in UTM 19 South	-22.00	-19.24 decimal degree
	ELE elevation	3540	6051	m
	SLO slope steepness	0	35	degree
	DPI daily potential insolation
	Nov-Dec	6900	9075	W m -²
	Jan-Feb	7040	8935	W m -²
	Mar-Apr	4604	7897	W m -²
	MPI morning potential insolation
	Nov-Dec	510	1861	W m -²
	Jan-Feb	466	1772	W m -²
	Mar-Apr	264	1240	W m -²
	API afternoon potential insolation
	Nov-Dec	2819	5075	W m -²
	Jan-Feb	2816	5000	W m -²
	Mar-Apr	2079	4308	W m -²
	LDS distance from salt lakes	0	5.04	Ln (km + 1)
	CTI compound topographic index	5.7	14.1	-

Table 1

 1

	for variables abbreviations)

Table 3 .

 3 Statistical comparison of 100-m frost occurrence predictions reaggregated at 1-km with

	observed 1-km frost occurrence values (n = 49353). B: bias; RMSE: root mean square error of
	prediction; R2: determination coefficient of the regression line between observed and predicted values.
	ND = November-December; JF = January-February; MA = March-April.
	Period	B	RMSE	R 2
	ND	0.0252	0.0925	0.74
	JF	0.0008	0.0644	0.87
	MA	0.0543	0.0945	0.80

Table 4 .

 4 Descriptive statistics of lapse rates of land surface temperature at night in three successive time periods.

	Period	Mean	Coefficient of	Coefficient of	Sample size
		(°C/100 m)	variation (%)	determination
	November -December	-0.609	35.5	0.50	366
	January -February	-0.642	24.3	0.68	355
	March -April	-0.598	29.8	0.56	361

Acknowledgements

The authors are most grateful to Danny Lo Seen (Cirad, France) for discussions on BRT methods.

They also aknowledge the comments made by the anonymous reviewers of the paper. This work was carried out with the financial support of the ANR (Agence Nationale de la Recherche -The French National Research Agency) under the Programme "Agriculture et Développement Durable", project "ANR-06-PADD-011-EQUECO".