Abstract

A numerical forecasting system is presented for automatic prediction of slippery road conditions at road station sites in Denmark. The system is based on a road conditions model forced by input from an operational atmospheric limited area model. Synoptic information on cloud cover and observations of temperature and humidity from the road station sites are taken into account in a sophisticated initialization procedure involving flux corrections for individual locations. The system is run operationally at the Danish Meteorological Institute with promising results. Currently, new forecasts 5 h ahead are produced every hour for 200 road station sites. The potential of the model system is illustrated by model forecasts from two selected periods and operational results. Problems related to the accurate forecasting of road conditions are discussed.

Introduction

An accurate prediction of meteorological parameters such as precipitation, temperature, and humidity close to the ground is of great importance for various applications. For example, warnings about slippery road conditions may be issued if snow, freezing rain, or rime can be forecast with sufficient accuracy. In addition, traffic delays and the risk of accidents may be significantly reduced by specific actions such as road salting.

An impressive amount of money is spent on winter road maintenance in many European countries. For example, it is estimated (Thornes 1996) that the total budget for winter road maintenance in the United Kingdom is about $200 million every year. For Denmark, being a smaller country, the corresponding budget is about half of this. The variability from year to year, however, is considerable.

Unnecessary road salting should be avoided for economic reasons and due to the risk of environmental damage. This means that optimal salting procedures should be sought. In this context, accurate road weather information is vital, which justifies the efforts that are spent on the development of advanced road conditions models. The present paper concerns the development of a numerical model system for operational forecasting of the road conditions in Denmark.

The prediction of the road conditions requires the production of accurate forecasts of temperature, humidity, and precipitation at the road surface. To provide this information, two strategies are possible. The first one relies on the manual work of a forecaster who issues warnings of slippery road conditions based on various meteorological tools—for example, synoptic observations and output from atmospheric models. The second possibility is based on automatic output from specialized models, which may involve statistical or deterministic methods. The two approaches can be used in combination if a forecaster supplies certain input data for the automatic system. This approach has been used in several operational environments—for example, using the U.K. Meteorological Office road surface temperature model (Rayer 1987) and in the Swedish Road Weather Information System (Bogren and Gustafsson1994). In addition, a forecaster may supply specific comments before the products are issued to the authorities responsible for road maintenance.

A manual nowcasting system for the prediction of slippery roads had been operational in Denmark until the end of the winter season 1993–94. This forecasting system was based on the monitoring of data from more than 200 road-weather stations throughout the country. New data—for example, measurements of air temperature and humidity at a level of 2 m—and road sensor measurements of surface temperature were submitted approximately every 10 min from the road-weather stations to the road masters and to the weather forecasters at the Danish Meteorological Institute (DMI). A weather forecaster issued every second hour short-range forecasts valid 3 h ahead. The forecasts were worked out on the basis of road weather data from all stations, together with the latest meteorological information available, for example, weather maps and satellite images. A number of parameters were forecasted, that is, road surface temperature, dewpoint, and air temperature (at 2-m height), cloudiness, precipitation, and wind velocity. The manual forecasts had to cover 13 regions of the country. In some weather situations this implied the production of individual forecasts for every region, which was quite a job for the weather forecaster.

Therefore DMI decided in 1992 to develop an automatic system for prediction of slippery roads, with financial support from the Danish road authorities. A first version of a road conditions model (RCM) has been described by Sass (1992). It is based on the solution of the heat budget equation at the interface between the atmosphere and the road. In an operational context, the model requires input of data from an atmospheric model and synoptic observations including road surface data. The fluxes of radiation, sensible heat, latent heat, and precipitation are computed from operational atmospheric forecasts of temperature, humidity, cloud cover, wind, and precipitation for the geographical locations of the road stations. A special feature of the model system is the advanced treatment of initial conditions including surface flux corrections for each station site. The atmospheric model supplying input to the RCM is the High Resolution Limited Area Model (HIRLAM) at DMI (Sass 1994; Källen 1996).

It is important to realize that the prediction of local road conditions such as road surface temperature and ice requires a careful treatment of initial conditions for the model computations and a realistic forcing from the atmospheric model providing input to the RCM. It appears to be a major problem that current atmospheric models in operational use do not analyze enough local features necessary for precise predictions. In addition, the models are not run with sufficient horizontal resolution (grid point density) to describe adequately the advection of temperature and humidity. At DMI the current model is run with a horizontal mesh size of approximately 22 km. It is well established from observational studies (e.g., Bogren and Gustafsson 1991) that local topography such as terrain and vegetation is responsible for variations in atmospheric variables down to a scale below 1 km. To obtain sufficiently accurate predictions it is thus relevant to include various local observational data in data assimilation. The data should be used to determine realistic initial conditions for the RCM forecasts. Also, the input from the atmospheric model should be corrected if it is inconsistent with observations of important parameters such as cloud cover.

A model system was tested semioperationally during the winter 1993–94 and became fully operational in 1995 for almost all (200) of the road station sites in Denmark. New forecasts are produced every hour and the forecast length is 5 h. The forecast variablesavailablein the manually based system are also produced by the automatic system and can be presented graphically. In addition, text describing the range of the forecast variables and their temporal evolution is generated automatically for the group of stations in each of the 13 geographical regions of Denmark. A weather forecaster at DMI is allowed to supply additional comments on the results of the automatic forecasts before the text, together with forecast data from individual stations, is sent to the road authorities.

The automatic forecasting system is described in the following sections (2 and 3). Section 2 is devoted to aspects related to data assimilation, and section 3 describes the individual components of the RCM forecast. The performance of the RCM system is addressed in section 4 describing experiments based on data from two selected periods. Also the performance of operational runs is discussed. Finally, some concluding remarks are presented in section 5.

Data assimilation

The RCM forecasting system requires input data from an atmospheric model and various observational data that need to be assimilated before the RCM forecast can operate properly. The different parts of the data processing are illustrated in Fig. 1. The data assimilation produces a model state at the forecast initial time tf and atmospheric input data to be supplied as a forcing to the RCM during the forecast. The atmospheric input data (HIRLAM data) are modified by observations, partly from the road station sites. The different observations and variables are explained in Fig. 1. The modified atmospheric data together with the primary model-produced data—that is, road temperatures, road water, and ice—may be considered as the complete model state since all these data contain prognostic information.

Fig. 1.

