- Email: [email protected]
Ensemble forecasts of road surface temperatures
Atmospheric Research 187 (2017) 33–41
Contents lists available at ScienceDirect
Atmospheric Research journal homepage: www.elsevier.com/locate/atmosres
Ensemble forecasts of road surface temperatures Zbyněk Sokol a,⁎, Vojtěch Bližňák a, Pavel Sedlák a, Petr Zacharov a, Petr Pešice a, Miroslav Škuthan b a b
Institute of Atmospheric Physics, Czech Academy of Sciences, Boční II 1401, 141 31, Prague, Czech Republic Czech Hydrometeorological Institute, Na Šabatce 2050/17, 143 06, Prague, Czech Republic
a r t i c l e
i n f o
Article history: Received 15 March 2016 Received in revised form 7 December 2016 Accepted 16 December 2016 Available online 22 December 2016 Keywords: Road weather forecast Road surface temperature Ensemble prediction
a b s t r a c t This paper describes a new ensemble technique for road surface temperature (RST) forecasting using an energy balance and heat conduction model. Compared to currently used deterministic forecasts, the proposed technique allows the estimation of forecast uncertainty and probabilistic forecasts. The ensemble technique is applied to the METRo-CZ model and stems from error covariance analyses of the forecasted air temperature and humidity 2 m above the ground, wind speed at 10 m and total cloud cover N in octas by the numerical weather prediction (NWP) model. N is used to estimate the shortwave and longwave radiation fluxes. These variables are used to calculate the boundary conditions in the METRo-CZ model. We found that the variable N is crucial for generating the ensembles. Nevertheless, the ensemble spread is too small and underestimates the uncertainty in the RST forecast. One of the reasons is not considering errors in the rain and snow forecast by the NWP model when generating ensembles. Technical issues, such as incorrect sky view factors and the current state of road surface conditions also contribute to errors. Although the ensemble technique underestimates the uncertainty in the RST forecasts, it provides additional information to road authorities who provide winter road maintenance. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Forecasts of road surface conditions are highly important for road maintenance services. These forecasts usually apply models based on solving energy balance and heat conduction equations (e.g., Bouilloud et al., 2009; Crevier and Delage, 2001; Greenfield and Takle, 2006; Jansson et al., 2006; Kangas et al., 2015; Rutz and Gibson, 2013; Shao and Lister, 1996). The energy balance equation calculates boundary conditions for the heat conduction equations, the accuracy of which has crucial impacts on forecasts of road surface temperature (RST) and road surface conditions. The boundary conditions depend on the current and future states of the atmosphere, primarily forecasts of shortwave radiation flux, longwave radiation flux, air temperature and humidity, wind speed and amount and type (rain, snow) of precipitation, which are required to calculate energy fluxes between the road surface and the near surface atmosphere. These prognostic variables are obtained from numerical weather prediction (NWP) model forecasts. The boundary conditions are influenced by various errors that are mainly caused by inaccurate NWP model forecasts and parameterizations applied to calculate surface fluxes from the NWP model variables. An improper consideration of shade effects may also contribute to errors. Errors also arise from the energy balance model assumption that
⁎ Corresponding author. E-mail address: [email protected] (Z. Sokol).
http://dx.doi.org/10.1016/j.atmosres.2016.12.010 0169-8095/© 2016 Elsevier B.V. All rights reserved.
the entire surface of the Earth is composed of a horizontally homogeneous road, which is not met in reality. The model that describes heat conduction in the road body also contributes to the overall error, particularly because the characteristics of the road layers (their thicknesses, thermal conductivities and capacities) may not be accurate. However, the contribution of these errors is not very important, especially for shorter lead times, because of the slow response of the surface temperature to the subsurface characteristics. Because of these uncertainties, deterministic forecasts do not provide sufficient information for road maintenance services to make optimal decisions. It is important to have additional information regarding the accuracy of the deterministic forecast in various forms, for example, the probability that the forecasted value will be in a given interval of values. The problem of the uncertainty of deterministic forecasts is also confronted in other meteorological forecasts. To solve this problem, an ensemble technique was introduced, developed and applied to weather forecasting over the last two decades. The ensemble technique consists of generating a set of deterministic forecasts that represent possible future states, and the problem of limited predictability and uncertainty can be addressed (e.g., Leutbecher and Palmer, 2008) by statistical processing of a set of obtained forecasts. The results of the ensemble technique depend fundamentally on how individual members of the ensemble are generated. In meteorological applications, ensembles are usually generated by varying the initial and/or boundary conditions; however, changes in model parameterizations and/or model parameter values have also been frequently used (e.g., Bouallegue et al., 2013).
34
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
In this paper, an ensemble technique for RST forecasting is presented using the METRo-CZ model (Sokol et al., 2014). The ensemble technique is used to develop a probabilistic forecast that the RST will be below a given threshold (0 °C), which is important information for road maintenance services. We applied a method for generating ensembles that simultaneously modifies forecasted air temperature and humidity 2 m above the ground (T2 m, H2 m), wind speed at 10 m (W10 m) and the total cloud cover N in octas obtained from the NWP model. N is used to estimate the shortwave and longwave radiation fluxes. To our knowledge, this is the first application of this or a similar approach. It is noteworthy that the NWP model used in this study (Section 2) does not provide an ensemble forecast. Therefore, an ensemble technique based on the deterministic forecast was applied. This paper is divided into four sections. Section 2 describes the METRo-CZ model, including the data that are used to prepare the initial and boundary conditions. The section also contains a description of the ensemble technique. An evaluation and discussion of the ensemble technique are presented in Section 3, and the conclusions are presented in Section 4. Details concerning the ensemble technique are provided in the Appendices.
2. The METRo-CZ model 2.1. Model description The Institute of Atmospheric Physics CAS in Prague and the Czech Hydrometeorological Institute (CHMI) developed and currently run the METRo-CZ model for forecasts of RST and road surface conditions. The model is based on solving the energy balance and heat conduction equations. It originated as the Model of the Environment and Temperature of Roads (METRo), which was originally developed by Environment Canada (Crevier and Delage, 2001). The METRo model documentation and source code are freely available at http:// home.gna.org/metro/. METRo-CZ is a 1D model that calculates independent forecasts at single points. It uses measured data from road weather stations (RWS) and forecasted data obtained by integration of the ALADIN NWP model as inputs to prepare initial and boundary conditions (Sokol et al., 2014). A detailed description of the data employed in this study is provided in Section 2.2. METRo-CZ and the original METRo comprise 3 basic steps (Crevier and Delage, 2001; Linden and Petty, 2008; Sokol et al., 2014). In the 1st step, which is called initialization, the vertical profile of the ground temperature is calculated. In the 2nd step, called coupling, radiative fluxes forecasted by the NWP model are corrected by an additive correction coefficient. The aim is to ensure that forecasted and measured RSTs are close at the beginning of the METRo-CZ forecast (TMET). The coefficient value is calculated iteratively. In each iteration, the correction coefficient is increased or decreased depending on the sign of the difference between the forecasted and measured RST and the magnitude of the change decreases in iteration steps. The process is stopped when the difference between the RSTs is less than 0.05 °C. The length of the coupling period depends on the beginning of the METRo-CZ forecast and is at least 4 h. If the forecast starts at TMET = 4, 5, …, 9 UTC, then the coupling begins at TNWP = 00 UTC using data from the ALADIN run starting at TNWP. The initialization, which only uses measured data from the RWSs, is run from TNWP-6 h to TNWP. Similarly, if the forecast starts at 10, …, 15 UTC, then the coupling begins at TNWP = 06 UTC, etc. The current version of METRo-CZ differs from METRo in various aspects. The original code was completely rewritten in Fortran 90, and a large majority of the model parameter values are now model inputs that can be easily changed. The data pre-processing was modified and tailored to the data used in the CHMI. The physical core of the model equations is the same as in METRo, and the basic differences between the METRo-CZ and the original METRo model will now be discussed.
METRo-CZ uses N to calculate the shortwave and longwave radiation fluxes. We found that the radiation fluxes calculated directly by the ALADIN NWP model yielded poorer RST forecasts than when N was utilized. We also found that the relationship between N and the shortwave radiation fluxes (Table 2 in (Crevier and Delage, 2001)) used in METRo overestimates global radiation in comparison with measurements at the Czech synoptic stations. Therefore, using measured global radiation and observed N we derived and applied a different relationship (Table 1). METRo-CZ focuses on nowcasting (i.e., very short-range forecasts) with lead times up to 6 h. It uses NWP data with time steps of 1 h and all available data from RWSs and calculates new forecasts every hour.
2.2. Model input data As previously stated, METRo-CZ uses prognostic data from the ALADIN NWP. Specifically, it employs the ALADIN-CZ version of the model, which is operated by the CHMI. ALADIN-CZ is a hydrostatic model that is integrated at a horizontal resolution of approximately 5 km; the integration begins four times per day (00, 06, 12 and 18 UTC), and its lead time is 54 h. The following forecasted data are required: air temperature T2 m and humidity H2 m at 2 m, wind speed W10 m at 10 m, model pressure at the model orography, cloud cover N, and accumulated precipitation and type (rain/snow). These model data are available 4 h after the start of NWP model run. METRo-CZ uses the following measurements from the RWSs: air temperature T2 m and humidity H2 m at 2 m, wind speed W10 m at 10 m, the road weather conditions (a code indicating dry or wet surface), RST, and temperature at 30 cm below the road surface (T30). If any T2 m, H2 m or W10 m data are not available, their values are estimated by linear interpolation of ALADIN-CZ data over time. If information on the road weather conditions is missing, a dry surface is assumed. If RST is not available the measurement term is not used. In this study, we used 27 RWSs, which were selected to cover various parts of the western and central Czech road network (highways and 1st class main roads) under various road conditions (lowlands and mountains). Their positions (marked with crosses) with RWS codes and a map of the highway network over the orography (blue lines indicate open highways) in the western and central Czech Republic are shown in Fig. 1. The RWSs performed measurements every 20 min. The data from each station were available for more than 96% of the measurement period, and the quality check of T2 m, H2 m, W10 m and RST removed less than 1% of the measurements. The quality check compared neighbouring measurements in the time sequences and removed unrealistically large differences in the measurements. Thus, less than 1% of the measurements were removed. The measured road weather conditions were not assessed. These data were available at each RWS for at least 94% of the measurement period, with the exceptions of RWSs J005 and P019, which had measurements for only 46% and 40% of cases, respectively.
Table 1 Original and new coefficients for the shortwave radiation fluxes (see Table 2 in [Crevier and Delage, 2001]). Cloud cover (octas)
D1 (original)
D1 (new)
0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 8/8
1.00 0.97 0.94 0.89 0.85 0.80 0.71 0.65 0.33
0.94 0.91 0.90 0.87 0.83 0.75 0.67 0.48 0.25
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
35
Fig. 1. The highway network (thick blue lines indicate open highways) in the western and central Czech Republic. The positions of the road weather stations used in this study are denoted by crosses with RWS codes. The orography is shown in the background. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
2.3. Ensemble generation 2.3.1. General information on the ensembles The purpose of the ensemble technique is to statistically describe the uncertainty in the deterministic forecasts and possibly refine the deterministic predictions using the obtained statistical characteristics. Forecast error is generally dependent on model errors (caused, for example, by the application of simplified equations and inaccurate numerical solutions) and errors in the input model data (initial and boundary conditions). In our experience, errors in forecasted atmospheric data that represent boundary conditions in METRo-CZ and describe the state of the near-surface layer of the atmosphere significantly influence RST forecasts (Sokol et al., 2014). Errors related to temperature, humidity, wind speed and, especially, radiation fluxes and the occurrence of precipitation (rain or snow) dominate model errors caused by the use of simplified equations. Therefore, we focused on creating ensembles using atmospheric temperature, humidity, wind speed and total cloud cover, and neglected errors and inaccuracies produced by the METRo-CZ equations. We also did not consider errors in the precipitation forecasts because their statistical characteristics are complicated and it is not clear how to use them in the context of the applied ensemble technique. In our calculations, we neglected the shading effect in spite of its importance due to lacking the required data at our disposal at present. We calculated an ensemble of forecasts by generating an ensemble set of boundary data and applied METRo-CZ to this data set. The set of
boundary data was prepared as follows: For forecasts starting at TMET and with projection time δ (in our case, δ = 6 h), we had available hourly NWP forecasts of T2 m, H2 m, W10 m and N, where the lead time t1 corresponded to TMET, and t2 corresponded to TMET + δ. To simplify the notation, we will introduce a vector T of length 4δ containing hourly forecasts for δ lead times and 4 variables: T ¼ ðT2
m ðt1 Þ; H2 m ðt1 Þ; W10 m ðt1 Þ; Nðt1 Þ; …; T2 m ðt2 Þ; H2 m ðt2 Þ; W10 m ðt2 Þ; Nðt2 ÞÞ:
ð1Þ We used an available set of historical NWP forecasts and observations, which we denote Ti,j and Oi,j, where i = 1, …,n refer to an individual forecast or observation for single hours, and j = 1, …,4δ (we repeat that δ = 6 h in our case) refer to components of the vector T. Using the historical data, we calculated the covariance matrix C: C ¼ ck;m ; k; m ¼ 1; …; 4δ;
ð2Þ
where ck;m ¼
n 1X T −Oi;k T i;m −Oi;m n i¼1 i;k
ð3Þ
and the indices of T and O stand for an individual forecast or observation and vector component, respectively. The matrix C was used to describe the error structure of the boundary conditions, i.e., of vector T. To obtain
36
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
the ensemble of forecasts, we generated an ensemble (set) of vectors T whose components were selected such that they fulfilled the covariance relationships described by Eq. (2). The procedure for generating these vectors is described in Appendix A. The covariance matrix C is the core of the applied method. Because METRo-CZ is a 1D model calculated separately for each RWS, it is natural to prepare matrices C for each station individually. Temperature, humidity and wind speed are measured at all RWSs; thus, the observational data in Eq. (3) are available. Because N is, however, not observed at those stations, we instead used satellite data from the Meteosat Second Generation and an algorithm based on the SAF algorithm PGE01 Cloud Mask to estimate N. The algorithm is described in Appendix B. Fig. 2 shows an example of the generated ensembles of T2 m, H2 m, W10 m and N. The ensemble had 50 members and was prepared for station P002 for 15 November 2012, 06 UTC. It should be mentioned that the humidity in (1–3) is expressed as a dew point temperature and then transformed into specific humidity, as required by METRo-CZ, using the perturbed temperature and forecasted pressure. Fig. 2 clearly demonstrates the sizes of the uncertainties in the forecasts of T2 m, H2 m, W10 m and N. N (subpanel d) is of particular concern; it seems to behave chaotically, but this is only partly true when considering the grey area, where 50% of the values are located.
2.3.2. Specific application of the ensemble technique and its limitations We prepared three versions of the ensemble application. The first version perturbed only T2 m, the second perturbed T2 m, H2 m and W10 m, and the third perturbed all the variables (T2 m, H2 m, W10 m and N). The construction of the matrix C differed among the versions. Because the first version always yielded the poorest scores in terms of Relative Operating Characteristics, Brier score and verification rank histogram (Wilks, 2006), which are used to evaluate ensembles in Section 3, it is not further evaluated in this paper. The second and third versions differed by the inclusion of N, which is an important factor in the uncertainty of the METRo-CZ prediction. However, the determination of N may be affected by significant errors that influence the covariance matrix C and the resulting ensembles, and thus may deteriorate ensemble attributes. Therefore, both versions with and without N were tested. In this paper, we used the covariance matrix C calculated using data from the winter season of 2012–2013 (October to March). Although this time period included the evaluation period, it had a small impact on the results, which we verified by gradually excluding just one month (from November to February) of data from the period used to calculate C and compared the METRo-CZ forecasts with those obtained using all data. The overall results were quite similar (not shown here). Therefore, we assume that the calculated covariances were sufficiently robust and the overlap of the time periods, caused by data availability, did not affect significantly the obtained results. In our subjective opinion, one winter season, typically the latest one, is the most suitable length of data for calculating C. The reason is that NWP models are frequently changed, and major changes, which may significantly influence the covariances, are usually performed once per year. As follows from Eq. (1), the values in matrix C depend on the difference between the start of the NWP forecasts and the METRo-CZ models (index j1). This dependence should also include the time of day because the NWP model forecast errors have diurnal variations. Therefore, we calculated C for each hour in a day. The basic shortcoming of the described ensemble technique is that the covariance error structure implicitly assumes that the errors in the T components are normally distributed. This assumption is reasonable for T2 m and H2 m but is not valid for W10 m and N. The violation of normality is not important for W10 m because the W10 m values have a small impact on the METRo-CZ forecasts. N, however, has a significant impact on the METRo-CZ forecasts. The distribution of the N errors is not normal but is approximately symmetric, with a significant peak near zero (Fig. B1); thus, we can expect that the violation of normality may not cause significant errors in the METRo-CZ forecasts. 3. Forecast evaluation 3.1. Implementation of METRo-CZ The RST and road weather conditions were calculated using METRoCZ for the period from 1 November 2012 to 28 February 2013. The forecasts were calculated at the 27 RWS introduced in Section 2.2. The forecasts started every hour and were projected out 6 h. For each hour, we performed an ensemble of 10, 30, 100, 250 and 500 forecasts. The first forecast used the T2 m, H2 m, W10 m and N data forecasted by ALADINCZ, and the remaining forecasts used the perturbed data. 3.2. Evaluation of the ensemble forecasts
Fig. 2. An example of the generated ensembles of T2 m (TA), H2 m (Q), W10 m (W) and N. The ensemble had 50 members and was prepared for station P002 for 15 November 2012, 06 UTC. The filled area in the subpanel d) contains 50% of the N values.
The evaluation was carried out via two closely related directions. First, we evaluated the properties of two types of ensembles generated using perturbations in (i) T2 m, H2 m, W10 m and N (ENS1) and (ii) T2 m, H2 m and W10 m (ENS2). The second direction focused on a discussion of the practical usefulness of the ensembles to express uncertainties in the predictions and forecasts in comparison with the deterministic forecast. Ideally, the true value should only be one additional member of the ensemble and indistinguishable from the members. This attribute of the
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
ensemble can be evaluated using the verification rank histogram (VRH; Wilks, 2006). The VRH for a scalar variable is calculated as follows: for each forecast, the ensemble values are ordered from the smallest to the largest value and the order of the observed values in the ensemble is determined. If the observed value is lower than the smallest ensemble value, the order is 1. If the observed value is between the nth and the (n + 1)th ensemble values, then the order is n + 1, etc. The VRH histogram shows frequencies of the calculated orders of observed values. The VRH of the ideal ensemble should have an approximately uniform distribution. In our plots, we used 10 classes in the histograms. The 1st and 10th classes contained the relative number of observed data values lower than the minimum ensemble values and greater than the maximum ensemble values, respectively. Classes 2 to 9 contained the relative numbers of observed data values in 8 equally wide intervals between ensemble minimum and maximum. Fig. 3 shows an example of the VRH of ensembles ENS1 and ENS2 for K013 RWS and for lead times of 60, 180 and 300 min for 300 ensemble members. The other RWSs had similar VRHs. The high rightmost and especially high leftmost columns indicate an under dispersion of the ensembles, meaning that the ensembles are overconfident and do not cover the entire variability in the observations. The high relative number of measurements that are lower than the minimum of the ensembles can be explained mainly by two factors that were not considered when generating the ensembles. First, rain and/or snow often lead to cooler surface temperatures than does a clean dry road surface, and second, the sky view factor, meaning a parameter that corresponds to the portion of the visible sky limited by relief (Zakšek et al., 2011), is not properly taken into account, which may result in shadowing and unexpected significant temperature drops during the day. The first factor is also related to the fact that the state of the road surface is frequently not known. The comparison of the VRHs for ENS1 and ENS2 clearly indicates that the cloud cover N is an important source of uncertainty in the forecasts. N mainly affected
37
the first class on the histogram, resulting in a higher RST forecast of METRo-CZ compared to the RWS measurements. For given probabilities P = 50%, 70% and 95%, RST intervals were calculated using the ensemble forecasts. For example, for P = 50%, the interval was determined as between the 25% and 75% percentiles of the forecasts for a given term. Fig. 4 shows the results obtained for all the stations as a function of the lead times for ENS1 and ENS2 with various ensemble sizes. Consistent with Fig. 3, it is clear that ensemble ENS1 underestimated the spread. The significantly poorer results for ENS2 confirm the important role of N in the construction of ensembles. From the viewpoint of ensemble sizes, it follows from Fig. 4 that the size should be at least 100 members for higher probabilities, e.g., P = 95%. The obtained results can be potentially utilized to extend the forecasted intervals by applying empirical corrections to better match the measurements. The Model Output Statistics method (MOS) e.g., (Sokol, 2003) is suitable for that task. To predict road surface conditions, it is important to forecast whether the surface temperature will be lower or higher than 0 °C; this forecast may be categorical or probabilistic, the latter of which is preferable. These types of forecasts are frequently evaluated using Relative Operating Characteristics (ROC) and derived characteristics (Wilks, 2006) that include the area under the ROC curve (Area) or by the Brier score (BS) (Wilks, 2006), which is similar to the standard mean-squareerror calculated using probabilities or binary values. For a perfect forecast, Area = 1; for a random forecast, Area = 0.5. The Area value shows the potential usefulness of the analysed forecast and allows a comparison of the deterministic and ensemble forecasts. Fig. 5 shows the mean BS (a) and the mean Area (b) over all the RWSs. It clearly demonstrates that the ensemble forecasts are preferable to the deterministic ones. The Area increases as the size of the ensemble increases, but the differences between the ensembles with 100 and 500 members were negligible. The evaluation using the BS agrees well with
Fig. 3. The VRHs of ensembles ENS1 and ENS2 with 100 members for station K013 and for lead times FT = 60 (1st row), 180 (2nd row) and 300 (3rd row) min. The data are divided into 10 classes, and the vertical axes show the frequencies of their occurrences. The left and right classes show the relative numbers of measured road surface temperatures below the ensemble minimum and above the ensemble maximum, respectively. The remaining data were homogeneously divided into 8 inner classes (columns).
38
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
Fig. 4. Probability Prb that the observed values are within the interval determined from the distribution of the ensemble values that should include P% of the values. The probabilities Prb are shown for the ENS1 and ENS2 ensembles and ensemble sizes of 10, 30, 100 and 500 for lead times up to 6 h.
the results obtained by the Area measure, but the decrease with ensemble size is smaller and is completely negligible for more than 100 members (not shown here). Fig. 5 again confirms that the inclusion of N into the ensemble is important and significantly increases the accuracy in terms of the Area and BS.
3.3. Example of an ensemble forecast Two examples of ensemble forecasts illustrating typical cases are shown in Fig. 6. The left column of Fig. 6 shows the sequence of forecasts for RWS K013 and the right column shows the forecasts for RWS P004.
Fig. 5. Subpanels (a) and (b) show the mean BS and the average ROC area (Area), respectively, over all the RWSs as a function of the lead time for the deterministic forecast (det) and ensemble forecasts. The results of ensembles ENS1 and ENS2 with 10 and 100 members for (a) and 10, 30, 100 and 500 members for (b) are denoted by solid lines (C in the legend) and dashed lines, respectively.
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
39
Fig. 6. The ensemble forecasts of RST in °C for RWSs K013 (left column) and P004 (right column) as a function of the lead times in h. The full and dashed curves show the deterministic forecasts and observed data, respectively. The grey areas denote the intervals for probabilities x = 99%, 90%, 70% and 50% that the RST will be within the interval denoted. The start times of the forecasts are indicated in the titles in a year, month, day and hour (2 digits) format.
The left column displays a case with a successful deterministic forecast that does not differ from the observed values by more than approximately 1 °C. The grey areas graphically depict the expected intervals of the RST values for the corresponding probabilities. For example, there is a probability of 50% that the RST will be in the darkest area for the given lead time. The uncertainty in the forecast naturally increases as lead time increases. A less successful forecast is presented in the right column. In that case, the RST forecast was difficult because a frontal system was passing over the Czech Republic with intermittent intensive rain, which created a water layer on the road surface that froze when the rain stopped, and the sky began to clear (Sokol et al., 2014). This was not correctly forecasted by the ALADIN-CZ NWP model, and because errors in the precipitation forecasts were not taken into account when the ensemble was generated, the ensemble spread was too narrow, and the observed values were well outside the borders of the 99% interval of values for longer lead times. The ensemble approach is especially suitable for predicting road surface conditions. The Metro model recognizes 8 road surface states (dry, wet, ice/snow, mixture of water and ice/snow, dew, melting ice/snow, frost and icing rain) and “not known”, which depend not only on the temperature but also on other characteristics of the road surface and near-surface air. Because these dependencies are very complicated, a deterministic forecast that provides just one state is insufficient to evaluate the possible risk of dangerous surface conditions. The ensemble forecast is very useful for this because each member of the ensemble predicts other conditions, and if the number of members is large enough, a simple statistical evaluation of the ensemble of predictions can be used to quantitatively evaluate the risk of the occurrence of specific road surface conditions. The
ensemble forecast is illustrated in Fig. 7. The displayed forecasts correspond to the RST forecasts in Fig. 6 (right column). Fig. 7 shows that single ensemble members yield various surface conditions when the RST is close to 0 °C and alternating rain and clear sky conditions occur. Currently, we are unable to objectively verify these ensemble forecasts of road surface conditions because we do not have observed data.
4. Conclusions This paper introduced the METRo-CZ model, which is used to forecast RST and road surface conditions. The basic innovation of the model compared to the original METRo model is the inclusion of ensemble forecasts that allows an expression of the uncertainties in the forecasts. To our knowledge, the ensemble approach based on the proposed method for generating ensembles using a deterministic NWP forecast has not been previously applied to models that forecast road surface temperatures. The derivation and description of a technique for generating ensembles is the core of the article, together with analyses of ensemble forecasts from the perspective of the ensembles' ability to represent real uncertainties in the forecast and the usefulness of the ensemble results for practical application. We plan to utilize ensemble NWP forecasts when they are provided by the CHMI and to study how the error covariance matrix C depends on various weather types, e.g., very changeable weather versus stable weather, and if classifying such weather can improve the forecast accuracy. The obtained results are summarized as follows.
40
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
Fig. 7. Statistically processed ensemble forecasts of road weather conditions. The vertical axis is probability, the colours (greyscale) indicate the road conditions, and the horizontal axis is the lead time. The Metro model recognizes 8 road surface states (dry, wet, ice/snow, mixture of water and ice/snow, dew, melting ice/snow, frost and icing rain) and “not known”.
1) The applied ensemble technique is based on covariance analyses of errors in ALADIN-CZ NWP forecasts. The technique is focused on variables with Gaussian error distributions, which is not the case for wind speed and total cloud cover; however, in our opinion, this is not an essential problem. The RST forecasts do not significantly depend on the wind speed, the errors in forecasts of cloud cover are quite symmetric, and the violation of normality is not very significant. 2) To perform the covariance analyses of errors in cloud cover forecasts, we used the SAF algorithm PGE01 Cloud Mask and found that it agrees well with subjective observations of cloud cover at synoptic stations. 3) Temperature, humidity and wind speed in the near-surface air layer are not sufficient to reasonably describe the uncertainties in RST forecasts. Radiation fluxes, which are derived from the total cloud cover, are crucial for RST development, and this quantity should be included in the ensemble generation. Even after doing so, however, the ensemble is too confident (the spread is not sufficient) and cannot describe the natural uncertainty in the forecasts. 4) There are several primary reasons for the insufficient ensemble spread and consequent underestimation of the forecast uncertainty. One reason is that NWP model errors in forecasted rain and snow are not considered when generating ensembles. It is clear that the road surface conditions (dry, wet, or covered by snow) significantly influence the development of the RST. However, we did not include these
characteristics in the ensembles owing to their distinct nonGaussian character of these characteristics. Other important reasons for errors are of a technical nature, for example, incorrect information on the sky view factor, unavailable or inaccurate information on the current road surface conditions, and inaccurate information on the technological parameters of roads. We can expect that these data shortcomings will at least be partly resolved in the future. 5) Although the ensemble technique underestimates the uncertainties in RST forecasts, it can be used in decision-making by road authorities who provide winter road maintenance. The outputs can be subjectively “corrected”, or a technique such as Model Output Statistics can be applied to objectively correct the forecasts. 6) In our opinion, the ensemble approach is especially suitable for predicting road surface conditions because it can express the probability of the occurrence of various road surface states, which deterministic forecasts do not allow. Unfortunately, we were unable to objectively verify the ensemble forecasts because we currently do not have required observation data. Preliminary results are promising. Acknowledgements This work was supported by the Czech Science Foundation grant GACR 13-34856S and TA01031509. The data were kindly provided by the Czech Hydrometeorological Institute. We acknowledge comments of reviewers which significantly improved the manuscript.
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
Appendix A. Generation of a Gaussian random vector u having a known covariance matrix C = cov(u, uT) A Gaussian random vector u that has a known covariance matrix C = cov(u, uT) can be calculated by u ¼ Lz;
ðA1Þ
where L is a lower triangle matrix obtained by the Cholesky decomposition of C, i.e., C ¼ LLT ;
ðA2Þ
and z = (z1, z2, …, zn)T is an independent Gaussian random vector with zero mean and unit variance (Huynh et al., 2008). In our application, we first constructed the matrix C defined by Eqs. (2) and (3) using historical data T (Eq. (1)) and calculated its Cholesky decomposition matrix L (Eq. (A2)). Second, we generated a Gaussian random number with zero mean and unit random z⁎ by a standard algorithm (e.g., Huynh et al., 2008; Wilks, 2006). Finally, we calculated T⁎ = Lz⁎ and obtained a forecast for one member of the ensemble by integrating the METRo-CZ model using data from vector T⁎ rather than the forecasted NWP model data. By repeatedly generating random z⁎ values, we obtained the whole ensemble of forecasts. In our experiments, the first ensemble member corresponded to the original NWP model data. Appendix B. Total cloud cover estimated from satellite data Measurements by Meteosat Second Generation (MSG) satellites were used to estimate the total cloud cover (N). MSG data are available every 15 min, and their horizontal resolution is approximately 5 × 5 km over Central Europe. To calculate N, the SAF algorithm PGE01 Cloud Mask is applied (SAFNWC, 2013; Badescu et al., 2016). It classifies cloud cover into the following 3 classes for each satellite grid point: cloud-free, cloud-contaminated (partly cloudy or semitransparent) and cloud-filled (opaque clouds) to which we assign 0, 4 and 8 octas of N, respectively. The obtained N is interpolated from the satellite grid into the new grid using the nearest neighbour technique with a horizontal resolution of 1 km. Finally, N is calculated for each new grid G as the mean value of the grid values of N from the surrounding
Fig. B1. A histogram of the differences between N observed by observers at synoptic stations and calculated from the MSG data using m = 7 and m = 9.
41
m × m grids, where m is an odd number, and grid G is in the centre of the m × m grids. This procedure simulates how N is subjectively estimated by observers. Based on tests, we used m = 9 in this study. The results of the algorithm were compared with N from SYNOP station reports, which contain the subjectively determined N. The comparison included data from 17 synoptic stations from the Czech Republic and from 4 winter months for the years 2012–2013. A histogram of the differences is shown in Fig. B1 and indicates that m = 9 and m = 7 yielded very similar results. Taking into account the uncertainty in the subjectively determined N, Fig. B1 indicates a good agreement and acceptable asymmetry between the two estimations. Similar column heights for values of x = 2, 3 and 4 that do not correspond to Gaussian distributions are caused by cases when the subjective observations reported N = 7 and N = 8, but the applied algorithm recognized semitransparent cloudiness, which was classified as N = 4. References Badescu, V., Paulescu, M., Brabec, M., 2016. Reconstruction of effective cloud fields geometry from series of Sunshine Number. Atmos. Res. (accepted). Bouallegue, Z.B., Theis, S.E., Gebhardt, C., 2013. Enhancing COSMO-DE ensemble forecasts by inexpensive techniques. Meteorol. Z. 22, 49–59. Bouilloud, L., Martin, E., Habets, F., Boone, A., Moigne, P.L., Livet, J., Marchetti, M., Foidart, A., Franchistéguy, P., Morel, S., Noilhan, J., Pettre, P., 2009. Road surface condition forecasting in France. J. Appl. Meteorol. Climatol. 48, 2513–2527. Crevier, L.-P., Delage, Y., 2001. METRo. A new model for road-condition forecasting in Canada. J. Appl. Meteorol. 40, 2026–2037. Greenfield, T.M., Takle, E.S., 2006. Bridge frost prediction by heat and mass transfer methods. J. Appl. Meteorol. Climatol. 45, 517–526. Huynh, H.T., Lai, V.S., Soumaré, I., 2008. Stochastic Simulation and Applications in Finance With MATLAB Programs. John Wiley & Sons Ltd, p. 338 (ISBN 978–0–470-72538-2). Jansson, C., Almkvist, E., Jansson, P.-E., 2006. Heat balance of an asphalt surface: observations and physically-based simulations. Meteorol. Appl. 13, 203–212. Kangas, M., Heikinheimo, M., Hippi, M., 2015. RoadSurf: a modelling system for predicting road weather and road surface conditions. Meteorol. Appl. 22, 544–553. Leutbecher, M., Palmer, T.N., 2008. Ensemble forecasting. Weather Forecast. 227, 3515–3539. Linden, S.K., Petty, K.R., 2008. The Use of METRo (Model of the Environment and Temperature of the Roads) in Roadway Operation Decision Support Systems. 24th Conference on IIPS. NCAR, Boulder, CO. Rutz, J.J., Gibson, C.V., 2013. Integration of a road surface model into NWS operations, B. Am. Meteor. Society 94. SAFNWC, 2013. Algorithm Theoretical Basis Document for “Cloud Products” (CMa PGE01 v3.2, CT-PGE02 v2.2 & CTTH-PGE03 v2.2). Shao, J., Lister, P.J., 1996. An automated nowcasting model of road surface temperature and state for winter road maintenance. J. Appl. Meteorol. 35, 1352–1361. Sokol, Z., 2003. MOS based precipitation forecasts for river basins. Wea. Forecast 18, 769–781. Sokol, Z., Zacharov, P., Sedlák, P., Hošek, J., Bližňák, V., Chládová, Z., Pešice, P., Škuthan, M., 2014. First experience with the application of the METRo model in the Czech Republic. Atmos. Res. 143, 1–16. Wilks, D.S., 2006. Statistical Methods in the Atmospheric Sciences. second ed. Academic Press, p. 627 (ISBN : 9780080456225). Zakšek, K., Oštir, K., Kokalj, Ž., 2011. Sky-view factor as a relief visualization technique. Remote Sens. 3, 398–415.
Contents lists available at ScienceDirect
Atmospheric Research journal homepage: www.elsevier.com/locate/atmosres
Ensemble forecasts of road surface temperatures Zbyněk Sokol a,⁎, Vojtěch Bližňák a, Pavel Sedlák a, Petr Zacharov a, Petr Pešice a, Miroslav Škuthan b a b
Institute of Atmospheric Physics, Czech Academy of Sciences, Boční II 1401, 141 31, Prague, Czech Republic Czech Hydrometeorological Institute, Na Šabatce 2050/17, 143 06, Prague, Czech Republic
a r t i c l e
i n f o
Article history: Received 15 March 2016 Received in revised form 7 December 2016 Accepted 16 December 2016 Available online 22 December 2016 Keywords: Road weather forecast Road surface temperature Ensemble prediction
a b s t r a c t This paper describes a new ensemble technique for road surface temperature (RST) forecasting using an energy balance and heat conduction model. Compared to currently used deterministic forecasts, the proposed technique allows the estimation of forecast uncertainty and probabilistic forecasts. The ensemble technique is applied to the METRo-CZ model and stems from error covariance analyses of the forecasted air temperature and humidity 2 m above the ground, wind speed at 10 m and total cloud cover N in octas by the numerical weather prediction (NWP) model. N is used to estimate the shortwave and longwave radiation fluxes. These variables are used to calculate the boundary conditions in the METRo-CZ model. We found that the variable N is crucial for generating the ensembles. Nevertheless, the ensemble spread is too small and underestimates the uncertainty in the RST forecast. One of the reasons is not considering errors in the rain and snow forecast by the NWP model when generating ensembles. Technical issues, such as incorrect sky view factors and the current state of road surface conditions also contribute to errors. Although the ensemble technique underestimates the uncertainty in the RST forecasts, it provides additional information to road authorities who provide winter road maintenance. © 2016 Elsevier B.V. All rights reserved.
1. Introduction Forecasts of road surface conditions are highly important for road maintenance services. These forecasts usually apply models based on solving energy balance and heat conduction equations (e.g., Bouilloud et al., 2009; Crevier and Delage, 2001; Greenfield and Takle, 2006; Jansson et al., 2006; Kangas et al., 2015; Rutz and Gibson, 2013; Shao and Lister, 1996). The energy balance equation calculates boundary conditions for the heat conduction equations, the accuracy of which has crucial impacts on forecasts of road surface temperature (RST) and road surface conditions. The boundary conditions depend on the current and future states of the atmosphere, primarily forecasts of shortwave radiation flux, longwave radiation flux, air temperature and humidity, wind speed and amount and type (rain, snow) of precipitation, which are required to calculate energy fluxes between the road surface and the near surface atmosphere. These prognostic variables are obtained from numerical weather prediction (NWP) model forecasts. The boundary conditions are influenced by various errors that are mainly caused by inaccurate NWP model forecasts and parameterizations applied to calculate surface fluxes from the NWP model variables. An improper consideration of shade effects may also contribute to errors. Errors also arise from the energy balance model assumption that
⁎ Corresponding author. E-mail address: [email protected] (Z. Sokol).
http://dx.doi.org/10.1016/j.atmosres.2016.12.010 0169-8095/© 2016 Elsevier B.V. All rights reserved.
the entire surface of the Earth is composed of a horizontally homogeneous road, which is not met in reality. The model that describes heat conduction in the road body also contributes to the overall error, particularly because the characteristics of the road layers (their thicknesses, thermal conductivities and capacities) may not be accurate. However, the contribution of these errors is not very important, especially for shorter lead times, because of the slow response of the surface temperature to the subsurface characteristics. Because of these uncertainties, deterministic forecasts do not provide sufficient information for road maintenance services to make optimal decisions. It is important to have additional information regarding the accuracy of the deterministic forecast in various forms, for example, the probability that the forecasted value will be in a given interval of values. The problem of the uncertainty of deterministic forecasts is also confronted in other meteorological forecasts. To solve this problem, an ensemble technique was introduced, developed and applied to weather forecasting over the last two decades. The ensemble technique consists of generating a set of deterministic forecasts that represent possible future states, and the problem of limited predictability and uncertainty can be addressed (e.g., Leutbecher and Palmer, 2008) by statistical processing of a set of obtained forecasts. The results of the ensemble technique depend fundamentally on how individual members of the ensemble are generated. In meteorological applications, ensembles are usually generated by varying the initial and/or boundary conditions; however, changes in model parameterizations and/or model parameter values have also been frequently used (e.g., Bouallegue et al., 2013).
34
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
In this paper, an ensemble technique for RST forecasting is presented using the METRo-CZ model (Sokol et al., 2014). The ensemble technique is used to develop a probabilistic forecast that the RST will be below a given threshold (0 °C), which is important information for road maintenance services. We applied a method for generating ensembles that simultaneously modifies forecasted air temperature and humidity 2 m above the ground (T2 m, H2 m), wind speed at 10 m (W10 m) and the total cloud cover N in octas obtained from the NWP model. N is used to estimate the shortwave and longwave radiation fluxes. To our knowledge, this is the first application of this or a similar approach. It is noteworthy that the NWP model used in this study (Section 2) does not provide an ensemble forecast. Therefore, an ensemble technique based on the deterministic forecast was applied. This paper is divided into four sections. Section 2 describes the METRo-CZ model, including the data that are used to prepare the initial and boundary conditions. The section also contains a description of the ensemble technique. An evaluation and discussion of the ensemble technique are presented in Section 3, and the conclusions are presented in Section 4. Details concerning the ensemble technique are provided in the Appendices.
2. The METRo-CZ model 2.1. Model description The Institute of Atmospheric Physics CAS in Prague and the Czech Hydrometeorological Institute (CHMI) developed and currently run the METRo-CZ model for forecasts of RST and road surface conditions. The model is based on solving the energy balance and heat conduction equations. It originated as the Model of the Environment and Temperature of Roads (METRo), which was originally developed by Environment Canada (Crevier and Delage, 2001). The METRo model documentation and source code are freely available at http:// home.gna.org/metro/. METRo-CZ is a 1D model that calculates independent forecasts at single points. It uses measured data from road weather stations (RWS) and forecasted data obtained by integration of the ALADIN NWP model as inputs to prepare initial and boundary conditions (Sokol et al., 2014). A detailed description of the data employed in this study is provided in Section 2.2. METRo-CZ and the original METRo comprise 3 basic steps (Crevier and Delage, 2001; Linden and Petty, 2008; Sokol et al., 2014). In the 1st step, which is called initialization, the vertical profile of the ground temperature is calculated. In the 2nd step, called coupling, radiative fluxes forecasted by the NWP model are corrected by an additive correction coefficient. The aim is to ensure that forecasted and measured RSTs are close at the beginning of the METRo-CZ forecast (TMET). The coefficient value is calculated iteratively. In each iteration, the correction coefficient is increased or decreased depending on the sign of the difference between the forecasted and measured RST and the magnitude of the change decreases in iteration steps. The process is stopped when the difference between the RSTs is less than 0.05 °C. The length of the coupling period depends on the beginning of the METRo-CZ forecast and is at least 4 h. If the forecast starts at TMET = 4, 5, …, 9 UTC, then the coupling begins at TNWP = 00 UTC using data from the ALADIN run starting at TNWP. The initialization, which only uses measured data from the RWSs, is run from TNWP-6 h to TNWP. Similarly, if the forecast starts at 10, …, 15 UTC, then the coupling begins at TNWP = 06 UTC, etc. The current version of METRo-CZ differs from METRo in various aspects. The original code was completely rewritten in Fortran 90, and a large majority of the model parameter values are now model inputs that can be easily changed. The data pre-processing was modified and tailored to the data used in the CHMI. The physical core of the model equations is the same as in METRo, and the basic differences between the METRo-CZ and the original METRo model will now be discussed.
METRo-CZ uses N to calculate the shortwave and longwave radiation fluxes. We found that the radiation fluxes calculated directly by the ALADIN NWP model yielded poorer RST forecasts than when N was utilized. We also found that the relationship between N and the shortwave radiation fluxes (Table 2 in (Crevier and Delage, 2001)) used in METRo overestimates global radiation in comparison with measurements at the Czech synoptic stations. Therefore, using measured global radiation and observed N we derived and applied a different relationship (Table 1). METRo-CZ focuses on nowcasting (i.e., very short-range forecasts) with lead times up to 6 h. It uses NWP data with time steps of 1 h and all available data from RWSs and calculates new forecasts every hour.
2.2. Model input data As previously stated, METRo-CZ uses prognostic data from the ALADIN NWP. Specifically, it employs the ALADIN-CZ version of the model, which is operated by the CHMI. ALADIN-CZ is a hydrostatic model that is integrated at a horizontal resolution of approximately 5 km; the integration begins four times per day (00, 06, 12 and 18 UTC), and its lead time is 54 h. The following forecasted data are required: air temperature T2 m and humidity H2 m at 2 m, wind speed W10 m at 10 m, model pressure at the model orography, cloud cover N, and accumulated precipitation and type (rain/snow). These model data are available 4 h after the start of NWP model run. METRo-CZ uses the following measurements from the RWSs: air temperature T2 m and humidity H2 m at 2 m, wind speed W10 m at 10 m, the road weather conditions (a code indicating dry or wet surface), RST, and temperature at 30 cm below the road surface (T30). If any T2 m, H2 m or W10 m data are not available, their values are estimated by linear interpolation of ALADIN-CZ data over time. If information on the road weather conditions is missing, a dry surface is assumed. If RST is not available the measurement term is not used. In this study, we used 27 RWSs, which were selected to cover various parts of the western and central Czech road network (highways and 1st class main roads) under various road conditions (lowlands and mountains). Their positions (marked with crosses) with RWS codes and a map of the highway network over the orography (blue lines indicate open highways) in the western and central Czech Republic are shown in Fig. 1. The RWSs performed measurements every 20 min. The data from each station were available for more than 96% of the measurement period, and the quality check of T2 m, H2 m, W10 m and RST removed less than 1% of the measurements. The quality check compared neighbouring measurements in the time sequences and removed unrealistically large differences in the measurements. Thus, less than 1% of the measurements were removed. The measured road weather conditions were not assessed. These data were available at each RWS for at least 94% of the measurement period, with the exceptions of RWSs J005 and P019, which had measurements for only 46% and 40% of cases, respectively.
Table 1 Original and new coefficients for the shortwave radiation fluxes (see Table 2 in [Crevier and Delage, 2001]). Cloud cover (octas)
D1 (original)
D1 (new)
0 1/8 2/8 3/8 4/8 5/8 6/8 7/8 8/8
1.00 0.97 0.94 0.89 0.85 0.80 0.71 0.65 0.33
0.94 0.91 0.90 0.87 0.83 0.75 0.67 0.48 0.25
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
35
Fig. 1. The highway network (thick blue lines indicate open highways) in the western and central Czech Republic. The positions of the road weather stations used in this study are denoted by crosses with RWS codes. The orography is shown in the background. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
2.3. Ensemble generation 2.3.1. General information on the ensembles The purpose of the ensemble technique is to statistically describe the uncertainty in the deterministic forecasts and possibly refine the deterministic predictions using the obtained statistical characteristics. Forecast error is generally dependent on model errors (caused, for example, by the application of simplified equations and inaccurate numerical solutions) and errors in the input model data (initial and boundary conditions). In our experience, errors in forecasted atmospheric data that represent boundary conditions in METRo-CZ and describe the state of the near-surface layer of the atmosphere significantly influence RST forecasts (Sokol et al., 2014). Errors related to temperature, humidity, wind speed and, especially, radiation fluxes and the occurrence of precipitation (rain or snow) dominate model errors caused by the use of simplified equations. Therefore, we focused on creating ensembles using atmospheric temperature, humidity, wind speed and total cloud cover, and neglected errors and inaccuracies produced by the METRo-CZ equations. We also did not consider errors in the precipitation forecasts because their statistical characteristics are complicated and it is not clear how to use them in the context of the applied ensemble technique. In our calculations, we neglected the shading effect in spite of its importance due to lacking the required data at our disposal at present. We calculated an ensemble of forecasts by generating an ensemble set of boundary data and applied METRo-CZ to this data set. The set of
boundary data was prepared as follows: For forecasts starting at TMET and with projection time δ (in our case, δ = 6 h), we had available hourly NWP forecasts of T2 m, H2 m, W10 m and N, where the lead time t1 corresponded to TMET, and t2 corresponded to TMET + δ. To simplify the notation, we will introduce a vector T of length 4δ containing hourly forecasts for δ lead times and 4 variables: T ¼ ðT2
m ðt1 Þ; H2 m ðt1 Þ; W10 m ðt1 Þ; Nðt1 Þ; …; T2 m ðt2 Þ; H2 m ðt2 Þ; W10 m ðt2 Þ; Nðt2 ÞÞ:
ð1Þ We used an available set of historical NWP forecasts and observations, which we denote Ti,j and Oi,j, where i = 1, …,n refer to an individual forecast or observation for single hours, and j = 1, …,4δ (we repeat that δ = 6 h in our case) refer to components of the vector T. Using the historical data, we calculated the covariance matrix C: C ¼ ck;m ; k; m ¼ 1; …; 4δ;
ð2Þ
where ck;m ¼
n 1X T −Oi;k T i;m −Oi;m n i¼1 i;k
ð3Þ
and the indices of T and O stand for an individual forecast or observation and vector component, respectively. The matrix C was used to describe the error structure of the boundary conditions, i.e., of vector T. To obtain
36
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
the ensemble of forecasts, we generated an ensemble (set) of vectors T whose components were selected such that they fulfilled the covariance relationships described by Eq. (2). The procedure for generating these vectors is described in Appendix A. The covariance matrix C is the core of the applied method. Because METRo-CZ is a 1D model calculated separately for each RWS, it is natural to prepare matrices C for each station individually. Temperature, humidity and wind speed are measured at all RWSs; thus, the observational data in Eq. (3) are available. Because N is, however, not observed at those stations, we instead used satellite data from the Meteosat Second Generation and an algorithm based on the SAF algorithm PGE01 Cloud Mask to estimate N. The algorithm is described in Appendix B. Fig. 2 shows an example of the generated ensembles of T2 m, H2 m, W10 m and N. The ensemble had 50 members and was prepared for station P002 for 15 November 2012, 06 UTC. It should be mentioned that the humidity in (1–3) is expressed as a dew point temperature and then transformed into specific humidity, as required by METRo-CZ, using the perturbed temperature and forecasted pressure. Fig. 2 clearly demonstrates the sizes of the uncertainties in the forecasts of T2 m, H2 m, W10 m and N. N (subpanel d) is of particular concern; it seems to behave chaotically, but this is only partly true when considering the grey area, where 50% of the values are located.
2.3.2. Specific application of the ensemble technique and its limitations We prepared three versions of the ensemble application. The first version perturbed only T2 m, the second perturbed T2 m, H2 m and W10 m, and the third perturbed all the variables (T2 m, H2 m, W10 m and N). The construction of the matrix C differed among the versions. Because the first version always yielded the poorest scores in terms of Relative Operating Characteristics, Brier score and verification rank histogram (Wilks, 2006), which are used to evaluate ensembles in Section 3, it is not further evaluated in this paper. The second and third versions differed by the inclusion of N, which is an important factor in the uncertainty of the METRo-CZ prediction. However, the determination of N may be affected by significant errors that influence the covariance matrix C and the resulting ensembles, and thus may deteriorate ensemble attributes. Therefore, both versions with and without N were tested. In this paper, we used the covariance matrix C calculated using data from the winter season of 2012–2013 (October to March). Although this time period included the evaluation period, it had a small impact on the results, which we verified by gradually excluding just one month (from November to February) of data from the period used to calculate C and compared the METRo-CZ forecasts with those obtained using all data. The overall results were quite similar (not shown here). Therefore, we assume that the calculated covariances were sufficiently robust and the overlap of the time periods, caused by data availability, did not affect significantly the obtained results. In our subjective opinion, one winter season, typically the latest one, is the most suitable length of data for calculating C. The reason is that NWP models are frequently changed, and major changes, which may significantly influence the covariances, are usually performed once per year. As follows from Eq. (1), the values in matrix C depend on the difference between the start of the NWP forecasts and the METRo-CZ models (index j1). This dependence should also include the time of day because the NWP model forecast errors have diurnal variations. Therefore, we calculated C for each hour in a day. The basic shortcoming of the described ensemble technique is that the covariance error structure implicitly assumes that the errors in the T components are normally distributed. This assumption is reasonable for T2 m and H2 m but is not valid for W10 m and N. The violation of normality is not important for W10 m because the W10 m values have a small impact on the METRo-CZ forecasts. N, however, has a significant impact on the METRo-CZ forecasts. The distribution of the N errors is not normal but is approximately symmetric, with a significant peak near zero (Fig. B1); thus, we can expect that the violation of normality may not cause significant errors in the METRo-CZ forecasts. 3. Forecast evaluation 3.1. Implementation of METRo-CZ The RST and road weather conditions were calculated using METRoCZ for the period from 1 November 2012 to 28 February 2013. The forecasts were calculated at the 27 RWS introduced in Section 2.2. The forecasts started every hour and were projected out 6 h. For each hour, we performed an ensemble of 10, 30, 100, 250 and 500 forecasts. The first forecast used the T2 m, H2 m, W10 m and N data forecasted by ALADINCZ, and the remaining forecasts used the perturbed data. 3.2. Evaluation of the ensemble forecasts
Fig. 2. An example of the generated ensembles of T2 m (TA), H2 m (Q), W10 m (W) and N. The ensemble had 50 members and was prepared for station P002 for 15 November 2012, 06 UTC. The filled area in the subpanel d) contains 50% of the N values.
The evaluation was carried out via two closely related directions. First, we evaluated the properties of two types of ensembles generated using perturbations in (i) T2 m, H2 m, W10 m and N (ENS1) and (ii) T2 m, H2 m and W10 m (ENS2). The second direction focused on a discussion of the practical usefulness of the ensembles to express uncertainties in the predictions and forecasts in comparison with the deterministic forecast. Ideally, the true value should only be one additional member of the ensemble and indistinguishable from the members. This attribute of the
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
ensemble can be evaluated using the verification rank histogram (VRH; Wilks, 2006). The VRH for a scalar variable is calculated as follows: for each forecast, the ensemble values are ordered from the smallest to the largest value and the order of the observed values in the ensemble is determined. If the observed value is lower than the smallest ensemble value, the order is 1. If the observed value is between the nth and the (n + 1)th ensemble values, then the order is n + 1, etc. The VRH histogram shows frequencies of the calculated orders of observed values. The VRH of the ideal ensemble should have an approximately uniform distribution. In our plots, we used 10 classes in the histograms. The 1st and 10th classes contained the relative number of observed data values lower than the minimum ensemble values and greater than the maximum ensemble values, respectively. Classes 2 to 9 contained the relative numbers of observed data values in 8 equally wide intervals between ensemble minimum and maximum. Fig. 3 shows an example of the VRH of ensembles ENS1 and ENS2 for K013 RWS and for lead times of 60, 180 and 300 min for 300 ensemble members. The other RWSs had similar VRHs. The high rightmost and especially high leftmost columns indicate an under dispersion of the ensembles, meaning that the ensembles are overconfident and do not cover the entire variability in the observations. The high relative number of measurements that are lower than the minimum of the ensembles can be explained mainly by two factors that were not considered when generating the ensembles. First, rain and/or snow often lead to cooler surface temperatures than does a clean dry road surface, and second, the sky view factor, meaning a parameter that corresponds to the portion of the visible sky limited by relief (Zakšek et al., 2011), is not properly taken into account, which may result in shadowing and unexpected significant temperature drops during the day. The first factor is also related to the fact that the state of the road surface is frequently not known. The comparison of the VRHs for ENS1 and ENS2 clearly indicates that the cloud cover N is an important source of uncertainty in the forecasts. N mainly affected
37
the first class on the histogram, resulting in a higher RST forecast of METRo-CZ compared to the RWS measurements. For given probabilities P = 50%, 70% and 95%, RST intervals were calculated using the ensemble forecasts. For example, for P = 50%, the interval was determined as between the 25% and 75% percentiles of the forecasts for a given term. Fig. 4 shows the results obtained for all the stations as a function of the lead times for ENS1 and ENS2 with various ensemble sizes. Consistent with Fig. 3, it is clear that ensemble ENS1 underestimated the spread. The significantly poorer results for ENS2 confirm the important role of N in the construction of ensembles. From the viewpoint of ensemble sizes, it follows from Fig. 4 that the size should be at least 100 members for higher probabilities, e.g., P = 95%. The obtained results can be potentially utilized to extend the forecasted intervals by applying empirical corrections to better match the measurements. The Model Output Statistics method (MOS) e.g., (Sokol, 2003) is suitable for that task. To predict road surface conditions, it is important to forecast whether the surface temperature will be lower or higher than 0 °C; this forecast may be categorical or probabilistic, the latter of which is preferable. These types of forecasts are frequently evaluated using Relative Operating Characteristics (ROC) and derived characteristics (Wilks, 2006) that include the area under the ROC curve (Area) or by the Brier score (BS) (Wilks, 2006), which is similar to the standard mean-squareerror calculated using probabilities or binary values. For a perfect forecast, Area = 1; for a random forecast, Area = 0.5. The Area value shows the potential usefulness of the analysed forecast and allows a comparison of the deterministic and ensemble forecasts. Fig. 5 shows the mean BS (a) and the mean Area (b) over all the RWSs. It clearly demonstrates that the ensemble forecasts are preferable to the deterministic ones. The Area increases as the size of the ensemble increases, but the differences between the ensembles with 100 and 500 members were negligible. The evaluation using the BS agrees well with
Fig. 3. The VRHs of ensembles ENS1 and ENS2 with 100 members for station K013 and for lead times FT = 60 (1st row), 180 (2nd row) and 300 (3rd row) min. The data are divided into 10 classes, and the vertical axes show the frequencies of their occurrences. The left and right classes show the relative numbers of measured road surface temperatures below the ensemble minimum and above the ensemble maximum, respectively. The remaining data were homogeneously divided into 8 inner classes (columns).
38
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
Fig. 4. Probability Prb that the observed values are within the interval determined from the distribution of the ensemble values that should include P% of the values. The probabilities Prb are shown for the ENS1 and ENS2 ensembles and ensemble sizes of 10, 30, 100 and 500 for lead times up to 6 h.
the results obtained by the Area measure, but the decrease with ensemble size is smaller and is completely negligible for more than 100 members (not shown here). Fig. 5 again confirms that the inclusion of N into the ensemble is important and significantly increases the accuracy in terms of the Area and BS.
3.3. Example of an ensemble forecast Two examples of ensemble forecasts illustrating typical cases are shown in Fig. 6. The left column of Fig. 6 shows the sequence of forecasts for RWS K013 and the right column shows the forecasts for RWS P004.
Fig. 5. Subpanels (a) and (b) show the mean BS and the average ROC area (Area), respectively, over all the RWSs as a function of the lead time for the deterministic forecast (det) and ensemble forecasts. The results of ensembles ENS1 and ENS2 with 10 and 100 members for (a) and 10, 30, 100 and 500 members for (b) are denoted by solid lines (C in the legend) and dashed lines, respectively.
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
39
Fig. 6. The ensemble forecasts of RST in °C for RWSs K013 (left column) and P004 (right column) as a function of the lead times in h. The full and dashed curves show the deterministic forecasts and observed data, respectively. The grey areas denote the intervals for probabilities x = 99%, 90%, 70% and 50% that the RST will be within the interval denoted. The start times of the forecasts are indicated in the titles in a year, month, day and hour (2 digits) format.
The left column displays a case with a successful deterministic forecast that does not differ from the observed values by more than approximately 1 °C. The grey areas graphically depict the expected intervals of the RST values for the corresponding probabilities. For example, there is a probability of 50% that the RST will be in the darkest area for the given lead time. The uncertainty in the forecast naturally increases as lead time increases. A less successful forecast is presented in the right column. In that case, the RST forecast was difficult because a frontal system was passing over the Czech Republic with intermittent intensive rain, which created a water layer on the road surface that froze when the rain stopped, and the sky began to clear (Sokol et al., 2014). This was not correctly forecasted by the ALADIN-CZ NWP model, and because errors in the precipitation forecasts were not taken into account when the ensemble was generated, the ensemble spread was too narrow, and the observed values were well outside the borders of the 99% interval of values for longer lead times. The ensemble approach is especially suitable for predicting road surface conditions. The Metro model recognizes 8 road surface states (dry, wet, ice/snow, mixture of water and ice/snow, dew, melting ice/snow, frost and icing rain) and “not known”, which depend not only on the temperature but also on other characteristics of the road surface and near-surface air. Because these dependencies are very complicated, a deterministic forecast that provides just one state is insufficient to evaluate the possible risk of dangerous surface conditions. The ensemble forecast is very useful for this because each member of the ensemble predicts other conditions, and if the number of members is large enough, a simple statistical evaluation of the ensemble of predictions can be used to quantitatively evaluate the risk of the occurrence of specific road surface conditions. The
ensemble forecast is illustrated in Fig. 7. The displayed forecasts correspond to the RST forecasts in Fig. 6 (right column). Fig. 7 shows that single ensemble members yield various surface conditions when the RST is close to 0 °C and alternating rain and clear sky conditions occur. Currently, we are unable to objectively verify these ensemble forecasts of road surface conditions because we do not have observed data.
4. Conclusions This paper introduced the METRo-CZ model, which is used to forecast RST and road surface conditions. The basic innovation of the model compared to the original METRo model is the inclusion of ensemble forecasts that allows an expression of the uncertainties in the forecasts. To our knowledge, the ensemble approach based on the proposed method for generating ensembles using a deterministic NWP forecast has not been previously applied to models that forecast road surface temperatures. The derivation and description of a technique for generating ensembles is the core of the article, together with analyses of ensemble forecasts from the perspective of the ensembles' ability to represent real uncertainties in the forecast and the usefulness of the ensemble results for practical application. We plan to utilize ensemble NWP forecasts when they are provided by the CHMI and to study how the error covariance matrix C depends on various weather types, e.g., very changeable weather versus stable weather, and if classifying such weather can improve the forecast accuracy. The obtained results are summarized as follows.
40
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
Fig. 7. Statistically processed ensemble forecasts of road weather conditions. The vertical axis is probability, the colours (greyscale) indicate the road conditions, and the horizontal axis is the lead time. The Metro model recognizes 8 road surface states (dry, wet, ice/snow, mixture of water and ice/snow, dew, melting ice/snow, frost and icing rain) and “not known”.
1) The applied ensemble technique is based on covariance analyses of errors in ALADIN-CZ NWP forecasts. The technique is focused on variables with Gaussian error distributions, which is not the case for wind speed and total cloud cover; however, in our opinion, this is not an essential problem. The RST forecasts do not significantly depend on the wind speed, the errors in forecasts of cloud cover are quite symmetric, and the violation of normality is not very significant. 2) To perform the covariance analyses of errors in cloud cover forecasts, we used the SAF algorithm PGE01 Cloud Mask and found that it agrees well with subjective observations of cloud cover at synoptic stations. 3) Temperature, humidity and wind speed in the near-surface air layer are not sufficient to reasonably describe the uncertainties in RST forecasts. Radiation fluxes, which are derived from the total cloud cover, are crucial for RST development, and this quantity should be included in the ensemble generation. Even after doing so, however, the ensemble is too confident (the spread is not sufficient) and cannot describe the natural uncertainty in the forecasts. 4) There are several primary reasons for the insufficient ensemble spread and consequent underestimation of the forecast uncertainty. One reason is that NWP model errors in forecasted rain and snow are not considered when generating ensembles. It is clear that the road surface conditions (dry, wet, or covered by snow) significantly influence the development of the RST. However, we did not include these
characteristics in the ensembles owing to their distinct nonGaussian character of these characteristics. Other important reasons for errors are of a technical nature, for example, incorrect information on the sky view factor, unavailable or inaccurate information on the current road surface conditions, and inaccurate information on the technological parameters of roads. We can expect that these data shortcomings will at least be partly resolved in the future. 5) Although the ensemble technique underestimates the uncertainties in RST forecasts, it can be used in decision-making by road authorities who provide winter road maintenance. The outputs can be subjectively “corrected”, or a technique such as Model Output Statistics can be applied to objectively correct the forecasts. 6) In our opinion, the ensemble approach is especially suitable for predicting road surface conditions because it can express the probability of the occurrence of various road surface states, which deterministic forecasts do not allow. Unfortunately, we were unable to objectively verify the ensemble forecasts because we currently do not have required observation data. Preliminary results are promising. Acknowledgements This work was supported by the Czech Science Foundation grant GACR 13-34856S and TA01031509. The data were kindly provided by the Czech Hydrometeorological Institute. We acknowledge comments of reviewers which significantly improved the manuscript.
Z. Sokol et al. / Atmospheric Research 187 (2017) 33–41
Appendix A. Generation of a Gaussian random vector u having a known covariance matrix C = cov(u, uT) A Gaussian random vector u that has a known covariance matrix C = cov(u, uT) can be calculated by u ¼ Lz;
ðA1Þ
where L is a lower triangle matrix obtained by the Cholesky decomposition of C, i.e., C ¼ LLT ;
ðA2Þ
and z = (z1, z2, …, zn)T is an independent Gaussian random vector with zero mean and unit variance (Huynh et al., 2008). In our application, we first constructed the matrix C defined by Eqs. (2) and (3) using historical data T (Eq. (1)) and calculated its Cholesky decomposition matrix L (Eq. (A2)). Second, we generated a Gaussian random number with zero mean and unit random z⁎ by a standard algorithm (e.g., Huynh et al., 2008; Wilks, 2006). Finally, we calculated T⁎ = Lz⁎ and obtained a forecast for one member of the ensemble by integrating the METRo-CZ model using data from vector T⁎ rather than the forecasted NWP model data. By repeatedly generating random z⁎ values, we obtained the whole ensemble of forecasts. In our experiments, the first ensemble member corresponded to the original NWP model data. Appendix B. Total cloud cover estimated from satellite data Measurements by Meteosat Second Generation (MSG) satellites were used to estimate the total cloud cover (N). MSG data are available every 15 min, and their horizontal resolution is approximately 5 × 5 km over Central Europe. To calculate N, the SAF algorithm PGE01 Cloud Mask is applied (SAFNWC, 2013; Badescu et al., 2016). It classifies cloud cover into the following 3 classes for each satellite grid point: cloud-free, cloud-contaminated (partly cloudy or semitransparent) and cloud-filled (opaque clouds) to which we assign 0, 4 and 8 octas of N, respectively. The obtained N is interpolated from the satellite grid into the new grid using the nearest neighbour technique with a horizontal resolution of 1 km. Finally, N is calculated for each new grid G as the mean value of the grid values of N from the surrounding
Fig. B1. A histogram of the differences between N observed by observers at synoptic stations and calculated from the MSG data using m = 7 and m = 9.
41
m × m grids, where m is an odd number, and grid G is in the centre of the m × m grids. This procedure simulates how N is subjectively estimated by observers. Based on tests, we used m = 9 in this study. The results of the algorithm were compared with N from SYNOP station reports, which contain the subjectively determined N. The comparison included data from 17 synoptic stations from the Czech Republic and from 4 winter months for the years 2012–2013. A histogram of the differences is shown in Fig. B1 and indicates that m = 9 and m = 7 yielded very similar results. Taking into account the uncertainty in the subjectively determined N, Fig. B1 indicates a good agreement and acceptable asymmetry between the two estimations. Similar column heights for values of x = 2, 3 and 4 that do not correspond to Gaussian distributions are caused by cases when the subjective observations reported N = 7 and N = 8, but the applied algorithm recognized semitransparent cloudiness, which was classified as N = 4. References Badescu, V., Paulescu, M., Brabec, M., 2016. Reconstruction of effective cloud fields geometry from series of Sunshine Number. Atmos. Res. (accepted). Bouallegue, Z.B., Theis, S.E., Gebhardt, C., 2013. Enhancing COSMO-DE ensemble forecasts by inexpensive techniques. Meteorol. Z. 22, 49–59. Bouilloud, L., Martin, E., Habets, F., Boone, A., Moigne, P.L., Livet, J., Marchetti, M., Foidart, A., Franchistéguy, P., Morel, S., Noilhan, J., Pettre, P., 2009. Road surface condition forecasting in France. J. Appl. Meteorol. Climatol. 48, 2513–2527. Crevier, L.-P., Delage, Y., 2001. METRo. A new model for road-condition forecasting in Canada. J. Appl. Meteorol. 40, 2026–2037. Greenfield, T.M., Takle, E.S., 2006. Bridge frost prediction by heat and mass transfer methods. J. Appl. Meteorol. Climatol. 45, 517–526. Huynh, H.T., Lai, V.S., Soumaré, I., 2008. Stochastic Simulation and Applications in Finance With MATLAB Programs. John Wiley & Sons Ltd, p. 338 (ISBN 978–0–470-72538-2). Jansson, C., Almkvist, E., Jansson, P.-E., 2006. Heat balance of an asphalt surface: observations and physically-based simulations. Meteorol. Appl. 13, 203–212. Kangas, M., Heikinheimo, M., Hippi, M., 2015. RoadSurf: a modelling system for predicting road weather and road surface conditions. Meteorol. Appl. 22, 544–553. Leutbecher, M., Palmer, T.N., 2008. Ensemble forecasting. Weather Forecast. 227, 3515–3539. Linden, S.K., Petty, K.R., 2008. The Use of METRo (Model of the Environment and Temperature of the Roads) in Roadway Operation Decision Support Systems. 24th Conference on IIPS. NCAR, Boulder, CO. Rutz, J.J., Gibson, C.V., 2013. Integration of a road surface model into NWS operations, B. Am. Meteor. Society 94. SAFNWC, 2013. Algorithm Theoretical Basis Document for “Cloud Products” (CMa PGE01 v3.2, CT-PGE02 v2.2 & CTTH-PGE03 v2.2). Shao, J., Lister, P.J., 1996. An automated nowcasting model of road surface temperature and state for winter road maintenance. J. Appl. Meteorol. 35, 1352–1361. Sokol, Z., 2003. MOS based precipitation forecasts for river basins. Wea. Forecast 18, 769–781. Sokol, Z., Zacharov, P., Sedlák, P., Hošek, J., Bližňák, V., Chládová, Z., Pešice, P., Škuthan, M., 2014. First experience with the application of the METRo model in the Czech Republic. Atmos. Res. 143, 1–16. Wilks, D.S., 2006. Statistical Methods in the Atmospheric Sciences. second ed. Academic Press, p. 627 (ISBN : 9780080456225). Zakšek, K., Oštir, K., Kokalj, Ž., 2011. Sky-view factor as a relief visualization technique. Remote Sens. 3, 398–415.