Introduction
The numerical parameterization of all processes which are involved in snow-atmosphere exchanges is one of the most promising ways to follow the state of snow cover accurately. With this aim the Centre d'Etudes de la Neige (CEN) has developed and validated a model called CROCUS (Reference Brun, Martin, Simon, Gendre and ColéouBrun and others, 1989 and 1992). This model is now being used for operational avalanche forecasting as well as for research and studies. The main purpose of the model is to simulate in real time the snow-cover evolution on numerous locations at different elevations, on different aspects and in different massifs of the French Alps. However, the model needs hourly meteorological information on the snow-atmosphere interface at every simulation point, most of which are not near meteorological stations. Therefore a meteorological application of objective analysis (Reference Brun, Durand, Guyomarc’h and MerindolBrun and others, 1991) was developed. (In this context objective means automatic, as opposed to subjective which means manual.) Its name is SAFRAN and its aim is to deduce the relevant parameters needed by the snow model at numerous given points by using all data available around these analysis points. SAFRAN’s specifications are to provide hourly data automatically for 300 m elevation increments on six aspects (Ν, E, W, S, SE, SW) of the 23 principal massifs in the French Alps. The weather data are for surface air temperature, wind speed, air humidity, cloudiness, precipitation, precipitation type, atmospheric radiation, and direct and scattered solar radiation.
The main assumption is that the homogeneity of each massif allows spatial interpolations for both altitude and aspect in the final analyzed result. However this does not imply that every massif is isolated; on the contrary, there are statistically modeled exchanges of information between all involved points (analyzed and observed) by means of different statistical parameters (variances and correlations) used in the analysis method described below.
The “real” spatial size of each massif is now rather variable, for they are the same as those used by the French for operational avalanche risk forecasting. As some are somewhat too large (e.g., Vanoise), it is planned to divide some of them to get homogeneous sizes of less than 1000 km2.
Our particular analysis problem is in fact very similar to the more general problem of the initialization of numerical weather prediction models, and we shall solve it with similar methods described by Reference GandinGandin (1963) and still in operational use (Reference DurandDurand, 1985).
Available Observation Data
A great variety of available observational data are used. These include the standard World Meteorological Organization (WMO) meteorological messages such as SYNOP or SYNOR (automatic but reduced version) from ground-based observations. There are about six SYNOP observatories in the vicinity of our massifs that are analyzed; they are situated at low altitude, but they provide a complete set of standard observation data at four times daily: 0,6,12,18 UTC. At higher altitudes we take the TEMP or PILOT messages for upper-air information (three or four elevations twice daily: 0,12 UTC) which are also very useful in those cases where outputs of the model are lacking. The French Met-Office and many ski resorts take part in an ancillary network of surface observations (NIVO-METEO) which provides pertinent information twice daily (6,12 UTC) during the winter ski period. We have about 100 observations of this kind over our total area (23 massifs, about 20 000 km2), but these data do not completely cover all aspects. For instance, the southern aspect has very sparse information and the altitude range is quite small (1000 – 2200m a.s.l). Other data sources from the climatological network can be used: some of its data can be obtained in real time when the NIVO-METEO network is not working, and the network’s daily precipitation measurements are in this case our only reference in middle-mountain areas.
Many other observed parameters are used for crosschecking and diagnostics concerning the accuracy or the time-consistency of the meteorological events, especially for precipitation and the time-sequentiality of their phases. If possible, we build bogus (not observed) data of partial cloudiness by combining the different elements of the SYNOP and NIVO-METEO messages as the “observed present weather”; the result is then adjusted to match the observed global cloudiness. Some similar processes are also used to make altitude profiles of relative humidity using the same observed elements. The 24 h observed precipitations are corrected with the observations of freshly fallen snow depth according to an estimation of snow density based on wind speed and air temperature (personal communication from E. Pahaut).
Every potential observation point has been examined and its representativity evaluated by avalanche fore-casters. Coefficients corresponding to the different aspects are applied to the data. Therefore every observation site can be used in the analysis with a different level of confidence according to its internal observation error and its relationship with the analyzed aspect.
Analysis Methods
Analysis Points
The massifs that we want to analyze are represented in our model by irregular pyramids whose sides (corresponding to the different aspects) have a precise geographical location. On each aspect of each massif, we first analyze a vertical profile of altitude variables (vertical increment of 300 m) which gives an estimate of the free atmosphere characteristics; we then analyze the surface air variables and precipitation at the same vertical location as the altitude variables. The analyzed altitude variables are used to compute the different incoming radiation components.
Altitude
The altitude parameters (air temperature, wind, humidity and cloudiness) are analyzed by the Optimal Analysis (OI) method which has been widely used for many years in meteorological cases and is decribed in a practical case in Reference DurandDurand (1985). It is in fact a linear regression which uses a flexible statistical scheme. All the correlations are homogeneously modeled by adapted functions which allow an estimate of every analyzed parameter with different possible adjacent data. To estimate the quantity Y a at location k
The number of observed parameters used is Ν and the quantities
An estimation of the variance of the analysis error of
Surface
The surface scheme is cruder, for it is performed after the application of the altitude scheme. Its own guess values can be first updated in order to take into account the information of the altitude analysis. The same formula as Equation (1) is used but the computation of Ρ is simpler and is only dependent upon the distance between the points i and k. (It is inspired by the formulation of the (SC) is iterative in order to treat different scales: the result after a scan is used as a guess of the following scan. In the final scan (the fifth in our case) the observations only influence their closest analysis points.
Precipitation
The analysis is performed only once a day because the 24 h precipitation is the most commonly observed parameter, especially by the NIVO-METEO network. We still use the same OI analysis package as in the altitude analysis but with different tunings: we can thus select observations farther away when the NIVO-METEO network is not available and there is no vertical correlation. An additional constraint is put in the linear system of Equation (2) which assures that the analysis weights P i are linearly combined by the guess value at the observation points
Guess-Field
The interpolation methods used here require a preliminary estimate of the desired quantities. This guess represents the best estimation of the final fields before taking into account the observed data, which will add supplementary information. For the temperature, wind and humidity altitude variables, this guess generally comes from an extraction of the French mesoscale-forecast model PERIDOT fields (grid size of about 35 km, 15 vertical sigma levels) (Reference ImbardImbard and others, 1987) with different scale operators in order to take into consideration the initial smoothed orography and our finer working scale. The variables in question are then interpolated horizontally by using Bichkard polynomials and vertically on analysis locations. Sometimes we have to estimate an extrapolated vertical structure under the orography of the meteorological model when it is higher than the desired analyzed altitude. If the PERIDOT model is lacking, it is still possible to use climatological profiles. In that case these profiles are first blended with the observations in order to be improved, but the final result is of low quality and our analysis problem is quite underdetermined.
The surface-boundary effects computed by the meteorological model are used to determine the preliminary surface values for temperature, wind and humidity. They are used as vertical differences (between surface and altitude values at the same altitude) corrected by the height and the wind speed in order to assure a good transition to the “free atmosphere” conditions from the surface values at the same elevations.
The guess cloudiness is diagnosed from the humidity-analyzed profile results. Three vertical cloud classes are computed when the vertically averaged humidity in the cloud layer is greater than an altitude-dependent critical value (Reference ImbardImbard and others, 1987). These three vertical classes correspond to the WMO cloud classification at low, medium and high altitude.
At the locations of the analyzed points the guess precipitation are obtained from climatologically analyzed fields which are available on a 5 km grid and which take into account orographic effects (Reference Bénichou and Le BretonBénichou and Le Breton, 1987). “Real” climatological values are computed at the observation points by using the observed series available in the CEN database. All of these computations are done according to seven different typical weather conditions deduced from the shape of the 500 hPa geopotential field making the seven different possible precipitation guess-fields spatially consistent with the observed series. Just before every analysis, the guess-field for each massif is slightly modified using observed values to fit daily meteorological conditions better. These conditions can exhibit larger differences between massifs than the corresponding smoothed guess-fields. For this purpose we keep the vertical guess gradients but the mean massif guess value is set to the observed one.
Analyses and Interpolation Processes
Temperature, Wind, Humidity, Cloudiness Analyses
These analyses are run each of the four synoptic hours because most of the data and the PERIDOT outputs are available only at this frequency. All previously mentioned observation types can be used, but the surface observations have weaker weights than the altitude data derived from the meteorological air soundings. The selection of the Ν elementary variables for the linear systems of Equation (2) is based on geographical and statistical criteria in order to assure that each analysis point is well surrounded with informative observations which can be at different altitudes and massifs. Currendy the maximal value of Ν is 13 but at some locations this number cannot be reached because of too few observations. All these selections are completely automatic and in the worst cases the system can run with very few data, providing a result of poorer quality and very close to the guess results. The cloudiness analysis is performed only on the three WMO layers.
The modifications of the altitude profiles are then used to modify the corresponding preliminary surface guess-field values before running the surface analysis; this brings some additional information and assures consistency between altitude and surface values. This second analysis step concerns only surface parameters, and is done with the SC scheme. Only mountainous ground-based observations are used and the single observations of wind speed of many automatic stations are completed with the guess-wind direction in order to be inserted in the surface-wind components analysis.
In order to eliminate automatically all erroneous data prior to their use in the different analysis schemes, all observed parameters are previously checked for their quality. This is done by performing an OI analysis at their own location without using them and comparing the obtained values to the observed one with a criteria based on the variance of the analysis error. If this value is low, it means that the analysis process has found partially redundant information, the analysis result is almost sure, and the test is rigorous; in the opposite case there is more laxity.
Time Interpolation
In order to correctly feed CROCUS, we must interpolate all our previous analysis results hourly on every analysis location. All quantities are available every 6 h except the precipitation result, which is daily. An hourly distribution function of precipitation is computed from the interpolated hourly specific humidity. It gives a probability of rain when it is greater than a threshold value corresponding roughly to 90% relative humidity. We can thus split our total 24 hours precipitation analyzed value into 24 hourly parts, but without having determined the phase at each hour. This method is quite rough, but it will be ameliorated with radar information. It is completed by an estimation of cloud heights based on observed averages, which helps to refine the vertical distribution of the three analyzed cloud layers.
All the analyzed altitude profiles and the surface wind are then linearly interpolated for each hour and the obtained profiles are used as input to a radiative model which provides the required hourly radiation values (Reference Ritter and GeleynRitter and Geleyn, 1992). The 12 UTC air surface temperature is slightly corrected, especially at the south aspect, by using an analyzed information relative to the daily maximum temperature and weighted by the analyzed cloudiness. For this purpose an analysis of the daily maximum temperature is performed by using the SC method with observed data and the 12 UTC air surface temperature analysis as guess. The idea is to force on these aspects a coherent diurnal effect. The main reason for this modification is to compensate for the clear weakness of the observation network at this aspect which leads to some biased values often due to an insufficient correction of the PERIDOT guess values on these aspects.
Once done, the surface temperature on every aspect-altitude is adjusted hourly by taking into account the forcing based on the previously computed solar radiation and a relaxation to an equilibrium temperature (Reference Martin and MainguyMartin and Mainguy, 1988). This method suggests two differential equations (day and night) in order to set these budgets which are solved with the five analyzed temperatures (6,12,18,0, 6UTC); it yields a complete analytical hourly solution fitting of the analyzed values. Some corrections to the surface humidity are then made in order to keep the air water content constant.
One of the main problems is the hourly determination of the precipitation phase. It is solved in three possible steps in each massif:
(1) Determination of the highest hourly zero isotherm in the temperature/altitude profiles. In the case of air temperature inversion with negative temperatures at low altitudes, the initial value is kept.
(2) Use of the past and present observed weather conditions at analysis times in combination with the analyzed vertical temperature gradient. We thus determine an averaged estimation of the snow-rain transition altitude. The previous estimation is used as initial value and control. We impose the final result to be in a safe vertical interval with limits determined by the altitudes of the different observation sites of the massifs observing rainfall or snowfall.
(3) At every observation site, the computation of the ratio of the equivalent water content of freshly fallen snow compared to the total water content received; the snow density is always computed with the Pahaut’s formula. For each site we evaluate also the sequence of the rainfall and snowfall events (i.e., temporal changes in precipitation type). This ratio is then combined with the results of the previous steps in order to adjust its vertical variation and is modeled by a bounded function which gives a value between 0 and 1 for each massif. This way we get a continual vertical indication of the proportion of the two events in the altitude range where both snow and rain occurred during the day. This altitude-dependent ratio is then imposed on each analysis altitude to the sequences of hourly interpolated precipitation.
Validation
Manual Meteorological Monitoring
Since the 1990–91 winter the results have been often verified by local forecasters of Meteo-France. Their monitoring and the following remarks allowed us to correct some preliminary errors and to implement and tune different algorithms such as the determination of the vertical distribution of snow—rain mentioned above. We can summarize in Figure 1 all the output variables of the analysis as they are sent to the snow model and are daily and automatically available for monitoring.
Fig. 1. Example of the results obtained on the northern aspect of the Vanoise massif on 20 June 1992 at 12 UTC. Surface humidity.... (on left) top left scale in %; wind speed → bottom left scale in ms−1; temperature — bottom scale in °C; cloudiness.... (in center) top scale in tenths. Precipitation is to right of central axis (top right scale
mean snow). All radiation values are plotted on right (direct, 
; scattered, 
; atmospheric, 
with their scales on top and bottom of right side in W m−2.
Objective Comparison
Two well-instrumented automatic observation sites, Col de Porte (CdP, 1340 m) and Col du Lac Blanc (CIB, 2800 m), are deliberately not included in the analysis system so that an objective comparisons can be made without any bias. Tables 1 and 2 show the comparisons during the winter of 1990–91 between all the observed values (known measurement errors and malfunctions in the instruments having been excluded) and those obtained in the analysis system after a double interpolation in height and aspect. Some general comments can be made:
Table 1 Objective comparisons on observed and analyzed quantities at Col de Porte (CdP, 1340m, fat) hourly (h), or at 12 UTC (hl2) or daily (d). Shown are means for observed and analyzed values, their difference, the standard deviation of their difference and the correlation. T, temperature; V, wind speed; H, relative humidity; IR, atmospheric radiation; So, global solar radiation; N, cloudiness; RR, precipitation
Total precipitation observed: 346 mm on 133 common days. Total precipitation analyzed: 270 mm on 133 common days.
Table 2 Same as Table 1 for Col du Lac Blanc (ClB, 2800m, flat).
Total precipitation observed: 297 mm on 103 common days. Total precipitation analyzed: 296 mm on 103 common days.
The observed cloudiness is roughly diagnosed with the direct sun radiation, and this can lead to computation problems, especially at dawn or twilight.
The daily precipitation values are quite well shown at CdP but their cumulative heights are poorly estimated because the site is less representative of its massif and suffers heavier precipitation than the rest of the massif. The results are opposite at CIB. Some current tests (not shown here) show better results at CdP after a more efficient adjustment between the precipitation guess-field and the NIVO-METEO data (see GUESS-FIELD).
The temperature seems to be the best parameter to evaluate; this is not the case with the wind in which the small scales are very badly represented by our scheme. The wind error originates mainly from the 0 and 18 UTC analysis times where no NIVO-METEO observation is available and where the PERIDOT model shows a large positive bias which is not corrected by data.
The humidity shows an important bias which is also variable with the altitude. This has some implications in the cloudiness determination.
The atmospheric radiation integrates this problem of cloudiness and humidity.
Preliminary Operational Runs
A completely automatic application of SAFRAN, CROCUS and some other applications was initialized on 1 October 1991 and has run daily since then without re-initialization (see Fig. 1). Its results have been successfully used during the 1991–92 winter, especially during the Olympic Games of Albertville, by the avalanche risk forecasters.
Climatological Validations
In order to validate the future operational system and also to verify the potential of the system as a climatological research tool, we studied its ability to reproduce the present snow climatology and its variability during the last ten years in the Alps. We took only the snow depth as the verification criterion (because it integrates several phenomena such as accumulation, settling, and melting). SAFRAN used as guess-field the European Center for Medium Range Weather Forecast (ECMWF) analyses instead of PERIDOT because its output fields are not stored. Some climatological stations (a few of them also NIVO-METEO sites, but all others at low elevation in mountainous areas) were also integrated into the system in order to supply precipitation measurements out of the winter season. CROCUS computed snow depths on 37 locations where snow depth measurements were available over the last ten years. Figure 2 shows the simulated snow depths over ten years for the ski resort of Tignes (2080 m) and corresponding observations. These results do not take into account snow drift by wind and localized small-scale phenomena: a bad estimation of snow depth can result where the micro-climatology of a site is different from the climatology of the massif. Furthermore the southern Alps are less well represented due to a lack of data and a more significant variability.
Fig. 2. Fig. 2. Ten years of observed (dotted line) and simulated (solid line) snow depths at Tignes (2080 m).
A recent simulation (not shown here) on these 37 locations has been performed using incomplete data. The results showed good stability; the modeled snow cover at the sites where data was incomplete had very similar characteristics to that of the previous complete run.
Conclusion
We have developed an automatic system for avalanche forecast which will soon be used operationally in the French avalanche forecasting services. Modeling is one of the most promising approaches for many other applications. This kind of system will soon be used as a climatologie tool to evaluate the response and sensitivity of the snow cover to different climatic change scenarios.
Some new developments are under study, such as the insertion of radar pictures and satellite data. The radar would provide an estimate of the time sequence of precipitation events, while the satellite would carry supplementary information on cloudiness.