Overview of the data processing in the RCM forecasting system. The diagram shows that various data are supplied during data assimilation from time td to tf and during the RCM forecast from time tf to tf + L. HIRLAM data (input): T(k) (kN) is temperature at the atmospheric model levels (°C), q(k) (kN) specific humidity at the atmospheric model levels (kg kg−1), T2 m temperature at 2-m height (°C), q2 m specific humidity at 2-m height (kg kg−1), V10 m wind speed (m s−1) at 10 m height, Q precipitation intensity [kg (m2 s)−1], ps surface pressure (Pa), and N number of atmospheric model levels. Synoptic data (input): Cob is observed total fractional cloud cover and Hob observed cloud-base height (m). Road station data (input): T2 m is observed temperature at 2-m height (°C), Td2 m observed dewpoint temperature at 2 m height (°C), and Ts(1) observed road surface temperature (°C). Additional variables (output): Ts(l) (lM) is road vertical temperature profile (°C), Ws road surface water (kg m−2), Is road surface ice (kg m−2), M number of model levels of the RCM, T10 m temperature at 10 m height from T2 m and T(N), q10 m specific humidity at 10 m height from q2 m and q(N), and C(k) fractional cloud cover at atmospheric model levels.

Fig. 1.

Overview of the data processing in the RCM forecasting system. The diagram shows that various data are supplied during data assimilation from time td to tf and during the RCM forecast from time tf to tf + L. HIRLAM data (input): T(k) (kN) is temperature at the atmospheric model levels (°C), q(k) (kN) specific humidity at the atmospheric model levels (kg kg−1), T2 m temperature at 2-m height (°C), q2 m specific humidity at 2-m height (kg kg−1), V10 m wind speed (m s−1) at 10 m height, Q precipitation intensity [kg (m2 s)−1], ps surface pressure (Pa), and N number of atmospheric model levels. Synoptic data (input): Cob is observed total fractional cloud cover and Hob observed cloud-base height (m). Road station data (input): T2 m is observed temperature at 2-m height (°C), Td2 m observed dewpoint temperature at 2 m height (°C), and Ts(1) observed road surface temperature (°C). Additional variables (output): Ts(l) (lM) is road vertical temperature profile (°C), Ws road surface water (kg m−2), Is road surface ice (kg m−2), M number of model levels of the RCM, T10 m temperature at 10 m height from T2 m and T(N), q10 m specific humidity at 10 m height from q2 m and q(N), and C(k) fractional cloud cover at atmospheric model levels.

Surface conditions

The temperature profile in the road is controlled through a careful initialization procedure. The equation of heat conduction in the road (see section 3) is solved during a 3-h data assimilation period prior to the forecast initial time tf. Observed surface temperatures (Ts) from road sensors are imposed as an upper boundary condition during this period to produce a realistic road temperature profile from observed surface temperatures and a temperature profile, which is initially obtained from the previous data assimilation cycle (RCM state 2 of Fig. 1). To reduce a possible drift of the predicted road surface temperature away from observed values in the early phase of the subsequent forecast, the following procedure is applied. A number of short simulations starting at time ta (see Fig. 1) are made for the last period (tfta = 20 min) of the 3-h data assimilation in order to reproduce the observed change of surface temperature during this short period. The forecast module described in section 3 is used, with an atmospheric forcing valid at the forecast initial time, plus a flux correction. The station-dependent flux correction necessary to reproduce the observed surface temperature change during this period is then applied during the subsequent forecast according to (1)

 
formula

−80 ≤ ϕ̂0 ≤ 80 W m−2 (see section 3). It is sufficient to use a flux increment of 5W m−2 to define a suitable number of simulations. The coefficient a1 is currently given a constant value (see appendix), allowing for incorporation of local phenomena declining during the forecast.

In addition, the initial values of surface water (Ws) and road ice (Is) need to be specified. These values are difficult to assess due to lack of reliable observations. Currently Ws is computed as a linearly varying function of the precipitation intensity Q, reaching a maximum value of 0.5 kg m−2 at a precipitation intensity of 0.5 mm h−1. If the observed road surface temperature is below the freezing point, only ice (snow) may be present initially.

Atmospheric input data

It is important to supply realistic input data valid for the atmosphere above the road station at the initial time tf and during the forecast. The original hourly HIRLAM data are interpolated linearly in the horizontal to values representing the atmospheric variables above each road station site. A modified state, denoted by primes, needs to be determined in order to improve the accuracy of the atmospheric data by means of observations.

The primary atmospheric variables are temperature and specific humidity from which a first estimate of cloud cover C (k, j) at the vertical level k and at time level j is made according to (2)

 
formula

This type of formula is often applied to represent cloud cover in atmospheric models, for example, Slingo (1987). In (2), RH is the relative humidity and RHc is a threshold relative humidity, which is 0.80 above a fixed height equal to 1000 m and increases linearly toward 1 at the ground. The cloud cover is zero if the relative humidity is below the threshold value. The first estimate according to (2) may be modified as cloud cover observations and convection are taken into account. The total cloud cover Cob and cloud-base height Hob are analyzed from (3), which is a weighted sum over the available Danish synoptic cloud cover observations:

 
formula

In (3), X may be either total cloud cover or cloud base height, I is the number of observations, Di is the distance between the road station and the site of the cloud cover observation, and D00 (=5 × 104 m) is an associated scaling distance. The exponential functions make it possible for the small-scale features of cloud cover, on a scale below 50 km, to be taken into account since the analyzed value will be dominated strongly by observations close to the road station. The fractional cloud cover C(k, j) derived in all vertical atmospheric levels and all time levels is computed as follows. The observed cloud cover Cob is inserted at the vertical level K closest to the observed cloud base (the lowest clouds observed). At other time levels decreasing weight should be put on the observation and increasing weight on the information from the atmospheric model. An exponential decay from analyzed cloud cover toward model cloud cover has been chosen, according to (4),

 
formula

where tj is the current time in the forecast and tob is the time of theobservation.The value of a2 (see appendix) should reflect the lifetime of clouds but also the accuracy of the atmospheric model. Hence, an optimal value depends on model tuning.

The cloud covers at other vertical levels than K are initially reduced to the observed value Cob if excess values exist. At other time levels, a formulation corresponding to (4) is applied:

 
formula

After the cloud-cover correction the humidity is adjusted according to (6) in order to be consistent with the cloud formula (2),

 
q′(k, j) = qsat (T, p){RHc + (1 − RHc)[C′ (k, j)]1/2},
(6)

where q′ is the corrected specific humidity and qsat is the saturation specific humidity at the current temperature and pressure. Finally, if convective precipitation is computed in the atmospheric model, the diagnosed cloud cover for tropospheric clouds above the lifting condensation level is corrected according to (7) similar to the formulation by Slingo (1987). The values of the constants a3 and a4 are given in the appendix:

 
Ccv = a3 + ln(1 + a4Q).
(7)

The remaining HIRLAM variables (V10 m, T2 m, q2 m, Q, and ps) explained in Fig. 1 are single level variables depending only on time.

The atmospheric fluxes of sensible and latent heat require the specification of wind, temperature, and humidity at a height of 10 m. The wind components require no interpolation since they are supplied at the computation level. The wind predicted by the atmospheric model is not modified by observations. Predicted temperature and specific humidity are available at 2 m. The initial values of temperature and dewpoint are taken from the observations at the individual road stations. The latest observations are usually no more than about 10 min old. The temperature and humidity at the 10-m level are obtained by linear interpolation between the values at 2 m and at the lowest atmospheric model level, currently at about 32 m.

A time variability is described for temperature and dewpoint in a way similar to (4)

 
y′( j) = y( j) + [yoby(0)] exp[−a5 (ttf)],
(8)

where tf is the initial time of the forecast. If a5 is set to zero, (8) expresses a bias correction equal to the difference between the atmospheric model value and the corresponding observation. For T2 m, the choice of a5 corresponds to an “e-folding” time of decay equal to 10 h. However, for dewpoint the modification is determined by the initial bias correction. This may be justified by a smaller diurnal variability of Td2 m and a less accurate forecast quality than for T2 m.

The precipitation intensity Q is in the form of rain if the lowest model level temperature is above 0°C, otherwise it is snow.Currently Q is not modified from any observations. The same applies to surface pressure that is used to compute air density.

The atmospheric data are linearly time interpolated between hourly values in order to supply a realistic input to the RCM every time step of the numerical integration.

Forecast

The basic problem is to carry out an accurate solution of the local energy budget at the road surface in order to predict the correct road conditions. All relevant processes, that is, radiation, surface sensible and latent heat flux, precipitation, ground heat conduction, and additional heat sources due to other effects—for example, the effects of traffic—need to be addressed in order to make accurate forecasts of the road conditions.

Solar and longwave radiative heat flux

The solar and infrared radiation are computed from the radiation scheme used in the atmospheric HIRLAM model (Sass et al. 1994; Savijärvi 1990; Sass and Christensen 1995). The net solar flux density ϕRs at the ground is computed according to (9) as a linear combination of a clear air part ϕRsa and a cloudy part ϕRsc:

 
ϕRs = (1 − α)[ϕRsa(1 − CM) + ϕRscCM],
(9)

where α is a surface albedo for solar radiation and CM is a total cloud cover determined from a maximum overlap assumption. The clear air term is parameterized according to Savijärvi (1990):

 
formula

In (10), S is the solar constant and p00 = 1000 hPa is a reference pressure. The first term in (10) depending on the zenith angle θ concerns the stratospheric absorption due to ozone. A major contribution to the extinction of solar radiation comes from tropospheric absorption due to water vapor, CO2, and O2. This is parameterized according to the second term in (10). Here, us is the vertically integrated water vapor path (in centimeters) throughout the atmosphere. It is linearly scaled by pressure and divided by cosθ. The last term involving two contributions describes the effect of scattering. The first contribution arises from scattering of the incoming solar beam, while the second one is a compensating effect due to reflected radiation, which is backscattered from the atmosphere above. The coefficients a6 and a7 (see appendix) that are larger than 1 represent a crude inclusion of effects due to aerosol absorption and scattering, respectively.

The cloudy contribution ϕRsc of (9) is given by (11):

 
formula

In (11), ϕRsH is the solar flux density at the top of the uppermost cloud layer. It is given by a formula corresponding to (10), and the surface albedo appearing in the backscattering term of (10) is replaced by an albedo representing the cloudy atmosphere below. Thetransmittance of flux density from top to bottom of the cloudy atmosphere is described by (pH, ps). A simple formula (12) has been fitted to results obtained with more detailed radiation schemes (Slingo et al. 1982; Stephens 1978; Liou and Wittmann 1979):

 
formula

where MT is the vertically integrated cloud water from the top of the uppermost cloud layer to a level below, in this case the bottom of the cloudy atmosphere. In (11) the denominator takes into account multiple reflections between ground and cloud. The constant a8 accounts for absorption in reflected beams. The values of a8a11 are given in the appendix.

Shading effects are taken into account in eight azimuth intervals. For each interval, a threshold shading angle representing the mean height of shading objects has been measured. For solar heights smaller than the threshold value, the flux density from direct solar beam is zero. For the clear-sky part of the radiation, it is assumed that the diffuse flux is a fixed fraction (10%) of the global radiation. This choice is reasonable according to measurements (e.g., Paltridge and Platt 1976, 116–118). The diffuse radiation is reduced by a sky-view factor KΩ representing the integrated effect of flux density originating from different directions within the visible part of the total solid angle 2π. It is assumed that radiation in shaded directions does not contribute to the net radiation. Hence, the net infrared radiation flux density ϕRt is computed according to (13):

 
ϕRt = ɛ0(ϕRta + ϕRtcσT4s) KΩ.
(13)

In (13), σ is the Stefan–Boltzmann constant and ɛ0 = 0.90 is the surface emissivity. The flux contribution ϕRta from the clear atmosphere is determined by

 
formula

In (14), the summation is over all vertical levels of the atmospheric model. Here, B(Tk) is the Planck blackbody radiation associated with the temperature Tk. The emissivity function E (Savijärvi 1990) takes into account the water vapor line spectrum:

 
formula

The emissivity E in (15) is defined from the top of the layer i to the bottom of layer j. The pressures pi−1/2 and pj+1/2 apply to model layer intersections. In (14), b1 represents the effect of radiation from other gases than water vapor, and the last term involving the specific humidity of the lowest atmospheric level N parameterizes the additional flux density due to water vapor continuum. The total downward infrared flux density from the atmosphere is obtained by adding a cloudy contribution ϕRtc defined as radiation transmitted to theground from an effective cloud cover Ce (Sass et al. 1994) according to the empirical expression (16)

 
formula

where B(Te) is the blackbody radiation associated with radiation from the effective cloud cover at temperature Te. The term in brackets is a fractional transmission estimated from the clear air contribution ϕRta and a low tropospheric effective temperature Ta. The values of the constants b1b8 are given in the appendix. The total net radiation ϕR at the road surface isthe sum of the net solar and infrared flux densities.

Sensible and latent heat flux

Traditional drag formulas as given by (17a) and (17b) are used to describe the fluxes of sensible and latent heat:

 
formula

where Z = 10 m, Cp is the specific heat of moist air at constant pressure, ρZ is the air density, and θZ and θs are potential temperatures at the computation level Z and at the surface, respectively. Similarly, qZ and qs are specific humidities at the same levels. The latter is determined by a surface wetness parameter Ws/Wc, where Wc = 0.5 kg m−2 (0 ≤ Ws/Wc ≤ 1), qsat(Ts) is the saturation specific humidity at the surface temperature Ts, and Lvs is the specific latent heat of evaporation if Ts > 0°C, otherwise it is the specific heat of sublimation. Here, k is the von Kármán constant.

In (17a) and (17b), the formulation of the dimensionless drag coefficients is an important issue since the road surface, which is characterized by a small roughness length, is normally embedded in a quite inhomogeneous environment different from a plane road. In recent years, several observational studies have shown that the heat and momentum fluxes over inhomogeneous surfaces are quite complex. For example, Beljaars and Holtslag (1991) show that the local roughness length to be applied for heat and moisture in these conditions may be several orders of magnitude smaller than the momentum roughness length. Currently, the roughness length used in (17a) and (17b) is z0 = 10−4 m, which is a typical value for very flat surfaces (e.g., Garrat 1977). The stability function f (Ri, Z/z0) depending on the Richardson number and Z/z0 is according to Louis (1979) for the unstable boundary layer. In stable stratification the neutral drag coefficients are retained in order to crudely account for the possibility that turbulence at the road surface may be influenced by mechanically generated turbulence in the environments characterized by a larger roughness length.

Ground heat flux

The solution ofthe parabolic equation of heat conduction (18) has been described by Sass (1992). An implicit numerical scheme is used (Smith 1978),

 
formula

The grid spacing is irregular with thin layers close to the surface. The temperature at the bottom layer (level 15) is determined by a climatological estimate depending on the time of year. The values of the heat conductivity λG, density ρG, and specific heat capacity CG are constant (see appendix).

Freezing or melting of precipitation and road water (ice)

The prediction equations for road water Ws (kg m−2) and ice (Is) are given by (19) and (20):

 
formula

In (19), Qrn is the precipitation flux density as rain. The second term represents formation of dew. The third term corresponds to melting if a preliminary updated road surface temperature is above 0°C. It is a sink term if the preliminary surface temperature is below the freezing point. In both cases, the magnitude is determined by the heat convergence in the uppermost road layer, due to ground heat flux from below and atmospheric heat fluxes. The last term in (19) is the runoff of road water. On a smooth road surface, water cannot accumulate without runoff at a rather low threshold value. It is assumed that the runoff ψ equals the sink necessary to prevent the road surface water to increase above the threshold value Wc. Similarly, in (20), the first term is the source term due to precipitation, and the second term a source term due to rime. The third term is the conversion term between ice and water or vice versa, corresponding to the third term of (19) but with opposite sign. In case of freezing or melting the temperatures in the road are recalculated with an upper boundary value Ts(1) = 0°C.

In addition to the already mentioned energy source due to latent heat flux, the energy release from the first and third terms of (19) and (20) are taken into account. This involves the energy release from freezing rain and melting snow.

Additional energy sources

The additional heat sources, such as influence from local topography and the effect of traffic, which gives rise to frictional heating of the road surface and modified turbulence conditions, are parameterized by the term ϕ̂(t) given by (1), mentioned in section 2a.

Experiments

To investigate the properties of the RCM system, it is useful to study the performance during extreme meteorological conditions. Two periods are selected for this purpose. The first one is from 21 to 22 March 1995. At this time of year, the energy input from the sun at a latitude of 55°N may amount to several hundreds of watts per square meter and the diurnal range of temperature may be large. Solar shading effects and a knowledge of the cloud distribution in the sky are sometimes important for a precise knowledge of the energy budget at the road surface. In the second period from 15 to 16 January 1995, it is difficult toobtain a precise knowledge of the cloud cover. The effects of cloud forcing and the surface flux correction term are investigated. Furthermore, a discussion is included on the results from operational runs based on a more extensive data sample. For verification, a direct comparison is made between forecasts and corresponding measurements at the road station sites. The difference in valid time between forecast and observation is normally no larger than 10 min.

Period 1 (21–22 March 1995)

On 21 March 1995 dry air blows across Denmark from the northwest and north, between a high pressure west of Denmark and a low pressure toward the northeast. At the coasts, wind speeds up to gale force are recorded. In the morning hours, the northern part of Jutland is partly cloudy, but otherwise the day is sunny and almost cloud free over Denmark. No precipitation occurs. Toward the evening the wind speeds decrease and there is a wind turning to a westerly direction.

During the morning on 22 March 1995, the cloud cover increases, and relatively mild and humid air is advected across the country. This leads transiently to the formation of rime. Figure 2 shows a synoptic chart of clouds, temperature, dewpoint, and wind at 0600 UTC 22 March 1995. The westerly winds increase during the day, becoming 10–15 m s−1. Light rain is observed (and forecasted) in the rather mild air with temperatures at a height of 2 m mainly between +4° and +7°C.

Fig. 2.

Synoptic chart of clouds, temperature, dewpoint, and wind in Denmark and surrounding areas at 0600 UTC 22 March 1995.

Fig. 2.

Synoptic chart of clouds, temperature, dewpoint, and wind in Denmark and surrounding areas at 0600 UTC 22 March 1995.

Of most relevance to the present study is the occurrence of road surface temperatures below 0°C on both nights and the fact that the advection of humid air gives rise to rime at most locations. Also, the clear skies on 21 March give the opportunity to study the effects of solar shading in rather extreme conditions.

The behavior of the RCM system as an average of 15 forecasts with initial times at 3-h intervals, the first one starting at 0000 UTC 21 March, is illustrated in Figs. 3a–d. These figures apply also to averaging of results for 200 operational road station sites. Figure 3a shows the average skill of frost prediction as measured by the percentage of frost that is predicted (hit rate) and the percentage of forecasted frost that is not observed (false alarm). Figure 3b shows the mean error (dashed line) and the mean absolute error (solid line) of surface temperature. The same measures are displayed in Fig. 3c and Fig. 3d for the predicted temperature and dewpoint at 2 m, respectively. Figures 4a–d show the corresponding results valid for a persistency forecast (see figure caption).

Fig. 3.

(a) Average skill of frost prediction for 200 road station sites during 21–22 March 1995. Solid line applies to the percentage of observed frost that is predicted (hit rate). Dashed line applies to the forecasted frost that is not observed (false alarm). (b) Mean absolute error (solid line) and mean error (dashed line) of the forecasted road surface temperature during the period 21–22 March 1995, as an average for 200 road station sites. (c) As in (b) but valid for temperature at a height of 2 m. (d) As in (b) but valid for dewpoint at a height of 2 m.

Fig. 3.

(a) Average skill of frost prediction for 200 road station sites during 21–22 March 1995. Solid line applies to the percentage of observed frost that is predicted (hit rate). Dashed line applies to the forecasted frost that is not observed (false alarm). (b) Mean absolute error (solid line) and mean error (dashed line) of the forecasted road surface temperature during the period 21–22 March 1995, as an average for 200 road station sites. (c) As in (b) but valid for temperature at a height of 2 m. (d) As in (b) but valid for dewpoint at a height of 2 m.

Fig. 4.

(a) As in Fig. 3a but valid for a persistency forecast, (b) As in Fig. 3b but valid for a persistency forecast. (c) As in Fig. 3c but valid for a persistency forecast. (d) As in Fig. 3d but valid for a persistency forecast.

Fig. 4.

(a) As in Fig. 3a but valid for a persistency forecast, (b) As in Fig. 3b but valid for a persistency forecast. (c) As in Fig. 3c but valid for a persistency forecast. (d) As in Fig. 3d but valid for a persistency forecast.

When comparing Figs. 3a–d with Figs. 4a–d, it is evident that the RCM forecasts perform well compared to persistency. This is most evident for the road surface temperature (Figs. 3b and 4b). The frequency of observed frost that is predicted (Figs. 3a and 4a) stays well above 80% after 5 h, whereas the result of the persistency forecast declines to a value below 40%. The percentage of forecasted frost that is not observed shows somewhat poorer results. It turns out that this is largely due to phase errors during transitions from positive to negative temperatures or vice versa. The mean error (bias) for surface temperature is less than 0.5°C for both persistency and RCM forecasts.

When evaluating the quality of the results it is necessary to take into account the rather extreme conditions during daytime on 21 March. For example, at 1200 UTC the standard deviation of the surfacetemperature computed on the basis of the 200 stations is 4.3°C. This large value occurs as a result of the large effects of shading. The lowest observed value is 2°C at a shaded location (57.1°N, 9.4°E) in the northern part of Jutland. The corresponding 3-h model prediction with initial time of 0900 UTC is 1°C higher. Similarly, the observed maximum surface temperature is about 20°C at a sun-exposed station (55.6°N, 9.7°E) in the eastern part of Jutland. A temperature of 19°C was predicted in the 3-h forecast, starting with a temperature of 10°C at 0900 UTC.

The forecasting quality of surface temperature can also be characterized by the tendency correlation expressing the correlation between forecasted and observed changes. By definition, the tendency correlation is zero for a persistency forecast. However, the correlation is high for the RCM forecasts. As an average over 200 stations, it is 0.90 for 3-h forecasts. It turns out that sun-exposed stations have a very high tendency correlation. An average computed for five such stations gives a correlation of 0.98. A lower value (0.80) is obtained as an average for five densely shaded stations. The smaller correlation for shaded stations may be explained by a generally smaller amplitude of the temperature changes. Small changes are often difficult to predict due to local effects.

The prediction of surface temperature on bridges is often considered to be particularly difficult. However, it turns out that the RCM is capable of producing an average tendency correlation equal to 0.96 for five bridges. This is almost as high as that for (other) sun-exposed stations. Although the computation of the tendency correlation is based on a small data sample, the numbers agree well with results based on data from other periods (see section 4c).

Also, the predictions for temperature and dewpoint at a height of 2 m (Figs. 3c and 3d) are superior to the corresponding results valid for persistency forecasts (Figs. 4c and 4d). For the dewpoint temperature, however, this applies only beyond a forecast range of 2 h. It turns out that the advection of humid air giving rise to dewpoint temperatures above the surface temperature in the morning hours of 22 March is rather well forecasted.

Period 2 (15–16 January 1995)

The morning hours of 15 January 1995 are dominated by cloudy weather and scattered light precipitation; that is, below 1 mm of precipitation is observed for the Danish area. Fog is reported at some places. Temperatures are mainly above the freezing point, but frost occurs at some places in clear air over the northern part of Jutland. A transient high pressure ridge forms during the day. The pressure gradient over Denmark becomes small. The westerly winds decrease and become variable in the evening. The cloud cover is very variable across the country during the daytime. In the late afternoon both temperature and dewpoint at 2 m, as well as the road surface temperature, drop significantly due to a small pressure gradient force and clear sky conditions at many places. This situation leads to the formation of rime and slippery roads in the evening and night, in particular in the eastern part of Denmark. The rime formation is reinforced during the night as winds start to increase from the south and southeast, with the advection of moist air across the country. This is illustrated in Fig. 5 showing synoptic observations of temperature, dewpoint, wind, weather, and clouds at 0000 UTC 16 January 1995 for the island of Zealand and the surroundings. The city of Copenhagen is marked by the letter “C.” In the Western part of the country the advection of cloudy, mild air starts earlier during the evening, and the slippery road conditions are of a shorter duration there. Fairly strongwinds from southerly directions are observed during the daytime of 16 January, at the coasts reaching a speed of 10–15 m s−1. It is cloudy, with temperatures well above the freezing point. However, toward the evening, clear skies appear again, first in the southern part of the country.

Fig. 5.

Synoptic observations of clouds, temperature, dewpoint, and wind at 0000 UTC 16 January 1995 on the island of Zealand and the surrounding smaller islands. The city of Copenhagen is marked by the letter “C.”

Fig. 5.

Synoptic observations of clouds, temperature, dewpoint, and wind at 0000 UTC 16 January 1995 on the island of Zealand and the surrounding smaller islands. The city of Copenhagen is marked by the letter “C.”

Also in this case, 15 subsequent forecasts are carried out for 200 road station sites, as mentioned previously. The accurate forecasting of road surface temperatures is rather difficult for this period due to the large spatial and time variability of clouds, including fog.

We may again compare the forecasting quality based on data from all stations with similar results of persistency forecasts. Figures 6a and 6b display the hit rate (solid lines) and false alarm (dashed lines) of the frost prediction that applies to the RCM forecasts and persistency, respectively. Also, the mean error (dashed lines) and the mean absolute error (solid lines) of the surface temperature forecasts are shown (Figs. 6c and 6d).

Fig. 6.

(a) Average skill of frostprediction for 200 road station sites during 15–16 January 1995. Solid line applies to the percentage of observed frost that is predicted (hit rate). Dashed line applies to the forecasted frost that is not observed (false alarm). (b) As in (a) but valid for a persistency forecast. (c) Mean absolute error (solid line) and mean error (dashed line) of the forecasted road surface temperature during the period 15–16 January 1995, as an average for 200 road station sites. (d) As in (c) but valid for a persistency forecast.

Fig. 6.

(a) Average skill of frostprediction for 200 road station sites during 15–16 January 1995. Solid line applies to the percentage of observed frost that is predicted (hit rate). Dashed line applies to the forecasted frost that is not observed (false alarm). (b) As in (a) but valid for a persistency forecast. (c) Mean absolute error (solid line) and mean error (dashed line) of the forecasted road surface temperature during the period 15–16 January 1995, as an average for 200 road station sites. (d) As in (c) but valid for a persistency forecast.

The mean absolute errors are considerably reduced when using the RCM system, in particular for the last part of the forecast. In addition, this applies to the forecasted 2-m temperature and dewpoint (not shown). The intensified rime formation during the evening and night, as a result of advection of moist air is reasonably well forecasted, with dewpoint temperatures well above the road surface temperature.

The hit rates are lower than the corresponding values obtained in the first period examined (Fig. 3a). This may be attributed to the inaccuracies arising from the difficult cloud prediction. It is possible to identify distinct phase errors. This is illustrated by data from stations in the area of Copenhagen. The time of transition to frost in the afternoon on 15 January is delayed by 1–2 h for three out of four stations. It turns out that an important reason for this phase error is a positive bias in the forecasted 2 m temperature. This is up to +4°C too high in the 3 h forecast valid at 1800 UTC 15 January. The positive bias is mainly due to the way of diagnosing 2-m temperature and humidity in the atmospheric model. In coastal regions, as considered here, the diagnosed values apply to conditions over sea rather than land if the physiographic data of the model indicate that the fraction of a model grid square covered by sea is larger than 50%. The phase delay of the transition to frost is considerably reduced in an experiment (not shown) where the observed temperatures at 2 m are imposed instead of the values strongly influenced by the sea.

The problem of producing accurate surface temperature forecasts is illustrated further by considering the results of mean error and mean absolute error for a single station at the airport of Copenhagen (55.60°N, 12.62°E). Figure 7a applies to a reference experiment including 15 subsequent forecasts as mentioned previously. The mean error is very low, but the mean absolute error exceeds 1°C at a forecast length of 5 h. The effect of using observed temperatures and cloud cover is shown in Fig. 7b. The observed temperatures at 2 m replace the corresponding forecast values, and it is demanded that a linearly time-interpolated observed total cloud cover should agree with the total model cloud cover derived from a maximum cloud overlap assumption. The cloud cover observations are available only every 3 h. It is clear from Fig. 7b that this modified forcing leads to a significant improvement. The mean absolute error does not exceed 0.5°C during the first 4 h. Even in this case, the information about clouds is not perfect due to inaccuracies imposed by the linear time interpolations and a lack of detailed knowledge about the vertical cloud structure.

Fig. 7.

(a) Mean absolute error (solid line) and mean error (dashed line) of the forecasted road surface temperature during the period 15–16 January 1995 for a road station at the airport of Copenhagen. (b) As in (a) but using prescribed cloud cover and 2-m temperature during the forecasts.

Fig. 7.

(a) Mean absolute error (solid line) and mean error (dashed line) of the forecasted road surface temperature during the period 15–16 January 1995 for a road station at the airport of Copenhagen. (b) As in (a) but using prescribed cloud cover and 2-m temperature during the forecasts.

It is furthermorerelevant to investigate the impact of the surface flux correction term (1). This is done by carrying out an experiment including all road station sites, with no surface flux correction to compensate for lacks such as insufficient cloud cover information and model deficiencies. The results are given in Fig. 8a, which shows a positive bias and a significant increase of absolute error compared to the corresponding reference experiment. A positive impact of the surface flux correction can also be demonstrated if the extreme condition is imposed that the cloud cover is forced to be zero. Figure 8b shows the results obtained when the surface flux correction is allowed to operate, while Fig. 8c shows the corresponding results with no flux correction included. It is clear from these figures that both bias and mean absolute error are substantially reduced when the flux correction term is allowed to operate.

Fig. 8.

(a) As in Fig. 6c but valid for an experiment with no surface flux correction term. (b) As in Fig. 6c but valid for an experiment with no clouds. (c) As in Fig. 6c but valid for an experiment with no clouds and no surface flux correction term.

Fig. 8.

(a) As in Fig. 6c but valid for an experiment with no surface flux correction term. (b) As in Fig. 6c but valid for an experiment with no clouds. (c) As in Fig. 6c but valid for an experiment with no clouds and no surface flux correction term.

In addition, a degradation of forecasting quality (relative to a reference run) occurs in an experiment where cloud cover observations are discarded (not shown). This confirms the expectation that observed cloud cover is beneficial to the forecasting skill.

Operational results

The verification results presented above are based on data from rather short periods of only two days. It is therefore relevant to investigate whether a similar forecasting skill can be achieved on a larger data sample. The operational production of data allows for this type of comparison. In general, it is found that the forecasting skill of operational runs verified on a monthly basis is no less than that presented in the case studies mentioned above. This is briefly illustrated in Figs. 9a and 9b showing verification results of operational surface temperature forecasts valid for February and March 1996, respectively. The mean absolute error of the RCM forecasts is shown by solid lines. The corresponding results valid for persistency forecasts are given by dashed lines. The dot–dashed curves apply to linear trend forecasts defined by linear time extrapolation of a mean observed surface temperature tendency during the last hour prior to the forecast. The verification is based on the averaging of results from all stations. New forecasts are made every hour. This is a disadvantage with respect to cloud cover observations, which are currently available only every 3 h.

Fig. 9.

(a) Operational monthly verification results for February 1996 as an average for 200 road station sites. The mean absolute error of the RCM forecasts is shown by the solid line. The dashed line and the dot–dashed line apply to persistency forecasts and linear trend forecasts, respectively (see text). (b) As in (a) but valid for operational forecasts in March 1996.

Fig. 9.

(a) Operational monthly verification results for February 1996 as an average for 200 road station sites. The mean absolute error of the RCM forecasts is shown by the solid line. The dashed line and the dot–dashed line apply to persistency forecasts and linear trend forecasts, respectively (see text). (b) As in (a) but valid for operational forecasts in March 1996.

The results shown for both months are consistent with the findings of the case studies. It is seen that the mean absolute error of persistency forecasts is significantly greater than the corresponding RCM value, except for the first hour where the different alternatives provide a very similar mean absolute error. Beyond a forecast length of 3 h the linear trend forecasts are poor compared to the persistency forecasts.

The operational results may be illustrated further by considering features such as road temperature forecast errors in specified error intervals, and the error as a function of different times of the day, that is, a diurnal error cycle. Figure 10a shows the forecast errors during February 1996 in six error intervals, as an average for all 200 station sites. The two extreme intervals cover the largest negative and positive errors, with absolute values larger than 2°C. The forecast range is 3 h. It is seen that almost 80% of the forecasts are within 1°C of the observed values. There is a slight predominance of small positive deviations over negative values. The remaining errors between 1° and 2°C and above 2°C, respectively, are rather symmetrically distributed around zero. The absolute errors above 2°C occur in approximately 5% of all forecasts. The error frequencies for both shaded station sites and sun-exposed locationsarequite similar to the average results valid for all stations. An error frequency diagram presenting the results from five shaded stations are shown in Fig. 10b.

Fig. 10.

(a) Error frequencies (%) of road surface temperature forecasts in February 1996 at 3 h, as an average for 200 station sites. The frequencies are in 1°C intervals. The two extreme intervals represent errors with absolute values larger than 2°C. (b) As in Fig. 10a but valid for five shaded station sites.

Fig. 10.

(a) Error frequencies (%) of road surface temperature forecasts in February 1996 at 3 h, as an average for 200 station sites. The frequencies are in 1°C intervals. The two extreme intervals represent errors with absolute values larger than 2°C. (b) As in Fig. 10a but valid for five shaded station sites.

The diurnal error cycle is illustrated in Figs. 11a,b. The first figure (11a) shows the 3-h surface temperature forecast errors valid at the hours displayed on the abscissa (all stations included). The bias is shown by the dashed line and the mean absolute error by the solid line. It is seen that the errors are small at night. During the daytime where the surface energy fluxes are larger due to solar radiation, the errors are bigger but still not much above 1°C in mean absolute error. The corresponding bias is below +0.5°C. The typical result for a sun-exposed station is quite similar to that of Fig. 11a, but the positive bias at 0900 UTC and 1200 UTC tends to be slightly larger. On the other hand, for a shaded location, the bias at 1200 UTC tends to be negative by up to about −0.5°C. Figure 11b applies to the same five shaded locations as in Fig. 10b.

Fig. 11.

(a) Mean error (dashed line) and mean absolute error (solid line) of forecast road surface temperature in February 1996 on different times of the day. Results apply to +3-h forecasts evaluated at the hours shown on the abscissa (average results for 200 stations). (b) As in (a) but valid for five shaded station sites (same stations as in Fig. 10b).

Fig. 11.

(a) Mean error (dashed line) and mean absolute error (solid line) of forecast road surface temperature in February 1996 on different times of the day. Results apply to +3-h forecasts evaluated at the hours shown on the abscissa (average results for 200 stations). (b) As in (a) but valid for five shaded station sites (same stations as in Fig. 10b).

Concluding remarks

Predictive skill has been demonstrated with the numerical RCM system as presented in sections 2 and 3. For example, this is illustrated by the comparisons of the RCM forecasted surface temperature with the corresponding results of persistency forecasts. A similar remark applies to temperature and dewpoint at 2 m. The verifications presented here are consistent with operational results, which tend to have slightly smaller errors in verifications on a monthly basis due to averaging over many weather situations. The slippery road conditions are in both case studies considered, dominated by rime deposition that occurs rather frequently in Denmark. In other weather situations, accurate precipitation forecasts—for example, snow or freezing rain, is a prerequisite for the prediction of slippery roads.

The results of the first period (21–22 March 1995) show that the RCM-system can describe the large observed differences between shaded and sun-exposed locations. The importance of accurate atmospheric input—for example, 2-m temperature and cloud cover—is illustrated by the sensitivity experiments of the second period (15–16 January 1995). It has been shown, both for cloudy and cloud-free conditions, that the flux correction term of the RCM system has the potential to improve the accuracy of the surface temperature forecasts.

The possible improvements of the results produced by the RCM system depend first on the additional quality that may be obtained in the uncorrected input data from the atmospheric model and, second, on the meteorologically significant refinements that can be put into the road conditions model itself.

The problem of specifying realistic temperature and humidity at 2 m can be solved if the atmospheric model is run with a higher resolution in the horizontal. This has already been done on a test basis (Nielsen et al. 1995). As the horizontal resolution becomes sufficient, the details of the associated land–sea mask will tend to give more realistic temperature and dewpoint even for road station sites close to the coast. If the horizontal resolution of the atmospheric model is insufficient to guarantee that a diagnosis valid for land is made for all road station sites, interpolated values should be computed from the neighboring grid points with a land diagnosis. The time variability of the temperature and dewpoint corrections (8) needs tuning to be optimal for a given atmospheric model.

A more crucial question concerns the possibility to predict clouds and precipitationaccurately. As indicated previously, the radiation scheme of the RCM should ideally include cloud geometry effects in order to describe properly the energy input to the road surface. This ultimately requires a cloud predictability in the atmospheric model on scales below 10 km. Currently, considerable efforts are being invested at various meteorological forecasting centers, to improve the forecast quality of high resolution operational atmospheric models (e.g., Ballard et al. 1993; Nielsen et al. 1995; Schrodin 1996). However, it is still an open question whether a substantial improvement of the cloud and precipitation prediction on scales below 10 km can be achieved within the near future. Meanwhile, a cloud correction procedure similar to that expressed by (4)–(6) can be used. It may be possible to improve the formulation by making the coefficient a2 in (4) and (5) variable with respect to the weather situation. An improvement can be expected to the extent that the lifetime of local cloud phenomena depend on the type of weather situation. A revision of this kind needs statistical verification on an extensive sample of cloud data. An increased quality of the cloud cover analyzed at the initial time of the forecast—for example, by using satellite data—is likely to have a positive impact as indicated by the result of the present study, that observed clouds have a positive impact on the forecast accuracy. Other data may be useful. For example, the present data assimilation used in the HIRLAM model does not make use of data from weather radar. The assimilation of such data in the data correction procedure of the RCM is therefore likely to improve the computation of precipitation.

The RCM system may be improved in other aspects; for example, there is still some uncertainty about the formulation of turbulent fluxes at the road surface, as mentioned in section 3. In addition, various physical constants such as the heat conductivity, the heat capacity, and the radiative properties of the road, can be made station dependent to be more realistic. However, the lacks are partly compensated by the flux correction term.

It should be borne in mind that a verification problem arises if an accuracy better than 0.5°C is demanded. This is approximately the level when measurement techniques such as the depth of measurement of the road surface temperature start to become important—for example, when a displacement of less than 1 cm is of importance. Also, the sampling frequency of the measured data is of significance. Even a data frequency as high as 10 min is not sufficient in all cases since the temperature may change by more than 0.5°C in 10 min. Hence, it is a challenge to achieve an operational accuracy better than that presented in Fig. 7b obtained with hindcasted clouds, which are not available in an operational context.

Acknowledgments

The author wishes to thank his colleague Jess U. Jørgensen for valuable technical assistance related to the operational environment at the Danish Meteorological Institute.

REFERENCES

REFERENCES
Ballard, S. P., R. Robinson, R. T. H. Barnes, S. Jackson, and S. Woltering, 1993: Development and performance of the new meso-scale model. Forecasting Research Tech. Rep. 40, U.K. Meteorological Office, 29 pp. [Available from U.K. Meteorological Office, London Road, Bracknell, Berkshire RG 12 2SZ, United Kingdom.]
.
Beljaars, A. C. M., and A. A. M. Holtslag, 1991: On flux parameterization over land surfaces foratmospheric models. J. Appl. Meteor., 30, 327–341.
Bogren, J., and T. Gustafsson, 1991: Nocturnal air and road surface temperature variations in complex terrain. Int. J. Climatol., 11, 443–455.
——, and ——, 1994: A combined statistical and energy balance model for prediction of road surface temperature. Proc. 7th Int. Road Weather Conf., Seefeld, Austria, Central Institute for Meteorology and Geodynamics, 43–51.
Garrat, J. R., 1977: Review of drag coefficients over oceans and continents. Mon. Wea. Rev., 105, 915–929.
Källen, E., 1996: HIRLAM Documentation Manual, System 2.5. Swedish Meteorological and Hydrological Institute, 178 pp. [Available from SMHI, S-60176 Norrköping, Sweden.]
.
Liou, K. N., and G. D. Wittmann, 1979: Parameterization of the radiative properties of clouds. J. Atmos. Sci., 36, 1261–1273.
Louis, J. -F., 1979: A parametric model of vertical eddy fluxes in the atmosphere. Bound.-Layer Meteor., 17, 187–202.
Nielsen, N. W., B. H. Sass, and J. Jørgensen, 1995: Meso-scale forecasts with an atmospheric Limited Area Model. Meteor. Appl., 2, 351–361.
Paltridge, P. W., and C. M. R. Platt, 1976: Radiative Processes in Meteorology and Climatology. Elsevier Scientific Publishing, 318 pp.
Rayer, P. J., 1987: The meteorological office forecast road surface temperature model. Meteor. Mag., 116, 180–191.
Sass, B. H., 1992: A numerical model for prediction of road temperature and ice. J. Appl. Meteor., 31, 1499–1506.
——, 1994: The DMI Operational HIRLAM Forecasting System Version 2.3—A Short Summary. DMI Tech. Rep. 94-8, 10 pp. [Available from the Danish Meteorological Institute, Zyngbyvej 100, 2100 Copenhagen, Denmark.]
.
——, and J. H. Christensen, 1995: A simple framework for testing the quality of atmospheric limited-area models. Mon. Wea. Rev., 123, 444–459.
——, L. Rontu, and P. Räisanen, 1994: HIRLAM-2 Radiation Scheme: Documentation and tests. HIRLAM Tech. Rep. 15, 42 pp. [Available from SMHI, S-60176 Norrköping, Sweden.]
.
Savijärvi, H., 1990: Fast radiation parameterization schemes for mesoscale and shortrange forecast models. J. Appl. Meteor., 29, 437–447.
Schrodin R., Ed., 1996: Quarterly Report of the Operational NWP-Models of the Deutscher Wetterdienst. Deutscher Wetterdienst Rep., 59 pp. [Available from DWD, Frankfurter Str. 135, D-63004, Offenbach a. M., Germany.]
.
Slingo, A., S. Nicholls, and J. Schmetz, 1982: Aircraft observations of marine stratocumulus during JASIN. Quart. J. Roy. Meteor. Soc., 108, 833–856.
Slingo, J. M., 1987: The development and verification of a cloud prediction scheme for the ECMWF model. Quart. J. Roy. Meteor. Soc., 113, 899–928.
Smith, G. D., 1978: Numerical Solutions of Partial Differential Equations: Finite Difference Methods. Oxford University Press, 304 pp.
Stephens, G. L., 1978: Radiative properties of extended water clouds. Part I. J. Atmos. Sci., 35, 2111–2122.
Thornes, J., 1996: The cost-benefit of winter road maintenance in the United Kingdom.Proc. Eighth Int. Road Weather Conf., Birmingham, United Kingdom, University of Birmingham, 1–10.

APPENDIX

Values of Constants

 
formula

Footnotes

Corresponding author address: Mr. Bent H. Sass, The Danish Meteorological Institute, Lyngbyvej 100, 2100 Copenhagen, Denmark.