**Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model**

Martina Tudor, Stjepan Ivatek-Šahdan, Antiono Stanešić, Kristian Horvath and Alica Bajić

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/55698

## **1. Introduction**

Numerical weather prediction (NWP) models are one of the factors contributing to the complex process of providing an accurate weather forecast. A number of global NWP models forecast weather over the whole Earth seven or more days in advance. These large scale global models do not provide high-resolution details that can be important for the weather we feel in a particular point (sensible weather). These details are often provided by limited area models (LAMs) that cover a particular area of interest as will be exemplified in this chapter. The Croatian Meteorological and Hydrological Service (CMHS) uses ALADIN (Aire Limitée Adaptation Dynamique développement InterNational, ALADIN International Team, 1997) limited area model for the operational weather forecast. This chapter describes the operational NWP aspects in CMHS and provides an insight into high impact weather phenomena that are typical for this region.

The quality of NWP model forecast depends on the NWP model that is used to compute the meteorological forecast fields. Model forecast should improve with the improvements in the model equations and parameterizations used to describe the atmospheric processes (Pielke 2002, Durran 1999, Werner et. al 1997). Another factor is the quality and availability of the measured atmospheric data used as input for the data assimilation procedure that creates the initial conditions for the model forecast. Limited area models also require the lateral boundary conditions data which are another factor contributing to the quality of the weather forecast. Both initial and lateral boundary conditions are usually taken from a larger scale and lower resolution NWP model covering larger geographical area. The initial conditions can be modified using a local initialization procedure.

© 2013 Tudor et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Tudor et al.; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

The operational forecast suite in CMHS uses the initial lateral boundary conditions from the ARPEGE (Action de Recherche Petite Echelle Grande Echelle, http://www.cnrm.meteo.fr/ gmapdoc/meshtml/guide\_ARP/arpege.html) global NWP model run operationally at Meteo-France. Another set of initial and boundary conditions available in CMHS is from the Inte‐ grated Forecast System (IFS, http://www.ecmwf.int/research/ifsdocs/CY36r1/ index.html) of the European Centre for Medium-Range Weather Forecast (ECMWF).

**4.** the NWP model forecast run that consists of:

**5.** post-processing of the model output data.

sets of input data, as will be described later.

**Figure 2.** The schedule and timing of the 00 UTC operational run in August 2011.

A considerable effort in the research and development of the ALADIN model has been devoted to the model optimization and various model options with an aim to provide the best possible forecast in the short time using available computer power (Yessad, 2011). To improve the 10

**b.** the NWP model forecast,

of atmospheric initial conditions for the model forecast run),

**a.** the data assimilation procedure that produces the analysis ("best possible" approximation

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

61

**c.** diagnostics of the derived parameters from the model variables during the forecast run,

Weather forecasts are time-critical applications. This means that the value of the products quickly degrades with time. In other words, the model forecast run has to finish before a certain time in the day, according to the requirements of the users. The choices of the input data and the model complexity that are used operationally also depend on the hour at which the weather forecast products should be available to the users. It is desired to keep the time span between the data measurements and the availability of the weather forecast as short as possible or within a reasonable interval usually determined by procedures in the forecast office or requirements of the other users. The time required for each forecast step in Figure 1 is determined by technical limitations such as the capacity of data transfer to the institution providing the weather forecast as well as the speed of the mainframe computer executing the NWP model software. The amount of data used in points (1) and (2) is limited by the capacity of data storages and speed of the data transfer lines. The complexity of the operational NWP model, as well as the size and grid resolution of the model domain is determined by this time span and the computer power. ARPEGE model data was the primary choice for the ALADIN operational model input in CMHS since it is available much earlier than the IFS model data (Fig 2). IFS is available later than ARPEGE because the data assimilation procedure starts later since it waits to assimilate more measured data. On the other hand, ARPEGE provides two

The files received from Meteo-France and ECMWF do not contain the global model fields covering the whole Earth on the native ARPEGE or IFS grid and model levels. To optimize the data transfer, the model data is interpolated to a limited-area Lambert projection grid in a resolution similar to the one of the global model in that area and it is covering a wider geographical area than the local LAM that uses it. Since ALADIN is a spectral model, initial and boundary conditions contain the meteorological data covering the whole model domain in the form of spectral coefficients, and not only the data on the lateral boundaries.

**Figure 1.** The scheme of the steps taken when making the weather forecast

The number of levels in the vertical is reduced (from 70 in ARPEGE or 91 in IFS to 37 in ALADIN Croatia), since LAM does not have to reach as high in the atmosphere as the global models. The global models have considerably longer forecast range (10 to 15 days in advance) that require knowledge of the atmospheric processes as high as the stratosphere top and mesosphere. These processes have negligible impact in the 72 hour forecast range used for the operational ALADIN forecast in CMHS..

The forecast process (Fig 1) consists of several steps:


The operational forecast suite in CMHS uses the initial lateral boundary conditions from the ARPEGE (Action de Recherche Petite Echelle Grande Echelle, http://www.cnrm.meteo.fr/ gmapdoc/meshtml/guide\_ARP/arpege.html) global NWP model run operationally at Meteo-France. Another set of initial and boundary conditions available in CMHS is from the Inte‐ grated Forecast System (IFS, http://www.ecmwf.int/research/ifsdocs/CY36r1/ index.html) of

The files received from Meteo-France and ECMWF do not contain the global model fields covering the whole Earth on the native ARPEGE or IFS grid and model levels. To optimize the data transfer, the model data is interpolated to a limited-area Lambert projection grid in a resolution similar to the one of the global model in that area and it is covering a wider geographical area than the local LAM that uses it. Since ALADIN is a spectral model, initial and boundary conditions contain the meteorological data covering the whole model domain

The number of levels in the vertical is reduced (from 70 in ARPEGE or 91 in IFS to 37 in ALADIN Croatia), since LAM does not have to reach as high in the atmosphere as the global models. The global models have considerably longer forecast range (10 to 15 days in advance) that require knowledge of the atmospheric processes as high as the stratosphere top and mesosphere. These processes have negligible impact in the 72 hour forecast range used for the

**1.** collect the measured data representative for the initial forecast time from upper-air soundings, surface stations, meteorological satellites, aircraft and radar data,

**2.** data quality control which removes measurements if values are too far from the first guess, removal of observations according to empirical (black)list and surface observations with

**3.** establishment of the initial and lateral boundary conditions using data obtained from a

too large height difference between model orography and real station height.

in the form of spectral coefficients, and not only the data on the lateral boundaries.

the European Centre for Medium-Range Weather Forecast (ECMWF).

60 Climate Change and Regional/Local Responses

**Figure 1.** The scheme of the steps taken when making the weather forecast

operational ALADIN forecast in CMHS..

larger scale (global) model,

The forecast process (Fig 1) consists of several steps:


Weather forecasts are time-critical applications. This means that the value of the products quickly degrades with time. In other words, the model forecast run has to finish before a certain time in the day, according to the requirements of the users. The choices of the input data and the model complexity that are used operationally also depend on the hour at which the weather forecast products should be available to the users. It is desired to keep the time span between the data measurements and the availability of the weather forecast as short as possible or within a reasonable interval usually determined by procedures in the forecast office or requirements of the other users. The time required for each forecast step in Figure 1 is determined by technical limitations such as the capacity of data transfer to the institution providing the weather forecast as well as the speed of the mainframe computer executing the NWP model software. The amount of data used in points (1) and (2) is limited by the capacity of data storages and speed of the data transfer lines. The complexity of the operational NWP model, as well as the size and grid resolution of the model domain is determined by this time span and the computer power. ARPEGE model data was the primary choice for the ALADIN operational model input in CMHS since it is available much earlier than the IFS model data (Fig 2). IFS is available later than ARPEGE because the data assimilation procedure starts later since it waits to assimilate more measured data. On the other hand, ARPEGE provides two sets of input data, as will be described later.

**Figure 2.** The schedule and timing of the 00 UTC operational run in August 2011.

A considerable effort in the research and development of the ALADIN model has been devoted to the model optimization and various model options with an aim to provide the best possible forecast in the short time using available computer power (Yessad, 2011). To improve the 10 m wind forecast in severe wind situations, like bura (local name for bora) windstorms, a procedure that provides a high resolution forecast of the wind 10 m above ground has been established. The high-resolution dynamical adaptation (HRDA) of the wind field using a hydrostatic version of the ALADIN model has been introduced to the operational suite soon after the 8 km resolution operational forecast runs have started in 2000 (Ivatek-Šahdan & Tudor, 2004).

## **2. Operational model characteristics**

ALADIN has been developed by a group of scientists from 16 countries and shares a consid‐ erable part with the global models IFS and ARPEGE. Both global models are a result of a coordinated effort of ECMWF and Meteo-France. Another LAM developed on the basis of IFS global model is HIRLAM (High Resolution Limited Area Model, Unden et al 2002). These models share the same source code in many parts of the model. A considerable part of the dynamics and data assimilation is in common too. ALADIN can use the same vertical discre‐ tization, grid-point dynamics and the various physics options as the global model ARPEGE. This section describes the operational model characteristics, the domains and other options used in the operational suite.

> ARPEGE global model and interpolated to the ALADIN grid. Afterwards, these fields are balanced using digital filter initialization (DFI) procedure (Lynch & Huang, 1994). However, the removal of the high-frequency wave energy from the initial conditions also affects the fast meteorological waves. Termonia (2008) showed how DFI can significantly reduce the depth

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

63

The second approach has not been used in CMHS, but blending has been used extensively in the ALADIN operational forecast suite in the Czech Meteorological and Hydrological Institute (CHMI, Brožkova et. al., 2001) as well as several other ALADIN member services (Hdidou, 2006). The blending procedure makes use of the fact that ALADIN is a spectral model. The initial model fields are constructed from the long waves coming from the low resolution model

The last approach is to initialize ALADIN model using data assimilation and currently it is used in experimental configuration in CMHS. Data assimilation is a procedure wherein information coming from measurements is combined with some a priori estimate (usually a short range forecast often called first-guess or background field) and their associated obser‐ vation and background model errors, respecting the model dynamical balances, to get the analysis - a maximal likelihood or "as best as possible" approximation of the true state of the atmosphere at a given time. There are different methods of solving this analysis problem (Hólm 2008) and at CMHS two of them are used: the optimal interpolation (OI, Courtier 1999) and the variational method (3DVAR, Hollingsworth et al, 1998). Optimal interpolation is used for updating land surface fields while 3DVAR is used for analysis of upper air fields. It was shown by Mahfouf (1991) that 2m analysis increments of temperature and humidity computed with OI can be used to update the land surface variables. As only the increments of 2m temperature

of the eye of the storm and proposes a solution in the scale-selective DFI (SSDFI).

(ARPEGE) and short waves from the high resolution model (ALADIN).

**Figure 3.** The ALADIN model domains and terrain height used operationally in CMHS.

#### **2.1. Model domains used in the operational suite**

The operational ALADIN model forecast is run on a Lambert-projection domain with 8 km horizontal resolution (Fig 3). The model fields are subsequently going through a dynamical adaptation procedure (Ivatek-Šahdan & Tudor, 2004) that produces 2 km resolution forecast of 10 m wind speed and gusts for a smaller domain shown in Figure 3. The procedure adapts the wind field of the 8 km resolution forecast and uses hydrostatic set-up of the ALADIN model with turbulence parameterization only. As an addition to the operational forecast, a 2 km resolution 24 hour forecast was established recently, that uses non-hydrostatic (NH) dynamics and the full parameterization set, including radiation, microphysics and convection schemes. This forecast is run for the same small domain (Fig 3).

## **2.2. Initialization and data assimilation**

In the ALADIN community, the initial conditions for the model forecast can be obtained using different approaches:


The first approach is used operationally in CMHS since the beginnings of the operational ALADIN forecast suite in 2000. There the upper-air and surface fields are taken from the

**Figure 3.** The ALADIN model domains and terrain height used operationally in CMHS.

m wind forecast in severe wind situations, like bura (local name for bora) windstorms, a procedure that provides a high resolution forecast of the wind 10 m above ground has been established. The high-resolution dynamical adaptation (HRDA) of the wind field using a hydrostatic version of the ALADIN model has been introduced to the operational suite soon after the 8 km resolution operational forecast runs have started in 2000 (Ivatek-Šahdan &

ALADIN has been developed by a group of scientists from 16 countries and shares a consid‐ erable part with the global models IFS and ARPEGE. Both global models are a result of a coordinated effort of ECMWF and Meteo-France. Another LAM developed on the basis of IFS global model is HIRLAM (High Resolution Limited Area Model, Unden et al 2002). These models share the same source code in many parts of the model. A considerable part of the dynamics and data assimilation is in common too. ALADIN can use the same vertical discre‐ tization, grid-point dynamics and the various physics options as the global model ARPEGE. This section describes the operational model characteristics, the domains and other options

The operational ALADIN model forecast is run on a Lambert-projection domain with 8 km horizontal resolution (Fig 3). The model fields are subsequently going through a dynamical adaptation procedure (Ivatek-Šahdan & Tudor, 2004) that produces 2 km resolution forecast of 10 m wind speed and gusts for a smaller domain shown in Figure 3. The procedure adapts the wind field of the 8 km resolution forecast and uses hydrostatic set-up of the ALADIN model with turbulence parameterization only. As an addition to the operational forecast, a 2 km resolution 24 hour forecast was established recently, that uses non-hydrostatic (NH) dynamics and the full parameterization set, including radiation, microphysics and convection schemes.

In the ALADIN community, the initial conditions for the model forecast can be obtained using

**2.** blending the large scale information from low resolution global model fields with the high

The first approach is used operationally in CMHS since the beginnings of the operational ALADIN forecast suite in 2000. There the upper-air and surface fields are taken from the

Tudor, 2004).

**2. Operational model characteristics**

62 Climate Change and Regional/Local Responses

**2.1. Model domains used in the operational suite**

This forecast is run for the same small domain (Fig 3).

**1.** using the large scale model data as initial conditions,

resolution features from the previous LAM forecast,

**2.2. Initialization and data assimilation**

**3.** using the data-assimilation procedure.

different approaches:

used in the operational suite.

ARPEGE global model and interpolated to the ALADIN grid. Afterwards, these fields are balanced using digital filter initialization (DFI) procedure (Lynch & Huang, 1994). However, the removal of the high-frequency wave energy from the initial conditions also affects the fast meteorological waves. Termonia (2008) showed how DFI can significantly reduce the depth of the eye of the storm and proposes a solution in the scale-selective DFI (SSDFI).

The second approach has not been used in CMHS, but blending has been used extensively in the ALADIN operational forecast suite in the Czech Meteorological and Hydrological Institute (CHMI, Brožkova et. al., 2001) as well as several other ALADIN member services (Hdidou, 2006). The blending procedure makes use of the fact that ALADIN is a spectral model. The initial model fields are constructed from the long waves coming from the low resolution model (ARPEGE) and short waves from the high resolution model (ALADIN).

The last approach is to initialize ALADIN model using data assimilation and currently it is used in experimental configuration in CMHS. Data assimilation is a procedure wherein information coming from measurements is combined with some a priori estimate (usually a short range forecast often called first-guess or background field) and their associated obser‐ vation and background model errors, respecting the model dynamical balances, to get the analysis - a maximal likelihood or "as best as possible" approximation of the true state of the atmosphere at a given time. There are different methods of solving this analysis problem (Hólm 2008) and at CMHS two of them are used: the optimal interpolation (OI, Courtier 1999) and the variational method (3DVAR, Hollingsworth et al, 1998). Optimal interpolation is used for updating land surface fields while 3DVAR is used for analysis of upper air fields. It was shown by Mahfouf (1991) that 2m analysis increments of temperature and humidity computed with OI can be used to update the land surface variables. As only the increments of 2m temperature

calculate increments of 2m temperature and relative humidity. These increments in a more or less sophisticated way are propagated into the land surface equivalents and used for updating the land surface variables. In the next step, the background upper air fields are analysed using 3DVAR. The standard NMC background error statistics is calculated by the aforementioned procedure over a 100-day period, from 15 Feb – 25 May 2008. The background error matrix was not tuned a posteriori such as in some other NWP systems (Bölöni and Horvath, 2010). The 3DVAR procedure includes a quality control of available data (screening) and minimiza‐ tion procedure that as output produces an analysis which is used for initiating the 6h forecast. The assimilation cycle is repeated every 6 hours, at 00, 06, 12 and 18 UTC. It is also run with a time delay big enough to enable usage of long cut off ARPEGE coupling files (the ARPEGE model is run later and in assimilation all data is used) for boundary conditions of 6h ALADIN forecast. Also long cut off data is used in cycle (data available after the time period needed to collect all observations of interest). Production from the assimilation cycle is done following same steps as in the assimilation cycle but at its end, a 72h forecast is performed. However, to perform a quasi-operational 72 hours forecast, timing constrain does not allow to use long cut off data and long cut off ARPEGE files; thus short cut off data and ARPEGE files are used.

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

65

The NWP model characteristics are usually described in the terms of the model numerics, dynamics and physics. Model numerics refers to the model computational domain, coordinate system, model resolution and grid as well as the mathematical methods used to solve the

Model dynamics refers to the resolvable processes that are resolved by the model grid and described by the set of model equations for horizontal and vertical momentum and conserva‐ tion of mass and thermodynamic properties. These processes encompass advection, pressure

The model dynamics is computed using semi-implicit time integration scheme (Robert, 1982). An implicit treatment of the gravity wave equation is absolutely stable (Durran, 1999). The semi-implicit scheme treats implicitly only a linearized form of the adjustment terms in the shallow water equations. The method results in solving the Helmholtz equation in spectral

ALADIN is a shallow-water spectral limited-area model. It applies Fourier spectral represen‐ tation of the model variables. That allows the advantage of fast Fourier transforms (FFTs) in both directions. An elliptic truncation that limits the Fourier series (Machenhauer & Haugen, 1987) ensures an isotropic horizontal resolution. The number of grid points in each horizontal direction of the whole integration area (N) is chosen so that N>3M+1, where M is the truncation wave-number in the same direction. This representation ensures that the nonlinear terms of the model equations are computed without aliasing. The model fields are transferred from

The accumulation of energy at the shortest wavelengths, due to spectral blocking, is reduced by a common 4th order numerical diffusion at the end of the time step. Semi-Lagrangian

**2.3. Model dynamics**

space.

system of the prognostic differential equations.

gradient force and adiabatic changes of heat/temperature.

spectral to grid-point space and back in each model time step.

**Figure 4.** Schematic view of the data assimilation cycle.

and relative humidity are needed, and only the surface measurements are used, rather simple OI method was implemented for computing the 2m analysis increments. For the upper air fields more control variables are defined, more sources of data are used, and more sophisti‐ cated variational algorithm (3DVAR) is implemented. 3DVAR belongs to the group of intermittent variational data assimilation strategies, and whereas the true state of the atmos‐ phere is never known, the principal challenge in determining the analysis is to estimate the background (model) errors. Namely, knowledge of the background error statistics is used to derive the analysis increments and transfer this information in horizontal and vertical directions in the atmosphere. The background error statistics is in CMHS estimated by using the so-called standard NMC method. In this approach, background error statistics is derived from a set of differences between forecasts that are valid at the same time, but at different forecast horizons, typically 12 and 36 hours. In addition, the multivariate formulation of the background error statistics that uses vorticity, divergence, temperature and surface pressure, and humidity as control variables (Berre, 2000), allows the propagation of the analysis increments from one to all the other variables, and thus increases the effect of each single observation. Finally, the potential use of the ensemble and seasonal formulations of the background error statistics is currently being investigated.

Practical implementation of data assimilation setup in CMHS is illustrated in Figure 4 that shows the scheme of one assimilation cycle. Because of limited computer resources the assimilation cycle is run in quasi-operational mode i.e. observation data is taken at the operational time, but the assimilation procedure is performed with certain time delay. An assimilation cycle is a sequence of a 6 hour forecast and analysis that is run on a daily basis. In the first step, SST in the 6h forecast from previous cycle (background) is replaced with SST taken from global model ARPEGE analysis. This replacement is done because no local SST analysis is implemented. In the second step, the background land surface variables are updated. As there are almost no or very little surface observations available, 2m observations are used to retrieve information about the surface temperature and soil water content. Therefore, after the quality control of the 2m observations, an optimal interpolation is used to calculate increments of 2m temperature and relative humidity. These increments in a more or less sophisticated way are propagated into the land surface equivalents and used for updating the land surface variables. In the next step, the background upper air fields are analysed using 3DVAR. The standard NMC background error statistics is calculated by the aforementioned procedure over a 100-day period, from 15 Feb – 25 May 2008. The background error matrix was not tuned a posteriori such as in some other NWP systems (Bölöni and Horvath, 2010). The 3DVAR procedure includes a quality control of available data (screening) and minimiza‐ tion procedure that as output produces an analysis which is used for initiating the 6h forecast. The assimilation cycle is repeated every 6 hours, at 00, 06, 12 and 18 UTC. It is also run with a time delay big enough to enable usage of long cut off ARPEGE coupling files (the ARPEGE model is run later and in assimilation all data is used) for boundary conditions of 6h ALADIN forecast. Also long cut off data is used in cycle (data available after the time period needed to collect all observations of interest). Production from the assimilation cycle is done following same steps as in the assimilation cycle but at its end, a 72h forecast is performed. However, to perform a quasi-operational 72 hours forecast, timing constrain does not allow to use long cut off data and long cut off ARPEGE files; thus short cut off data and ARPEGE files are used.

## **2.3. Model dynamics**

and relative humidity are needed, and only the surface measurements are used, rather simple OI method was implemented for computing the 2m analysis increments. For the upper air fields more control variables are defined, more sources of data are used, and more sophisti‐ cated variational algorithm (3DVAR) is implemented. 3DVAR belongs to the group of intermittent variational data assimilation strategies, and whereas the true state of the atmos‐ phere is never known, the principal challenge in determining the analysis is to estimate the background (model) errors. Namely, knowledge of the background error statistics is used to derive the analysis increments and transfer this information in horizontal and vertical directions in the atmosphere. The background error statistics is in CMHS estimated by using the so-called standard NMC method. In this approach, background error statistics is derived from a set of differences between forecasts that are valid at the same time, but at different forecast horizons, typically 12 and 36 hours. In addition, the multivariate formulation of the background error statistics that uses vorticity, divergence, temperature and surface pressure, and humidity as control variables (Berre, 2000), allows the propagation of the analysis increments from one to all the other variables, and thus increases the effect of each single observation. Finally, the potential use of the ensemble and seasonal formulations of the

Practical implementation of data assimilation setup in CMHS is illustrated in Figure 4 that shows the scheme of one assimilation cycle. Because of limited computer resources the assimilation cycle is run in quasi-operational mode i.e. observation data is taken at the operational time, but the assimilation procedure is performed with certain time delay. An assimilation cycle is a sequence of a 6 hour forecast and analysis that is run on a daily basis. In the first step, SST in the 6h forecast from previous cycle (background) is replaced with SST taken from global model ARPEGE analysis. This replacement is done because no local SST analysis is implemented. In the second step, the background land surface variables are updated. As there are almost no or very little surface observations available, 2m observations are used to retrieve information about the surface temperature and soil water content. Therefore, after the quality control of the 2m observations, an optimal interpolation is used to

background error statistics is currently being investigated.

**Figure 4.** Schematic view of the data assimilation cycle.

64 Climate Change and Regional/Local Responses

The NWP model characteristics are usually described in the terms of the model numerics, dynamics and physics. Model numerics refers to the model computational domain, coordinate system, model resolution and grid as well as the mathematical methods used to solve the system of the prognostic differential equations.

Model dynamics refers to the resolvable processes that are resolved by the model grid and described by the set of model equations for horizontal and vertical momentum and conserva‐ tion of mass and thermodynamic properties. These processes encompass advection, pressure gradient force and adiabatic changes of heat/temperature.

The model dynamics is computed using semi-implicit time integration scheme (Robert, 1982). An implicit treatment of the gravity wave equation is absolutely stable (Durran, 1999). The semi-implicit scheme treats implicitly only a linearized form of the adjustment terms in the shallow water equations. The method results in solving the Helmholtz equation in spectral space.

ALADIN is a shallow-water spectral limited-area model. It applies Fourier spectral represen‐ tation of the model variables. That allows the advantage of fast Fourier transforms (FFTs) in both directions. An elliptic truncation that limits the Fourier series (Machenhauer & Haugen, 1987) ensures an isotropic horizontal resolution. The number of grid points in each horizontal direction of the whole integration area (N) is chosen so that N>3M+1, where M is the truncation wave-number in the same direction. This representation ensures that the nonlinear terms of the model equations are computed without aliasing. The model fields are transferred from spectral to grid-point space and back in each model time step.

The accumulation of energy at the shortest wavelengths, due to spectral blocking, is reduced by a common 4th order numerical diffusion at the end of the time step. Semi-Lagrangian interpolators can be more or less diffusive (Staniforth and Cote, 1991). Combining two interpolators of different diffusivity with the flow deformation as a weighting factor yields more physical horizontal diffusion scheme that is based on the physical properties of the flow (Vana et al., 2008). This semi-Lagrangian horizontal diffusion (SLHD) is combined with numerical diffusion that removes short waves from the high layers of the atmosphere.

model data in the coupling zone is interpolated linearly in time. The time span between the available LBC data of 3 hours and the coupling zone width is 8 grid-points in 8 km resolution. Three hours can be sufficiently long and 64 km can be narrow so that a particular meteoro‐ logical feature can cross the coupling zone within the coupling interval (Fig 5). Solutions to this problem have been proposed, that include an alternative coupling scheme or interpolation in time (Tudor & Termonia, 2010), but have not been implemented to the operational model so far. The operational 2km resolution forecast uses LBC data from the operational 8 km

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

67

In the two-way coupling approach, the flow information about the state of the atmosphere goes from the small-scale model to the large scale one, as well as in the opposite direction. This method requires simultaneous integration of both high and low resolution runs. For that reason, the operational 8km and 2 km resolution runs use one-way coupling that allows transfer of information only from the large-scale to the small-scale model run. The high

Model physics describe the processes unresolved by the model grid, radiation, cloud micro‐ physics and surface processes. The processes described by the model physics are said to be parameterized (Pielke, 2002). The sub-grid-scale and diabatic processes are parameterized so their contribution to the changes in the large-scale state of the atmosphere can be considered.

The convective processes redistribute momentum, heat and moisture in the vertical. The deep convection parameterization (Fig 6) used in ALADIN (Gerard & Geleyn, 2005) is a prognostic mass-flux scheme (Gerard, 2007) where convective processes are treated with the use of prognostic variables for updraft and downdraft vertical velocities and mesh fractions (Gerard

resolution forecast starts only after the low resolution one is finished.

resolution forecast with a 1 h interval.

**Figure 5.** Schematic view illustrating the coupling problem.

**2.5. Model physics**

et al., 2009).

The advection of the prognostic variables in the model is computed using two-time-level semi-Lagrangian scheme. The method takes the model grid-points as the arrival points of the trajectory. The trajectories are computed one time step backwards to the origin points.

The model prognostic variables that are involved in the semi-implicit computations in the hydrostatic version of the model are surface pressure, the horizontal wind components, temperature and water vapour. The non-hydrostatic dynamics involves two additional model variables: pressure departure and vertical divergence that are treated by the semi-implicit computations. The developments in physics have introduced more prognostic variables to the model, such as cloud water and ice, rain and snow, as well as convective updraft and down‐ draft vertical velocities and mesh fractions. These quantities can be (optionally) advected by the semi-Lagrangian scheme and diffused by SLHD but do not enter the semi-implicit computations. In the operational forecast run, all these variables are advected, but SLHD is applied only to water vapour, cloud water and ice. This configuration supports modelling the advection in the atmospheric front, but may spoil the forecast of rainfall due to orography or other local feature that does not move with the flow,

The finite difference method is used to solve the model equations in the vertical, on 37 levels of hybrid pressure type eta coordinate (Simmons & Burridge, 1981). The primitive prognostic equations are solved for the prognostic variables using the two time level, semi-implicit, semi-Lagrangian advection scheme with a second-order accurate treatment of the nonlinear residual (Gospodinov et al. 2001).

## **2.4. Lateral boundary conditions and coupling**

LAM uses the large scale model data in a narrow coupling zone on the lateral boundaries and at discrete time intervals. The coupling of the model variables is done using Davies (1976) relaxation scheme in a narrow zone on the lateral boundaries of the LAM domain (Fig 5). The model dynamics requires usage of time dependent and periodic LBCs (Haugen & Machenha‐ uer, 1993). The coupling procedure has to be applied at the very beginning or end of the gridpoint computations (Radnoti, 1995) due to constraints imposed by the model dynamics. Various schemes for the lateral boundary treatment are associated to different problems (Davies, 1983). Werner et al. (1997) give an overview of weaknesses of the LAM forecast caused by the LBCs.

The lateral boundary conditions (LBC) are operationally obtained with a 3 hour interval. This interval is a compromise between the need to reduce the amount of data that have to be transferred and stored and the need to capture fast cyclones, fronts and other meteorological phenomena that can enter the LAM domain through its lateral boundaries. The coupling scheme requires large scale data in the coupling zone for each LAM time step. The large scale model data in the coupling zone is interpolated linearly in time. The time span between the available LBC data of 3 hours and the coupling zone width is 8 grid-points in 8 km resolution. Three hours can be sufficiently long and 64 km can be narrow so that a particular meteoro‐ logical feature can cross the coupling zone within the coupling interval (Fig 5). Solutions to this problem have been proposed, that include an alternative coupling scheme or interpolation in time (Tudor & Termonia, 2010), but have not been implemented to the operational model so far. The operational 2km resolution forecast uses LBC data from the operational 8 km resolution forecast with a 1 h interval.

In the two-way coupling approach, the flow information about the state of the atmosphere goes from the small-scale model to the large scale one, as well as in the opposite direction. This method requires simultaneous integration of both high and low resolution runs. For that reason, the operational 8km and 2 km resolution runs use one-way coupling that allows transfer of information only from the large-scale to the small-scale model run. The high resolution forecast starts only after the low resolution one is finished.

**Figure 5.** Schematic view illustrating the coupling problem.

#### **2.5. Model physics**

interpolators can be more or less diffusive (Staniforth and Cote, 1991). Combining two interpolators of different diffusivity with the flow deformation as a weighting factor yields more physical horizontal diffusion scheme that is based on the physical properties of the flow (Vana et al., 2008). This semi-Lagrangian horizontal diffusion (SLHD) is combined with numerical diffusion that removes short waves from the high layers of the atmosphere.

The advection of the prognostic variables in the model is computed using two-time-level semi-Lagrangian scheme. The method takes the model grid-points as the arrival points of the trajectory. The trajectories are computed one time step backwards to the origin points.

The model prognostic variables that are involved in the semi-implicit computations in the hydrostatic version of the model are surface pressure, the horizontal wind components, temperature and water vapour. The non-hydrostatic dynamics involves two additional model variables: pressure departure and vertical divergence that are treated by the semi-implicit computations. The developments in physics have introduced more prognostic variables to the model, such as cloud water and ice, rain and snow, as well as convective updraft and down‐ draft vertical velocities and mesh fractions. These quantities can be (optionally) advected by the semi-Lagrangian scheme and diffused by SLHD but do not enter the semi-implicit computations. In the operational forecast run, all these variables are advected, but SLHD is applied only to water vapour, cloud water and ice. This configuration supports modelling the advection in the atmospheric front, but may spoil the forecast of rainfall due to orography or

The finite difference method is used to solve the model equations in the vertical, on 37 levels of hybrid pressure type eta coordinate (Simmons & Burridge, 1981). The primitive prognostic equations are solved for the prognostic variables using the two time level, semi-implicit, semi-Lagrangian advection scheme with a second-order accurate treatment of the nonlinear residual

LAM uses the large scale model data in a narrow coupling zone on the lateral boundaries and at discrete time intervals. The coupling of the model variables is done using Davies (1976) relaxation scheme in a narrow zone on the lateral boundaries of the LAM domain (Fig 5). The model dynamics requires usage of time dependent and periodic LBCs (Haugen & Machenha‐ uer, 1993). The coupling procedure has to be applied at the very beginning or end of the gridpoint computations (Radnoti, 1995) due to constraints imposed by the model dynamics. Various schemes for the lateral boundary treatment are associated to different problems (Davies, 1983). Werner et al. (1997) give an overview of weaknesses of the LAM forecast caused

The lateral boundary conditions (LBC) are operationally obtained with a 3 hour interval. This interval is a compromise between the need to reduce the amount of data that have to be transferred and stored and the need to capture fast cyclones, fronts and other meteorological phenomena that can enter the LAM domain through its lateral boundaries. The coupling scheme requires large scale data in the coupling zone for each LAM time step. The large scale

other local feature that does not move with the flow,

**2.4. Lateral boundary conditions and coupling**

(Gospodinov et al. 2001).

66 Climate Change and Regional/Local Responses

by the LBCs.

Model physics describe the processes unresolved by the model grid, radiation, cloud micro‐ physics and surface processes. The processes described by the model physics are said to be parameterized (Pielke, 2002). The sub-grid-scale and diabatic processes are parameterized so their contribution to the changes in the large-scale state of the atmosphere can be considered.

The convective processes redistribute momentum, heat and moisture in the vertical. The deep convection parameterization (Fig 6) used in ALADIN (Gerard & Geleyn, 2005) is a prognostic mass-flux scheme (Gerard, 2007) where convective processes are treated with the use of prognostic variables for updraft and downdraft vertical velocities and mesh fractions (Gerard et al., 2009).

Cloud microphysics describe the processes of condensation, evaporation, freezing and melting as well as the processes that transform the cloud water droplets and ice crystals into rain and snow (Fig 6). ALADIN uses a simple microphysics scheme with prognostic cloud water and ice, rain and snow (Catry et al., 2007) and a statistical approach for sedimentation of precipi‐ tation (Geleyn et al., 2008). The microphysics scheme of the resolved processes is kept as close as possible to the original scheme that was without the prognostic condensates. The original scheme is the Kessler (1969) type scheme with modifications (Geleyn et. al. 1994) that include melting and freezing and imposing a brutal transition from ice to water at the temperature triple point.

tic values of TKE (turbulent kinetic energy) according to Geleyn et al. (2006) and Redelsperg‐

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

69

The gravity wave drag and lift parameterization scheme describes the vertical momentum flux due to the atmospheric waves generated by the unresolved topographic features. The surface parameterization schemes describe the soil properties and their impact on the meteorological model variables. ISBA (Interaction Soil Biosphere Atmosphere) is the surface scheme used in the operational forecast (Noilhan & Planton, 1989) as well as in the surface data assimilation (Giard & Bazile, 2000). The scheme describes the exchange of heat and moisture between land surface and air using thermal and hydrological properties of particular soil types (described with dominant land-use type, percentages of clay, sand and silt, useful soil depth, thermal roughness length) and vegetation changes during the year (described by vegetation fraction, leaf-area index, surface resistance to evapotranspiration). The prognostic variables are temperature, liquid and solid water contents are computed on two layers representing the soil surface and deep soil properties, and additional variables describe surface snow reservoir, density and albedo and ice and water on leaves. Snow budget depends on snow precipitation,

Each parameterization scheme obtains the information on the state of the atmosphere from the model variables and possibly output from other parameterizations. Then it uses a set of closure assumptions that relate the parameterized process to the state of the atmosphere. The param‐ eterization schemes have the largest impact on the prediction of the sensible weather at the Earth's surface (Pielke, 2002). These schemes became more complex over time and interact with

Model diagnostics produces the numerical values of the sensible weather fields from the model variables, such as the accumulated precipitation, instantaneous cloudiness or maximum wind gusts. Wind, temperature and humidity are interpolated from the model levels to the standard meteorological measurement heights (10 and 2 meters above surface) using a parameterized

The Croatian mountains are relatively small in extend, many of them have width close to 10 km and length less than 50 km. But the mountain peaks reach over 1 km inland, and several mountain tops reach over 1.5 km very close to the coastline (Fig 3). These mountains are separated by deep and inhabited valleys while important roads and connections go through

Weather in Croatia is in many ways modified or even controlled by orography (Zaninović et. al. 2008), but the large scale models do not resolve many of these local weather patterns. These weather patterns result from interaction of synoptic forcing with orography. Improved representation of orography in a meteorological model is expected to improve the forecast of weather phenomena that are strongly influenced by local orography (Horvath et. al., 2011,

evaporation and melting as well as snow accumulated on vegetation.

each other, the numerical and dynamical parts of the model.

er et al. (2001).

profile (Geleyn, 1988).

the mountain passes.

Bajić, 2007, Branković et. al., 2008).

**2.6. High resolution dynamical adaptation**

The radiation processes described in the model encompass the transfer, scattering, absorption and reflection of the shortwave solar radiation and long-wave thermal radiation of the Earth's surface and clouds. There are several radiation schemes available in the ALADIN model (Morcrette, 1989, Mlawer et. al. 1997), but the simplest one is used in the operational version. The operational scheme (Ritter and Geleyn 1992) is based on Geleyn and Hollingsworth (1979) scheme. It is simple and computationally cheap since it uses only one spectral band for longwave and one for short-wave radiation computations. The scheme has been enhanced recently (Geleyn et. al. 2005a, 2005b) but these modifications did not improve the cloud and temperature forecast in the stratus and fog case (Tudor, 2010).

**Figure 6.** Schematic view of the scheme of the prognostic variables describing microphysics and convection and asso‐ ciated processes in the version of the ALADIN model used in CMHS.

The turbulence parametrization scheme describes the impact of the unresolved motion and surface roughness on the vertical transfer of heat, momentum and moisture. The turbulent exchange coefficients are modified from Louis et al. (1982) and are computed using prognos‐

tic values of TKE (turbulent kinetic energy) according to Geleyn et al. (2006) and Redelsperg‐ er et al. (2001).

The gravity wave drag and lift parameterization scheme describes the vertical momentum flux due to the atmospheric waves generated by the unresolved topographic features. The surface parameterization schemes describe the soil properties and their impact on the meteorological model variables. ISBA (Interaction Soil Biosphere Atmosphere) is the surface scheme used in the operational forecast (Noilhan & Planton, 1989) as well as in the surface data assimilation (Giard & Bazile, 2000). The scheme describes the exchange of heat and moisture between land surface and air using thermal and hydrological properties of particular soil types (described with dominant land-use type, percentages of clay, sand and silt, useful soil depth, thermal roughness length) and vegetation changes during the year (described by vegetation fraction, leaf-area index, surface resistance to evapotranspiration). The prognostic variables are temperature, liquid and solid water contents are computed on two layers representing the soil surface and deep soil properties, and additional variables describe surface snow reservoir, density and albedo and ice and water on leaves. Snow budget depends on snow precipitation, evaporation and melting as well as snow accumulated on vegetation.

Each parameterization scheme obtains the information on the state of the atmosphere from the model variables and possibly output from other parameterizations. Then it uses a set of closure assumptions that relate the parameterized process to the state of the atmosphere. The param‐ eterization schemes have the largest impact on the prediction of the sensible weather at the Earth's surface (Pielke, 2002). These schemes became more complex over time and interact with each other, the numerical and dynamical parts of the model.

Model diagnostics produces the numerical values of the sensible weather fields from the model variables, such as the accumulated precipitation, instantaneous cloudiness or maximum wind gusts. Wind, temperature and humidity are interpolated from the model levels to the standard meteorological measurement heights (10 and 2 meters above surface) using a parameterized profile (Geleyn, 1988).

## **2.6. High resolution dynamical adaptation**

Cloud microphysics describe the processes of condensation, evaporation, freezing and melting as well as the processes that transform the cloud water droplets and ice crystals into rain and snow (Fig 6). ALADIN uses a simple microphysics scheme with prognostic cloud water and ice, rain and snow (Catry et al., 2007) and a statistical approach for sedimentation of precipi‐ tation (Geleyn et al., 2008). The microphysics scheme of the resolved processes is kept as close as possible to the original scheme that was without the prognostic condensates. The original scheme is the Kessler (1969) type scheme with modifications (Geleyn et. al. 1994) that include melting and freezing and imposing a brutal transition from ice to water at the temperature

The radiation processes described in the model encompass the transfer, scattering, absorption and reflection of the shortwave solar radiation and long-wave thermal radiation of the Earth's surface and clouds. There are several radiation schemes available in the ALADIN model (Morcrette, 1989, Mlawer et. al. 1997), but the simplest one is used in the operational version. The operational scheme (Ritter and Geleyn 1992) is based on Geleyn and Hollingsworth (1979) scheme. It is simple and computationally cheap since it uses only one spectral band for longwave and one for short-wave radiation computations. The scheme has been enhanced recently (Geleyn et. al. 2005a, 2005b) but these modifications did not improve the cloud and temperature

**Figure 6.** Schematic view of the scheme of the prognostic variables describing microphysics and convection and asso‐

The turbulence parametrization scheme describes the impact of the unresolved motion and surface roughness on the vertical transfer of heat, momentum and moisture. The turbulent exchange coefficients are modified from Louis et al. (1982) and are computed using prognos‐

triple point.

68 Climate Change and Regional/Local Responses

forecast in the stratus and fog case (Tudor, 2010).

ciated processes in the version of the ALADIN model used in CMHS.

The Croatian mountains are relatively small in extend, many of them have width close to 10 km and length less than 50 km. But the mountain peaks reach over 1 km inland, and several mountain tops reach over 1.5 km very close to the coastline (Fig 3). These mountains are separated by deep and inhabited valleys while important roads and connections go through the mountain passes.

Weather in Croatia is in many ways modified or even controlled by orography (Zaninović et. al. 2008), but the large scale models do not resolve many of these local weather patterns. These weather patterns result from interaction of synoptic forcing with orography. Improved representation of orography in a meteorological model is expected to improve the forecast of weather phenomena that are strongly influenced by local orography (Horvath et. al., 2011, Bajić, 2007, Branković et. al., 2008).

**Figure 7.** m wind forecast with 8 km (a) and 2 km resolution (b) for 6 UTC 5th Feb 2003.

The dynamical adaptation method of Žagar & Rakovec (1999) has been adapted for the purpose of operational forecast of the 10 m wind. The method provided successful operational forecast of the 10 m wind (Ivatek-Šahdan & Tudor, 2004) an has been used extensively in research impact studies (Bajić et al., 2007) as well as case studies (Tudor & Ivatek-Šahdan, 2002) of severe wind.

The meteorological model fields are first interpolated from the low resolution (8km in this case) to a higher resolution (2km) grid, but on considerably lower number of model levels, from 37 to 15 levels in the vertical. The number of vertical levels is reduced to minimize the computational cost. The levels close to the ground are of similar density, but the levels higher in the atmosphere are mostly omitted. Then a hydrostatic version of the ALADIN model is run for 30 time steps with a 60 second time step. The same large scale model data is used for the initial and LBC data, therefore the fields on the lateral boundaries do not change during the adaptation procedure. Turbulence is the only parameterization scheme used. Contributions from the moist and radiation processes are not computed to accelerate the model run.

full complexity of the physical parameterization package. These conclusions were based on analysis of a number cases of bura associated with a small scale pressure disturbance. Once a pressure disturbance is modelled, the wind field acts according to the measurements. The disturbances have to be detected first in the pressure measurements from the automatic stations that are performed on a less dense network than the wind measurements. The operational forecast using the latter configuration has been established only recently. The

**Figure 8.** Vertical cross-section of wind speed (shaded), direction (vectors) and potential temperature (black isolines)

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

71

In CMHS, the operational ALADIN forecast is run twice per day, starting from 00 and 12 UTC analyses. The operational ALADIN model forecast is run 72 hours in advance on a Lambertprojection domain with 8 km horizontal resolution on 37 hybrid sigma-pressure levels in the vertical. The model fields are subsequently going through a dynamical adaptation procedure (Ivatek-Šahdan and Tudor, 2004) that produces 2 km resolution forecast of 10 m wind speed and gusts. There are two more sets of initial and LBC data available from both ARPEGE and IFS. The analyses at 06 and 18 UTC are used for cycling of the data-assimilation procedure. Only 6 hour forecasts are produced from these analyses and they are used as a first guess in the data assimilation procedure for the next analysis time (12 and 00 UTC) that initiates the 72 hour forecast. Alternative 8 km resolution 72 h forecasts run from the initial fields created in 00 and 12 UTC data-assimilation. These runs provide an alternative forecast fields that differs slightly from the first operational run. As an addition to the operational forecast, 2 km resolution 24 hour forecast was established recently, that uses the non-hydrostatic (NH)

number of cases is still too short to make a thorough statistical analysis.

**2.7. The operational forecast schedule**

for 6 UTC 5th Feb 2003.

HRDA improves 10 m wind forecast (by 15%), especially in weather situation with strong bura wind, since wind-speed and direction depend on the terrain configuration upstream. Figures 7a and 7b illustrates the impact of high resolution orography on the 10m wind forecast. HRDA wind is much weaker (reduces from 10-15 m/s to less than 5m/s) and changes the direction (up to 180°, depending on location) from the low resolution forecast in the valley upstream of the mountain. Downstream of the mountain, wind is much stronger (increases from 15 m/s in low resolution model to more than 30 m/s) in the high resolution model run as mountain wave breaks and windstorm reaches the mountain slope (Fig 8).

The 10 m wind forecast from HRDA has been used in several applications, and improved the safety of the traffic in the air, on the sea as well as on land. Unfortunately, there were also events of short bura episodes connected to a local pressure disturbance (Tudor & Ivatek-Šahdan, 2010) that were not predicted by HRDA. These were used in a further study to assess the impact of the neglected effects to these bura episodes. This procedure is able to forecast bura onset, duration and strength if it is a consequence of a synoptic forcing. It misses those cases initiated by small scale disturbances in the pressure fields. These disturbances have to be modelled first, using the non-hydrostatic model version on more levels in the vertical with Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model http://dx.doi.org/10.5772/55698 71

**Figure 8.** Vertical cross-section of wind speed (shaded), direction (vectors) and potential temperature (black isolines) for 6 UTC 5th Feb 2003.

full complexity of the physical parameterization package. These conclusions were based on analysis of a number cases of bura associated with a small scale pressure disturbance. Once a pressure disturbance is modelled, the wind field acts according to the measurements. The disturbances have to be detected first in the pressure measurements from the automatic stations that are performed on a less dense network than the wind measurements. The operational forecast using the latter configuration has been established only recently. The number of cases is still too short to make a thorough statistical analysis.

#### **2.7. The operational forecast schedule**

(a) (b)

70 Climate Change and Regional/Local Responses

**Figure 7.** m wind forecast with 8 km (a) and 2 km resolution (b) for 6 UTC 5th Feb 2003.

breaks and windstorm reaches the mountain slope (Fig 8).

wind.

The dynamical adaptation method of Žagar & Rakovec (1999) has been adapted for the purpose of operational forecast of the 10 m wind. The method provided successful operational forecast of the 10 m wind (Ivatek-Šahdan & Tudor, 2004) an has been used extensively in research impact studies (Bajić et al., 2007) as well as case studies (Tudor & Ivatek-Šahdan, 2002) of severe

The meteorological model fields are first interpolated from the low resolution (8km in this case) to a higher resolution (2km) grid, but on considerably lower number of model levels, from 37 to 15 levels in the vertical. The number of vertical levels is reduced to minimize the computational cost. The levels close to the ground are of similar density, but the levels higher in the atmosphere are mostly omitted. Then a hydrostatic version of the ALADIN model is run for 30 time steps with a 60 second time step. The same large scale model data is used for the initial and LBC data, therefore the fields on the lateral boundaries do not change during the adaptation procedure. Turbulence is the only parameterization scheme used. Contributions from the moist and radiation processes are not computed to accelerate the model run.

HRDA improves 10 m wind forecast (by 15%), especially in weather situation with strong bura wind, since wind-speed and direction depend on the terrain configuration upstream. Figures 7a and 7b illustrates the impact of high resolution orography on the 10m wind forecast. HRDA wind is much weaker (reduces from 10-15 m/s to less than 5m/s) and changes the direction (up to 180°, depending on location) from the low resolution forecast in the valley upstream of the mountain. Downstream of the mountain, wind is much stronger (increases from 15 m/s in low resolution model to more than 30 m/s) in the high resolution model run as mountain wave

The 10 m wind forecast from HRDA has been used in several applications, and improved the safety of the traffic in the air, on the sea as well as on land. Unfortunately, there were also events of short bura episodes connected to a local pressure disturbance (Tudor & Ivatek-Šahdan, 2010) that were not predicted by HRDA. These were used in a further study to assess the impact of the neglected effects to these bura episodes. This procedure is able to forecast bura onset, duration and strength if it is a consequence of a synoptic forcing. It misses those cases initiated by small scale disturbances in the pressure fields. These disturbances have to be modelled first, using the non-hydrostatic model version on more levels in the vertical with

In CMHS, the operational ALADIN forecast is run twice per day, starting from 00 and 12 UTC analyses. The operational ALADIN model forecast is run 72 hours in advance on a Lambertprojection domain with 8 km horizontal resolution on 37 hybrid sigma-pressure levels in the vertical. The model fields are subsequently going through a dynamical adaptation procedure (Ivatek-Šahdan and Tudor, 2004) that produces 2 km resolution forecast of 10 m wind speed and gusts. There are two more sets of initial and LBC data available from both ARPEGE and IFS. The analyses at 06 and 18 UTC are used for cycling of the data-assimilation procedure. Only 6 hour forecasts are produced from these analyses and they are used as a first guess in the data assimilation procedure for the next analysis time (12 and 00 UTC) that initiates the 72 hour forecast. Alternative 8 km resolution 72 h forecasts run from the initial fields created in 00 and 12 UTC data-assimilation. These runs provide an alternative forecast fields that differs slightly from the first operational run. As an addition to the operational forecast, 2 km resolution 24 hour forecast was established recently, that uses the non-hydrostatic (NH) dynamics in the ALADIN model and the ALADIN's full parameterization set, including the convection scheme. This forecast runs is performed once per day, following the 00 UTC operational 8 km resolution forecast. It uses the 6 hour forecast from the 8 km resolution operational run as input initial file and runs with SSDFI. This high-resolution forecast is run for 24 hours, until 6 UTC on the next day. This procedure allows covering the 24 hour period used to collect precipitation data from the rain-gauges.

## **2.8. Predictability**

High impact weather events are often forecasted by means of their predictability. ALADIN has been used in several predictability studies (e.g. Branković et al., 2007) as well as the Limited Area Ensemble Forecasting (Wang et al., 2011) system. In most of these studies, ECMWF ensemble forecasts were dynamically downscaled for several severe weather cases. Therefore, the initial conditions come from perturbations generated by singular vectors of ECMWF. Wang et al. (2011) combined different initial and boundary conditions from the perturbed global model with other aspects of forecast uncertainty, such as blending, multi-physics approach and breeding.

(a) (b)

**Figure 9.** m wind forecast with 8 km (a) and 2 km resolution (b), for 15 UTC 7th Nov 1999.

maximum of the wind component normal to the flight direction by 3 m/s.

(a) (b)

15-16 UTC 7th Nov 1999.

**Figure 10.** Comparison of measured and model normal wind components on flight levels 660 m (a) and 330 m, for

The strongest wind gust ever measured in Croatia was 69.0 m/s on Maslenica station (Bajić 2003) that was still operational during December 2003. During severe bura in that area, the 10 m wind speed is usually underestimated by the 8 km resolution forecast by 30 m/s on average

MAP Special Observing Period 15 (SOP 15), from 7th to 9th Nov 1999, was characterized by severe bura on the eastern Adriatic Sea (Fig 9). The bura was associated with an intensive cyclone that was slowly moving over central Italy. A cold air outbreak from northeast intensified the pressure gradient over the coastal mountains and intensified wind in the lower layers of the atmosphere. During SOP 15, special measurements were taken from ELECTRA aircraft, through a northwest-southeast transect, parallel to the eastern Adriatic coast (Grubi‐ šić, 2004). The measurements were taken on levels 330 m and 660 m above the sea level and in the period from 15 to 16 UTC on 7th Nov 1999 (Fig 10). The 8km resolution forecast underes‐ timated wind speed in the vicinity of the coastal mountains by 30% on average (Fig 9a). The wind forecast improved considerably (up to 10% over or under estimation of wind speed) using HRDA (Fig 9b). But the wind above the open sea (where the aircraft flew) was similar, both at surface (Fig 9) and at flight levels (Fig 10). The wind maximum at 660 m level is well predicted, as well as the northern minimum. The southern minimum is weaker by 5 m/s than measured and moved further south. On the lower level, at 330 m, model overestimates the

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

73

Predictability studies of bura cases have found little sensitivity to different initial conditions (Ivatek-Šahdan & Ivančan-Picek, 2006). Horvath et al. (2009) suggested a small effect of uncertainties in the upstream initial conditions on the bura events in the southern Adriatic Sea. Another study of severe bura (Branković et al., 2007) exemplified a case of gale force bura that was predicted with a probability exceeding 95%. Dynamical downscaling of ECMWF fields from 40 to 8 km resolution improves the precipitation rate and pattern in the intensive precipitation cases (Branković et al., 2008) too, the predicted precipitation is roughly double in the high resolution forecasts, but still reaches only 20% of the observed precipitation maximum.

## **3. High impact weather events**

The capability of the ALADIN model to predict high-impact weather events in Croatia is illustrated by several exemplary cases that encompass various weather types.

## **3.1. Bura (bora)**

Bura is a local word describing a downslope windstorm on the eastern Adriatic coast. It is a northeasterly wind with high gustiness that is mostly controlled by the upstream topography of Dinaric Alps. During a windstorm episode, several roads and ferry lines to the islands get closed. The strength and variability of bura wind varies from case to case (e.g. Horvath et al., 2009) and reaches up to 69 m/s and varies in space, although several places became famous for it. The processes that control bura occurrence and strength are briefly illustrated. Grisogono & Belušić (2009) give a comprehensive review of recent advances in bura research.

**Figure 9.** m wind forecast with 8 km (a) and 2 km resolution (b), for 15 UTC 7th Nov 1999.

dynamics in the ALADIN model and the ALADIN's full parameterization set, including the convection scheme. This forecast runs is performed once per day, following the 00 UTC operational 8 km resolution forecast. It uses the 6 hour forecast from the 8 km resolution operational run as input initial file and runs with SSDFI. This high-resolution forecast is run for 24 hours, until 6 UTC on the next day. This procedure allows covering the 24 hour period

High impact weather events are often forecasted by means of their predictability. ALADIN has been used in several predictability studies (e.g. Branković et al., 2007) as well as the Limited Area Ensemble Forecasting (Wang et al., 2011) system. In most of these studies, ECMWF ensemble forecasts were dynamically downscaled for several severe weather cases. Therefore, the initial conditions come from perturbations generated by singular vectors of ECMWF. Wang et al. (2011) combined different initial and boundary conditions from the perturbed global model with other aspects of forecast uncertainty, such as blending, multi-physics approach

Predictability studies of bura cases have found little sensitivity to different initial conditions (Ivatek-Šahdan & Ivančan-Picek, 2006). Horvath et al. (2009) suggested a small effect of uncertainties in the upstream initial conditions on the bura events in the southern Adriatic Sea. Another study of severe bura (Branković et al., 2007) exemplified a case of gale force bura that was predicted with a probability exceeding 95%. Dynamical downscaling of ECMWF fields from 40 to 8 km resolution improves the precipitation rate and pattern in the intensive precipitation cases (Branković et al., 2008) too, the predicted precipitation is roughly double in the high resolution forecasts, but still reaches only 20% of the observed precipitation

The capability of the ALADIN model to predict high-impact weather events in Croatia is

Bura is a local word describing a downslope windstorm on the eastern Adriatic coast. It is a northeasterly wind with high gustiness that is mostly controlled by the upstream topography of Dinaric Alps. During a windstorm episode, several roads and ferry lines to the islands get closed. The strength and variability of bura wind varies from case to case (e.g. Horvath et al., 2009) and reaches up to 69 m/s and varies in space, although several places became famous for it. The processes that control bura occurrence and strength are briefly illustrated. Grisogono

illustrated by several exemplary cases that encompass various weather types.

& Belušić (2009) give a comprehensive review of recent advances in bura research.

used to collect precipitation data from the rain-gauges.

**2.8. Predictability**

72 Climate Change and Regional/Local Responses

and breeding.

maximum.

**3.1. Bura (bora)**

**3. High impact weather events**

MAP Special Observing Period 15 (SOP 15), from 7th to 9th Nov 1999, was characterized by severe bura on the eastern Adriatic Sea (Fig 9). The bura was associated with an intensive cyclone that was slowly moving over central Italy. A cold air outbreak from northeast intensified the pressure gradient over the coastal mountains and intensified wind in the lower layers of the atmosphere. During SOP 15, special measurements were taken from ELECTRA aircraft, through a northwest-southeast transect, parallel to the eastern Adriatic coast (Grubi‐ šić, 2004). The measurements were taken on levels 330 m and 660 m above the sea level and in the period from 15 to 16 UTC on 7th Nov 1999 (Fig 10). The 8km resolution forecast underes‐ timated wind speed in the vicinity of the coastal mountains by 30% on average (Fig 9a). The wind forecast improved considerably (up to 10% over or under estimation of wind speed) using HRDA (Fig 9b). But the wind above the open sea (where the aircraft flew) was similar, both at surface (Fig 9) and at flight levels (Fig 10). The wind maximum at 660 m level is well predicted, as well as the northern minimum. The southern minimum is weaker by 5 m/s than measured and moved further south. On the lower level, at 330 m, model overestimates the maximum of the wind component normal to the flight direction by 3 m/s.

**Figure 10.** Comparison of measured and model normal wind components on flight levels 660 m (a) and 330 m, for 15-16 UTC 7th Nov 1999.

The strongest wind gust ever measured in Croatia was 69.0 m/s on Maslenica station (Bajić 2003) that was still operational during December 2003. During severe bura in that area, the 10 m wind speed is usually underestimated by the 8 km resolution forecast by 30 m/s on average (Fig 11a). The HRDA forecast improves the surface wind forecast as it puts the wind speed maximum closer to the surface (Fig 11b), especially to the slope of the mountain as the mountain wave breaks. The wind maximum is still underestimated by 10-12%.

**3.2. Forecast of the road conditions**

The bura variability in space and time has a pronounced influence on road traffic. Therefore, knowing bura characteristics is a necessary condition for road transport safety. To properly organize the traffic safety system, special emphasis should be given to the quality of measured long term wind speed and direction data and low resolution atmospheric forecast models.

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

75

**Figure 13.** The ANEMO-ALARM user interface. The measured (thick lines) and modelled (thin line) mean wind speed is

An application, named ANEMO-ALARM (Bajić et. al., 2008), has been developed that assists road authorities in managing the traffic on the roads affected by strong wind and turbulence. The application is based on measured and forecasted wind speed and gusts for a choice of locations on Croatian roads that are most affected by severe wind. The application communi‐ cates with the user through a graphical user interface (Fig 13). The interface shows current and expected alarm status for road traffic safety conditions for any of the three categories of vehicles (green is for open road, yellow for preparedness status and expected road closure and red

Prognosed cloud water and ice, rain and snow mixing ratios are valuable output of the operational model forecast. Their three-dimensional distribution helps in identifying the areas that are a potential threat to aviation. This data can provide important information on the state of the atmosphere inside a cloud, especially when in-situ and remote measurements are sparse or insufficient to describe more detailed structure of the cloud. This value is briefly illustrated by means of a case when weather contributed to the crash of an airplane on the Velebit

On 5th February 2009, a Cesna 303, crashed 500m north of the Vagan peak of the Velebit Mountain. The small airplane flew from Zagreb to Zadar, first at 2600 m, then at 2000 m feet.

shown on a graph as blue line, wind gusts are red and direction is green.

indicated that the road is closed).

**3.3. Air crash investigation**

Mountain.

**Figure 11.** Vertical cross-section of wind speed (shaded), direction (vectors) and potential temperature (black isolines) from 8 km (a) and 2km (b) runs, for 3 UTC 24th Dec 2003.

There were several severe wind episodes in December 2003 (Fig 12), particularly in the period from December 22 to 26 2003, during which infrastructure was damaged on the highway going through the area hit by the severe bura, including the Baričević meteorological station (Fig 12a). The 2 km resolution forecast is in better agreement to the measured data for strong bura events (Fig 12a) as well as the most severe one (Fig 12b). In the absence of atmospheric soundings, these results confirm that the vertical structure of the meteorological fields in the lower troposphere also improves in high resolution for this case.

**Figure 12.** Comparison of measured 10 minute maximum (light blue) and average (blue), 8 km resolution mean wind forecast (red), 2 km resolution mean wind (orange) and wind gust forecast (yellow) at Baričević (a) and Maslenica (b) for December 2003.

## **3.2. Forecast of the road conditions**

(Fig 11a). The HRDA forecast improves the surface wind forecast as it puts the wind speed maximum closer to the surface (Fig 11b), especially to the slope of the mountain as the

**Figure 11.** Vertical cross-section of wind speed (shaded), direction (vectors) and potential temperature (black isolines)

There were several severe wind episodes in December 2003 (Fig 12), particularly in the period from December 22 to 26 2003, during which infrastructure was damaged on the highway going through the area hit by the severe bura, including the Baričević meteorological station (Fig 12a). The 2 km resolution forecast is in better agreement to the measured data for strong bura events (Fig 12a) as well as the most severe one (Fig 12b). In the absence of atmospheric soundings, these results confirm that the vertical structure of the meteorological fields in the

**Figure 12.** Comparison of measured 10 minute maximum (light blue) and average (blue), 8 km resolution mean wind forecast (red), 2 km resolution mean wind (orange) and wind gust forecast (yellow) at Baričević (a) and Maslenica (b)

mountain wave breaks. The wind maximum is still underestimated by 10-12%.

(a) (b)

lower troposphere also improves in high resolution for this case.

(a) (b)

for December 2003.

from 8 km (a) and 2km (b) runs, for 3 UTC 24th Dec 2003.

74 Climate Change and Regional/Local Responses

The bura variability in space and time has a pronounced influence on road traffic. Therefore, knowing bura characteristics is a necessary condition for road transport safety. To properly organize the traffic safety system, special emphasis should be given to the quality of measured long term wind speed and direction data and low resolution atmospheric forecast models.

**Figure 13.** The ANEMO-ALARM user interface. The measured (thick lines) and modelled (thin line) mean wind speed is shown on a graph as blue line, wind gusts are red and direction is green.

An application, named ANEMO-ALARM (Bajić et. al., 2008), has been developed that assists road authorities in managing the traffic on the roads affected by strong wind and turbulence. The application is based on measured and forecasted wind speed and gusts for a choice of locations on Croatian roads that are most affected by severe wind. The application communi‐ cates with the user through a graphical user interface (Fig 13). The interface shows current and expected alarm status for road traffic safety conditions for any of the three categories of vehicles (green is for open road, yellow for preparedness status and expected road closure and red indicated that the road is closed).

#### **3.3. Air crash investigation**

Prognosed cloud water and ice, rain and snow mixing ratios are valuable output of the operational model forecast. Their three-dimensional distribution helps in identifying the areas that are a potential threat to aviation. This data can provide important information on the state of the atmosphere inside a cloud, especially when in-situ and remote measurements are sparse or insufficient to describe more detailed structure of the cloud. This value is briefly illustrated by means of a case when weather contributed to the crash of an airplane on the Velebit Mountain.

On 5th February 2009, a Cesna 303, crashed 500m north of the Vagan peak of the Velebit Mountain. The small airplane flew from Zagreb to Zadar, first at 2600 m, then at 2000 m feet. Before the crash, the airplane descended from 1900 to 1600 m. The time of the accident was recorded as 13:53 UTC when air traffic control lost contact with the pilot.

Strong southwest wind prevailed throughout the troposphere. Wind speed increased with height, upstream of the mountain, in Zadar, the maximum wind speed of 20 knots was reached at 750m reaching 45 knots just below 3 km in the Zagreb radiosonde measurements. Temper‐ atures on Zadar airport (88 m) and in Zagreb (128 m) were 14 to 15 °C. On Zavižan peak (1594 m) of the Velebit mountain temperature was 0.3°C, relative humidity was 100% and there was very low visibility, indicating that the mountain top was in a cloud. The satellite images show Velebit mountain covered by a thick cloud and mountain waves in the lee (Fig 14a). The cloudtop temperature was below 0°C (Fig 14b).

**Figure 14.** Satellite picture of clouds (a) and cloud top temperature (b) for 14 UTC 5th Feb 2009.

The model vertical cross-section of vertical velocity and potential temperature (Fig 15a) shows strong downward motion above and north of Vagan peak and variable vertical velocity downstream. Changes in the vertical velocity are accompanied by waves in the isolines of potential temperature (blue lines in Fig 15a). As the airplane flew from Zagreb to Zadar, it faced strong opposing wind and changing upward and downward motions that were stronger as the airplane approached the Vagan peak.

The mountain was covered by layers of clouds consisting of both cloud water and ice (Fig 15b) and high relative humidity existed in the cloud free area. The 0o C-isoline was above the 6500 feet flight level over the valley close to Zagreb, but the 0°C isoline descends below that height above the mountains close to Vagan peak (Fig 15b). Turbulent kinetic energy had high values above and downstream of the mountain (more than 2.5\*104 m2 /s2 ).

#### **3.4. Jugo (sirocco)**

The Adriatic region is often affected by another severe wind that is rather familiar to the local population. Jugo is a local word describing the southeasterly wind in the eastern Adriatic. It interupts the ferry lines to the islands as well as other activities on the sea and the coast since it generates large sea waves due to the long fetch. It is often associated with heavy rainfall. The strength of jugo wind is often underestimated in southern Adriatic, especially in Dubrovnik and on Palagruža, Lastovo and Vis islands by 15% (2 km resolution) to 30% (8km resolution) on average. The time-variability of jugo wind speed and direction is underestimated over the open sea (at least in some cases), as was revealed in the comparison with the measurements in the open sea during the DART (Dynamics of the Adriatic in Real-Time) oceanic research cruise in March 2006 (Tudor, 2011). The model error of the wind speed forecast is not uniform

**Figure 15.** Vertical cross-sections of (a) vertical velocity (shaded), wind vectors parallel to the cross-section (vectors), potential temperature (K, blue lines) and (b) cloud water and ice (shaded), relative humidity (blue lines) and tempera‐

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

77

and can exceed 50% in some intervals for both 8 and 2 km resolution forecast.

ture (°C, red lines) for 14 UTC February 5, 2009. The flight levels are shown as horizontal green lines.

(a)

(b)

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model http://dx.doi.org/10.5772/55698 77

Before the crash, the airplane descended from 1900 to 1600 m. The time of the accident was

Strong southwest wind prevailed throughout the troposphere. Wind speed increased with height, upstream of the mountain, in Zadar, the maximum wind speed of 20 knots was reached at 750m reaching 45 knots just below 3 km in the Zagreb radiosonde measurements. Temper‐ atures on Zadar airport (88 m) and in Zagreb (128 m) were 14 to 15 °C. On Zavižan peak (1594 m) of the Velebit mountain temperature was 0.3°C, relative humidity was 100% and there was very low visibility, indicating that the mountain top was in a cloud. The satellite images show Velebit mountain covered by a thick cloud and mountain waves in the lee (Fig 14a). The cloud-

recorded as 13:53 UTC when air traffic control lost contact with the pilot.

(a) (b)

and high relative humidity existed in the cloud free area. The 0o

above and downstream of the mountain (more than 2.5\*104

**Figure 14.** Satellite picture of clouds (a) and cloud top temperature (b) for 14 UTC 5th Feb 2009.

The model vertical cross-section of vertical velocity and potential temperature (Fig 15a) shows strong downward motion above and north of Vagan peak and variable vertical velocity downstream. Changes in the vertical velocity are accompanied by waves in the isolines of potential temperature (blue lines in Fig 15a). As the airplane flew from Zagreb to Zadar, it faced strong opposing wind and changing upward and downward motions that were stronger

The mountain was covered by layers of clouds consisting of both cloud water and ice (Fig 15b)

feet flight level over the valley close to Zagreb, but the 0°C isoline descends below that height above the mountains close to Vagan peak (Fig 15b). Turbulent kinetic energy had high values

The Adriatic region is often affected by another severe wind that is rather familiar to the local population. Jugo is a local word describing the southeasterly wind in the eastern Adriatic. It interupts the ferry lines to the islands as well as other activities on the sea and the coast since it generates large sea waves due to the long fetch. It is often associated with heavy rainfall. The strength of jugo wind is often underestimated in southern Adriatic, especially in Dubrovnik

 m2 /s2 ).

C-isoline was above the 6500

top temperature was below 0°C (Fig 14b).

76 Climate Change and Regional/Local Responses

as the airplane approached the Vagan peak.

**3.4. Jugo (sirocco)**

**Figure 15.** Vertical cross-sections of (a) vertical velocity (shaded), wind vectors parallel to the cross-section (vectors), potential temperature (K, blue lines) and (b) cloud water and ice (shaded), relative humidity (blue lines) and tempera‐ ture (°C, red lines) for 14 UTC February 5, 2009. The flight levels are shown as horizontal green lines.

and on Palagruža, Lastovo and Vis islands by 15% (2 km resolution) to 30% (8km resolution) on average. The time-variability of jugo wind speed and direction is underestimated over the open sea (at least in some cases), as was revealed in the comparison with the measurements in the open sea during the DART (Dynamics of the Adriatic in Real-Time) oceanic research cruise in March 2006 (Tudor, 2011). The model error of the wind speed forecast is not uniform and can exceed 50% in some intervals for both 8 and 2 km resolution forecast.

or without 3Dvar and using ARPEGE and IFS initial and boundary conditions) One run

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

79

**Figure 17.** Accumulated precipitation from the operational (a) and data assimilation (b) forecast for 25th Sep 2010.

The flash flood hit Dubrovnik, Croatia, on November 22, 2010. The rain gauge measurements exceeded 100 mm/24h in the area and the one in Dubrovnik measured 161.4 mm/24h (the highest on record by that time), with a peak intensity of 71.5 mm/h (Fig 18b). The flood water

The operational and parallel suite model results forecast large 24 hour accumulated rainfall amounts in the area around Dubrovnik, with the maxima over 100 mm located above the surrounding areas of Montenegro and Bosnia and Herzegovina. The results of the highresolution (2 km) non-hydrostatic ALADIN model run showed a dependency of the forecasts on the input data from initial and lateral boundary conditions. The position and time of the precipitation maxima in the high resolution ALADIN output fields (Fig 18a) was similar to the lower resolution (8 km) run used for initial and lateral boundary conditions. The peak intensity of precipitation in all model runs was late (3 hours in runs coupled to Arpege, 6 hours for IFS). The reason for this delay is yet unknown and could be attributed to the absence of measure‐

**Figure 18.** Accumulated 24 h precipitation from the 2 km resolution run (a) and measured (b) for 23rd Nov 2010.

generated additional intensive rain band over an island 25 km east of Pula.

(a) (b)

level reached 1.5 m caused damage and endangered lives.

(a) (b)

ments over the southern Adriatic.

**Figure 16.** Comparison of measured (purple), 2 km (full lines) and 8 km resolution wind speed forecast (dashed) for Palagruža (a), Dubrovnik (b), Hvar (c) and Mljet (d) for 7th-9th Nov 2010.

#### **3.5. Heavy rain**

Large amounts of rain (> 100 mm) that fall during a short time period (< 6 hours) can cause flash floods, especially in places where the terrain supports accumulation of water as on low slope coast below a mountain. The ability of the operational ALADIN model to predict such cases is exemplarily illustrated by two cases.

On September 25, 2010, just after midnight an intensive rain hit Pula city on the southern part of the Istria Peninsula, Croatia. The rain was intensive for several hours and the rainfall rate measured by the ombrograph reached 43.9 mm per hour. Several rain-gauges in the area measured more than 150 mm/24h. The operational ALADIN forecast severely underestimated (the forecast maximum was below 10mm) the rainfall over the Istria peninsula during the night from September 24 to 25, 2010 (Fig 17a). The wind field is shown for 00 UTC on September 25, 2010 the time that is close to the period of the maximum rain intensity. The parallel suite rainfall structures were slightly better (Fig 17b) than the operational one, with second maximum of rainfall over the Istria peninsula reaching more than 35mm, but the predicted rainfall amount was far below the measured one. It is assumed that the observed severe precipitation was caused by convective activity supported by the synoptic conditions and/or local conditions that were not represented correctly in the initial conditions (as suggested by the result of the data assimilation forecast, Fig 17b) or the model was not able to represent its development. The last hypothesis has been tested using the non-hydrostatic set-ups of the ALADIN model and in higher resolution using different initial and boundary conditions (from 8 km runs with or without 3Dvar and using ARPEGE and IFS initial and boundary conditions) One run generated additional intensive rain band over an island 25 km east of Pula.

(a) (b)

78 Climate Change and Regional/Local Responses

(c) (d)

Palagruža (a), Dubrovnik (b), Hvar (c) and Mljet (d) for 7th-9th Nov 2010.

cases is exemplarily illustrated by two cases.

**3.5. Heavy rain**

)

**Figure 16.** Comparison of measured (purple), 2 km (full lines) and 8 km resolution wind speed forecast (dashed) for

Large amounts of rain (> 100 mm) that fall during a short time period (< 6 hours) can cause flash floods, especially in places where the terrain supports accumulation of water as on low slope coast below a mountain. The ability of the operational ALADIN model to predict such

On September 25, 2010, just after midnight an intensive rain hit Pula city on the southern part of the Istria Peninsula, Croatia. The rain was intensive for several hours and the rainfall rate measured by the ombrograph reached 43.9 mm per hour. Several rain-gauges in the area measured more than 150 mm/24h. The operational ALADIN forecast severely underestimated (the forecast maximum was below 10mm) the rainfall over the Istria peninsula during the night from September 24 to 25, 2010 (Fig 17a). The wind field is shown for 00 UTC on September 25, 2010 the time that is close to the period of the maximum rain intensity. The parallel suite rainfall structures were slightly better (Fig 17b) than the operational one, with second maximum of rainfall over the Istria peninsula reaching more than 35mm, but the predicted rainfall amount was far below the measured one. It is assumed that the observed severe precipitation was caused by convective activity supported by the synoptic conditions and/or local conditions that were not represented correctly in the initial conditions (as suggested by the result of the data assimilation forecast, Fig 17b) or the model was not able to represent its development. The last hypothesis has been tested using the non-hydrostatic set-ups of the ALADIN model and in higher resolution using different initial and boundary conditions (from 8 km runs with

**Figure 17.** Accumulated precipitation from the operational (a) and data assimilation (b) forecast for 25th Sep 2010.

The flash flood hit Dubrovnik, Croatia, on November 22, 2010. The rain gauge measurements exceeded 100 mm/24h in the area and the one in Dubrovnik measured 161.4 mm/24h (the highest on record by that time), with a peak intensity of 71.5 mm/h (Fig 18b). The flood water level reached 1.5 m caused damage and endangered lives.

The operational and parallel suite model results forecast large 24 hour accumulated rainfall amounts in the area around Dubrovnik, with the maxima over 100 mm located above the surrounding areas of Montenegro and Bosnia and Herzegovina. The results of the highresolution (2 km) non-hydrostatic ALADIN model run showed a dependency of the forecasts on the input data from initial and lateral boundary conditions. The position and time of the precipitation maxima in the high resolution ALADIN output fields (Fig 18a) was similar to the lower resolution (8 km) run used for initial and lateral boundary conditions. The peak intensity of precipitation in all model runs was late (3 hours in runs coupled to Arpege, 6 hours for IFS). The reason for this delay is yet unknown and could be attributed to the absence of measure‐ ments over the southern Adriatic.

**Figure 18.** Accumulated 24 h precipitation from the 2 km resolution run (a) and measured (b) for 23rd Nov 2010.

#### **3.6. Fog**

Although far less dramatic than the previous ones, fog is also a weather event that can have a large impact since it is a nuisance to the traffic. The results of the operational model forecasts were not satisfactory for 2m temperature and cloudiness in a case of fog and low stratus accompanied with a temperature inversion in December 2004. The problems remained through several weeks of December 2004, as fog and low stratus remained covering the inland parts of Croatia, as well as Hungary and other valleys surrounding eastern Alps.

measured (Fig 19). Actually, the measured 2m temperature changed only slightly during the day (0-2°C), as well as from one day to another. This behaviour has made the 2m temperature

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

81

**Figure 21.** As Fig 20, but obtained from a simulation using the Xu-Randall cloudiness (Xu & Randall, 1996) scheme.

These results have encouraged an extensive case study to find a model set-up that would produce an acceptable forecast of the low clouds as well as 2m temperature diurnal cycle (Tudor, 2010). The model initially recognizes the existence of the temperature inversion and the layer of air adjacent to the surface is almost saturated. But the cloud scheme diagnoses little cloudiness or fog, far less than exists in reality (Fig 20). Consequently, the radiation scheme heats the ground and breaks the inversion, making the situation even worse. In reality, a fog or low stratus layer keeps the longwave radiation from ground and reflects the shortwave radiation keeping the temperature of the surface layer almost constant. Longwave radiation flux divergence and shortwave reflection at the top of the cloud supports cooling and forma‐

The study (Tudor 2010) has shown that the combination of Xu-Randall cloud diagnosing scheme (Xu & Randall, 1996), new critical relative humidity profile, maximum overlap and operational (Ritter and Geleyn, 1982) radiation scheme produces forecast of low clouds and 2 m temperature of 1-3°C that is closest to the observed values. This result stayed nearly the same even when the more sophisticated radiation schemes were used. The introduction of prognostic microphysics improved the results too by increasing cloud coverage by 5-10% and reducing the temperature amplitude only slightly (up to 0.3°C). The introduction of more physical semi-Lagrangian horizontal diffusion (SLHD) has improved the fog and low stratus forecast close to terrain slopes by increasing cloud coverage by up to 50%. The numerical diffusion mixes the temperature and moisture along model levels that follow orography. Consequently, in the model the cloud layer from the valley is mixed with the clear air on the mountain ridge and the cloud in the valley reduces. Replacement of numerical diffusion with

SLHD reduces this mixing and supports cloud formation on slopes of the valleys.

forecast as useless as the cloudiness forecast.

tion of a temperature inversion.

**Figure 19.** Comparison of the model 2m temperature evolution with the data measured (violet) in Križevci meteoro‐ logical station. Exp1 to exp5 are obtained with the operational radiation scheme, and in exp6 to exp10 a new, en‐ hanced, radiation scheme has been used. Experiments with random maximum overlap are shown with circles and those with random overlap are shown with squares. Open marks show experiments with the operational critical rela‐ tive humidity profile and full marks show the experiments with the new critical relative humidity profile.

**Figure 20.** The operational forecast of low (shades of red), medium (shades of violet) and high (shades of blue) cloudi‐ ness for 06 UTC on December 15, 2004.

Subsequent operational model forecasts kept underestimating the cloud coverage with 30-60% on average when compared to observations by the crew from the synoptic stations. Conse‐ quently, the amplitude of the 2 m temperature diurnal cycle was much larger (5-10°C) than

**3.6. Fog**

80 Climate Change and Regional/Local Responses

Although far less dramatic than the previous ones, fog is also a weather event that can have a large impact since it is a nuisance to the traffic. The results of the operational model forecasts were not satisfactory for 2m temperature and cloudiness in a case of fog and low stratus accompanied with a temperature inversion in December 2004. The problems remained through several weeks of December 2004, as fog and low stratus remained covering the inland

**Figure 19.** Comparison of the model 2m temperature evolution with the data measured (violet) in Križevci meteoro‐ logical station. Exp1 to exp5 are obtained with the operational radiation scheme, and in exp6 to exp10 a new, en‐ hanced, radiation scheme has been used. Experiments with random maximum overlap are shown with circles and those with random overlap are shown with squares. Open marks show experiments with the operational critical rela‐

**Figure 20.** The operational forecast of low (shades of red), medium (shades of violet) and high (shades of blue) cloudi‐

Subsequent operational model forecasts kept underestimating the cloud coverage with 30-60% on average when compared to observations by the crew from the synoptic stations. Conse‐ quently, the amplitude of the 2 m temperature diurnal cycle was much larger (5-10°C) than

ness for 06 UTC on December 15, 2004.

tive humidity profile and full marks show the experiments with the new critical relative humidity profile.

parts of Croatia, as well as Hungary and other valleys surrounding eastern Alps.

**Figure 21.** As Fig 20, but obtained from a simulation using the Xu-Randall cloudiness (Xu & Randall, 1996) scheme.

measured (Fig 19). Actually, the measured 2m temperature changed only slightly during the day (0-2°C), as well as from one day to another. This behaviour has made the 2m temperature forecast as useless as the cloudiness forecast.

These results have encouraged an extensive case study to find a model set-up that would produce an acceptable forecast of the low clouds as well as 2m temperature diurnal cycle (Tudor, 2010). The model initially recognizes the existence of the temperature inversion and the layer of air adjacent to the surface is almost saturated. But the cloud scheme diagnoses little cloudiness or fog, far less than exists in reality (Fig 20). Consequently, the radiation scheme heats the ground and breaks the inversion, making the situation even worse. In reality, a fog or low stratus layer keeps the longwave radiation from ground and reflects the shortwave radiation keeping the temperature of the surface layer almost constant. Longwave radiation flux divergence and shortwave reflection at the top of the cloud supports cooling and forma‐ tion of a temperature inversion.

The study (Tudor 2010) has shown that the combination of Xu-Randall cloud diagnosing scheme (Xu & Randall, 1996), new critical relative humidity profile, maximum overlap and operational (Ritter and Geleyn, 1982) radiation scheme produces forecast of low clouds and 2 m temperature of 1-3°C that is closest to the observed values. This result stayed nearly the same even when the more sophisticated radiation schemes were used. The introduction of prognostic microphysics improved the results too by increasing cloud coverage by 5-10% and reducing the temperature amplitude only slightly (up to 0.3°C). The introduction of more physical semi-Lagrangian horizontal diffusion (SLHD) has improved the fog and low stratus forecast close to terrain slopes by increasing cloud coverage by up to 50%. The numerical diffusion mixes the temperature and moisture along model levels that follow orography. Consequently, in the model the cloud layer from the valley is mixed with the clear air on the mountain ridge and the cloud in the valley reduces. Replacement of numerical diffusion with SLHD reduces this mixing and supports cloud formation on slopes of the valleys.

## **4. Conclusion**

This chapter gives an overview on weather forecasting using the set-up of the NWP model ALADIN that is used for operational weather forecast in CMHS as an example for operation weather forecasting. ALADIN is a state-of-the-art modern NWP model. Using ALADIN we exemplarily discuss short-comings and challenges in modern operational weather forecasting. A high-resolution LAM is intended to predict the sub-synoptic weather features forced by topography or other local characteristics that can be absent in the main synoptic pattern. Successful prediction of these small-scale features enables usage of the LAM forecast in predicting the conditions important for the flight safety, vehicle road safety or navigation at sea. The operational suite has to be tuned in order to predict the high-impact weather events of local character that could be missing in the large scale forecast. The domain properties as well as the forecast model complexity are formed according to the needs of the forecast users and the computing capabilities.

useful in forecasting other phenomena could benefit from the high resolution, such as local convective storms. The windstorm that is a consequence of a local pressure disturbance requires a full forecast run in 2 km resolution using more complex NH ALADIN model setup (Tudor & Ivatek-Šahdan, 2010). Several air crash investigations have revealed that this "fullrun" high resolution ALADIN forecast enables the prediction of lee waves and zones of

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

83

These flaws of the operational suite have encouraged the introduction of the 2 km resolution 24 hour forecast with NH ALADIN set-up using the complete set of physics parameterizations to the operational suite at the beginning of summer 2011. Several case studies of high-impact weather phenomena have been used to set-up the model configuration used operationally.

Martina Tudor, Stjepan Ivatek-Šahdan, Antiono Stanešić, Kristian Horvath and Alica Bajić

[1] ALADIN International Team, (1997) The ALADIN project: Mesoscale modelling seen as a basic tool for weather forecasting and atmospheric research. *WMO Bull*., 46, 317–

[2] Bajić, A. (2003) Očekivani režim strujanja vjetra na autocesti Sv. Rok (jug) – Masleni‐

[3] Bajić, A., Ivatek-Šahdan, S. & Horvath, K. (2007) Spatial distribution of wind speed in

[4] Bajić, A., Ivatek-Šahdan, S., Žibrat, Z. (2008) ANEMO-ALARM iskustva operativne

[5] Berre L. (2000) Estimation of Synoptic and Mesoscale Forecast Error Covariances in a

[6] Bölöni, G., Horvath K. (2010) Diagnosis and tuning of the background error statistics

[7] Branković, Č., Matjačić, B., Ivatek-Šahdan, S. & Buizza, R. (2007) Dynamical down‐ scaling of ECMWF EPS forecasts applied to cases of severe weather in Croatia.

Croatia obtained using the ALADIN model. *Cro. Met. J.* 42 , 67–77.

primjene prognoze smjera i brzine vjetra. *GIU Hrvatski cestar.* 109-114

Limited-Area Model, *Mon. Wea. Rev.* 128, 644-667.

*ECMWF RD Technical Memorandum* 507, 38 pp.

in a variational data assimilation system. *Időjárás* 114, 1-19.

increased turbulence as well as icing zones.

**Author details**

**References**

324.

ca. *Građevinar*, 55, 149-158.

The results from these studies have provided encouraging results.

Croatian Meteorological and Hydrological Service, Croatia

High impact weather events can be of local origin, or determined by small scale characteristics particular to a certain region. Croatia is not spared from severe weather events, such as severe windstorms (bura andjugowind),torrentialrainandflashfloods. Suchevents areoftenforecast in terms of probability via the ensemble prediction system of a global model. The severity of the event can be controlled by mechanisms not resolved by the global model, and dynamical adaptation of the individual EPS members improves the predictability of the severe weather events (Branković et al., 2007). The NWP section of CMHS has invested most effort in the deterministichigh-resolutionforecastbasedonanalysisofnumeroushighimpactweathercases.

The operational 8 km resolution forecast is run using DFI with the initial file from the global model ARPEGE since most of the measured data is assimilated in the global model. The dataassimilation suite has been established and runs parallel to the operational suite forecast. This parallel run has revealed the benefits of doing the data-assimilation in higher resolution (8 km) than in the global model (see the Pula flash flood case).

Forecasting the formation and distribution of fog is controlled by subtle variations in humidity and temperature fields and require maintenance of a sensitive balance between the advection, radiation and turbulent fluxes. The profound influence of topography on the atmospheric flow reflects itself in a number of features that affect local weather, such as upslope, downslope and gap winds, coastal barrier jets as well as land sea breezes.

The operational high-resolution (2 km) forecast of the 10 m wind is achieved through the dynamical adaptation method adapted from Žagar & Rakovec (1999). The prediction of severe wind variability and strength is improved in weather situations with strong flow over complex topography, like bura in Croatia (Ivatek-Šahdan & Tudor, 2004). This high resolution wind forecast has provided warnings in numerous cases severe wind as well as to forecast the road conditions, assess the wind energy potential (Horvath et al., 2011) and has been used in numerous applications.

The forecast abilities of HRDA and similar packages are limited to the severe windstorms related to a synoptic forcing. Omitting moist and radiation processes prevents it from being useful in forecasting other phenomena could benefit from the high resolution, such as local convective storms. The windstorm that is a consequence of a local pressure disturbance requires a full forecast run in 2 km resolution using more complex NH ALADIN model setup (Tudor & Ivatek-Šahdan, 2010). Several air crash investigations have revealed that this "fullrun" high resolution ALADIN forecast enables the prediction of lee waves and zones of increased turbulence as well as icing zones.

These flaws of the operational suite have encouraged the introduction of the 2 km resolution 24 hour forecast with NH ALADIN set-up using the complete set of physics parameterizations to the operational suite at the beginning of summer 2011. Several case studies of high-impact weather phenomena have been used to set-up the model configuration used operationally. The results from these studies have provided encouraging results.

## **Author details**

**4. Conclusion**

82 Climate Change and Regional/Local Responses

and the computing capabilities.

numerous applications.

than in the global model (see the Pula flash flood case).

gap winds, coastal barrier jets as well as land sea breezes.

This chapter gives an overview on weather forecasting using the set-up of the NWP model ALADIN that is used for operational weather forecast in CMHS as an example for operation weather forecasting. ALADIN is a state-of-the-art modern NWP model. Using ALADIN we exemplarily discuss short-comings and challenges in modern operational weather forecasting. A high-resolution LAM is intended to predict the sub-synoptic weather features forced by topography or other local characteristics that can be absent in the main synoptic pattern. Successful prediction of these small-scale features enables usage of the LAM forecast in predicting the conditions important for the flight safety, vehicle road safety or navigation at sea. The operational suite has to be tuned in order to predict the high-impact weather events of local character that could be missing in the large scale forecast. The domain properties as well as the forecast model complexity are formed according to the needs of the forecast users

High impact weather events can be of local origin, or determined by small scale characteristics particular to a certain region. Croatia is not spared from severe weather events, such as severe windstorms (bura andjugowind),torrentialrainandflashfloods. Suchevents areoftenforecast in terms of probability via the ensemble prediction system of a global model. The severity of the event can be controlled by mechanisms not resolved by the global model, and dynamical adaptation of the individual EPS members improves the predictability of the severe weather events (Branković et al., 2007). The NWP section of CMHS has invested most effort in the deterministichigh-resolutionforecastbasedonanalysisofnumeroushighimpactweathercases. The operational 8 km resolution forecast is run using DFI with the initial file from the global model ARPEGE since most of the measured data is assimilated in the global model. The dataassimilation suite has been established and runs parallel to the operational suite forecast. This parallel run has revealed the benefits of doing the data-assimilation in higher resolution (8 km)

Forecasting the formation and distribution of fog is controlled by subtle variations in humidity and temperature fields and require maintenance of a sensitive balance between the advection, radiation and turbulent fluxes. The profound influence of topography on the atmospheric flow reflects itself in a number of features that affect local weather, such as upslope, downslope and

The operational high-resolution (2 km) forecast of the 10 m wind is achieved through the dynamical adaptation method adapted from Žagar & Rakovec (1999). The prediction of severe wind variability and strength is improved in weather situations with strong flow over complex topography, like bura in Croatia (Ivatek-Šahdan & Tudor, 2004). This high resolution wind forecast has provided warnings in numerous cases severe wind as well as to forecast the road conditions, assess the wind energy potential (Horvath et al., 2011) and has been used in

The forecast abilities of HRDA and similar packages are limited to the severe windstorms related to a synoptic forcing. Omitting moist and radiation processes prevents it from being Martina Tudor, Stjepan Ivatek-Šahdan, Antiono Stanešić, Kristian Horvath and Alica Bajić

Croatian Meteorological and Hydrological Service, Croatia

## **References**


[8] Branković, Č., Matjačić, B., Ivatek-Šahdan, S. & Buizza, R. (2008) Downscaling of ECMWF Ensemble Forecasts for Cases of Severe Weather: Ensemble Statistics and Cluster Analysis. *Mon.Wea. Rev.* 136, 3323–3342.

[22] Gerard, L. (2007) An integrated package for subgrid convection, clouds and precipi‐ tation compatible with the meso-gamma scales. *Quart. J. Roy. Meteor. Soc*., 133, 711–

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

85

[23] Gerard, L. & Geleyn, J.-F. (2005) Evolution of a subgrid deep convection parametriza‐ tion in a limited area model with increasing resolution. *Quart. J. Roy. Meteor. Soc.*,

[24] Gerard, L., Piriou, J.-M., Brožková, R., Geleyn, J.-F. & Banciu, D. (2009) Cloud and Precipitation Parameterization in a Meso-Gamma-Scale Operational Weather Predic‐

[25] Giard, D. and Bazile, E. (2000) Implementation of a new assimilation scheme for soil and surface variables in a global NWP model. *Mon. Wea. Rev.* 128, 997-1015.

[26] Gospodinov, I., Spiridonov, V., & Geleyn, J.-F., (2001) Second order accuracy of twotime-level semi-Lagrangian schemes. *Quart. J. Roy. Meteor. Soc.*, 127, 1017–1033. [27] Grisogono, B., Belušić, D. (2009) A review of recent advances in understanding the meso- and microscale properties of the severe Bora wind. *Tellus* 61A, 1–16.

[28] Grubišić, V. (2004) Bora-driven potential vorticity banners over the Adriatic*. Quart. J.*

[29] Haugen, J. E. & Machenhauer, B. (1993) A spectral limited-area model formulation with time-dependent boundary conditions applied to the shallow-water equations.

[30] Hdidou, F.Z., (2006) Test of an assimilation suite for ALBACHIR based on the varia‐

[31] Hólm, E.V., 2008: Lecture notes on assimilation algorithms. *ECMWF, European Centre*

[32] Hollingsworth, F. Rabier and M. Fisher, 1998: The ECMWF implementation of threedimensional variational assimilation (3d-Var). I: Formulation. *Quart. J. Roy. Meteor.*

[33] Horvath, K., Bajić, A., & Ivatek-Šahdan, S. (2011) Dynamical Downscaling of Wind Speed in Complex Terrain Prone To Bora-Type Flows. *J. Appl. Meteor. Climatol*., 50,

[34] Horvath, K., Ivatek-Šahdan, S., Ivančan-Picek, B. & Grubišić, V. (2009) Evolution and structure of two severe cyclonic bora events: contrast between the northern and

[35] Ivatek-Šahdan, S. & Ivančan-Picek, B. (2006) Effects of different initial and boundary conditions in ALADIN/HR simulations during MAP IOPs. Meteorol. Z. 15, 187–197.

730.

131, 2293–2312.

tion Model. *Mon. Wea. Rev.,* 137, 3960–3977.

*Roy. Meteor. Soc. 130*, 2571-2603.

*Mon. Wea. Rev.*, 121, 2618–2630.

*Soc.,* 124, 1783-1807.

1676–1691.

tional technique. available at www.wmo.int

*for Medium-Range Weather Forecasts,* Reading, England, 30 pp.

southern Adriatic. *Weather and forecasting* 24, 946-964.


[22] Gerard, L. (2007) An integrated package for subgrid convection, clouds and precipi‐ tation compatible with the meso-gamma scales. *Quart. J. Roy. Meteor. Soc*., 133, 711– 730.

[8] Branković, Č., Matjačić, B., Ivatek-Šahdan, S. & Buizza, R. (2008) Downscaling of ECMWF Ensemble Forecasts for Cases of Severe Weather: Ensemble Statistics and

[9] Brozkova R, Klarić, D., Ivatek-Šahdan, S., Geleyn, J.F., Casse, V., Siroka, M., Radnoti, G., Janousek, M., Stadlbacher, K., Seidl, H. (2001) DFI blending: an alternative tool

for preparation of the initial conditions for LAM. WGNE *Blue Book.* 31, 1.7-1.8.

*pean Centre for Medium-Range Weather Forecasts*, Reading, England, 59 pp.

[10] Catry B., Geleyn J.-F., Tudor M., Bénard P. & Trojakova A. (2007). Flux conservative thermodynamic equations in a mass-weighted framework. *Tellus* 59A, pp 71–79 [11] Courtier, 1999: Data assimilation concepts and methods, *ECMWF lecture notes, Euro‐*

[12] Davies, H. C. (1976) A lateral boundary formulation for multilevel prediction mod‐

[13] Davies, H. (1983) Limitations of some common lateral boundary schemes used in re‐

[14] Durran, D.R. (1999) Numerical methods for wave equations in geophysical fluid dy‐

[15] Geleyn J.-F. (1988). Interpolation of wind, temperature and humidity values from

[16] Geleyn, J.-F., Bazile, E., Bougeault, P., Déqué, M., Ivanovici, V., Joly, A., Labbé, L., Piedélièvre, J.-P., Piriou, J.-M., Royer, J.-F. (1994) Atmospheric parametrizations schemes in Meteo-France's ARPEGE NWP model. *ECMWF seminar proceedings on*

[17] Geleyn J.-F. & Hollingsworth A (1979) An economical analytical method for the com‐ putation of the interaction between scattering and line absorption of radiation. *Beitr*

[18] Geleyn J.-F., Benard P. & Fournier, R. (2005a) A general-purpose extension of the Malkmus band-model average equivalent width to the case of the Voigt line profile.

[19] Geleyn J.-F., Fournier R., Hello G., Pristov N. (2005b) A new 'bracketing' technique for a flexible and economical computation of thermal radiative fluxes, scattering ef‐ fects included, on the basis the Net Exchanged Rate (NER) formalism. WGNE Blue

[20] Geleyn J.-F., Vana F., Cedilnik J., Tudor M. & Catry B. (2006). An intermediate solu‐ tion between diagnostic exchange coefficients and prognostic TKE methods for verti‐

[21] Geleyn J.-F., Catry B., Bouteloup Y. & Brožkova, R. (2008). A statistical aproach for sedimentation inside a microphysical precipitation scheme. *Tellus* 60A,649–662

model levels to the height of measurement. *Tellus,* 40A, pp.347–351

*Parametrization of sub-grid scale physical processes,* pp. 385-402.

Cluster Analysis. *Mon.Wea. Rev.* 136, 3323–3342.

els. *Quart. J. Roy. Meteor. Soc.*, 102, 405–418.

namics. *Springer*, pp 465.

84 Climate Change and Regional/Local Responses

*Phys Atmos* 52:1–16.

Book

*Quart. J. Roy. Meteor. Soc.* 131:2757–2768

cal turbulent transport. *WGNE Blue Book*

gional NWP models. *Mon. Wea. Rev.*, 111, 1002–1012.


[36] Ivatek-Šahdan, S. & Tudor M. (2004) Use of high-resolution dynamical adaptation in operational suite and research impact studies. *Meteorol Z* 13(2):1–10

[51] Robert, A. (1982) A semi-Lagrangian and semi-implicit numerical integration scheme

Forecasting Weather in Croatia Using ALADIN Numerical Weather Prediction Model

http://dx.doi.org/10.5772/55698

87

[52] Simmons, A.J., Burridge, D.M. 1981: An Energy and Angular-Momentum Conserv‐ ing Vertical Finite-Difference Scheme and Hybrid Vertical Coordinates. *Mon. Wea.*

[53] Staniforth A, Côté J. (1991) Semi-Lagrangian integration schemes for atmospheric

[54] Termonia, P. (2008) Scale-selective digital filtering initialization. *Mon. Wea. Rev.*, 136,

[55] Tudor M. (2010). Impact of horizontal diffusion, radiation and cloudiness parameter‐ ization schemes on fog forecasting in valleys. *Met. Atm. Phy.* Vol.108, pp. 57-70.

[56] Tudor, M. (2011). The meteorological aspects of the DART field experiment and pre‐

[57] Tudor, M. & Ivatek-Šahdan, S. (2002) The MAP-IOP 15 case study. *Cro. Met. J*. 37,

[58] Tudor, M. & Ivatek-Šahdan, S. (2010) The case study of bura of 1st and 3rd February

[59] Tudor, M., Termonia, P., (2010) Alternative formulations for incorporating lateral boundary data into limited-area models. *Mon. Wea. Rev.* 138, pp. 2867-2882.

[60] Unden, P., Rontu, L., Jarvinen, H., Lynch, P., Calvo, J., Cats, G., Cuxart, J., Eerola, K., Fortelius, C., Garcia-Moya, J.A., Jones, C., Lenderlink, G., McDonald, A., McGrath, R., Navascues, B., Nielsen, N., Odegaard, V., Rodriguez, E., Rummukainen, M., Room, R., Sattler, K., Hansen Sass, B., Savijarvi, H., Schreur, B., Sigg, R., The, H., Tijm, S., (2002) HIRLAM-5 scientific documentation. *SMHI, S-601 76 Norrkoping, Swe‐*

[61] Váňa F., Bénard P., Geleyn J.-F., Simon A. & Seity Y. (2008). Semi-Lagrangian advec‐ tion scheme with controlled damping–an alternative way to nonlinear horizontal dif‐ fusion in a numerical weather prediction model. *Quart. J. Roy. Meteor. Soc*, Vol.134,

[62] Wang, Y., Bellus, M., Wittmann, C., Steinheimer, M., Weidle, F., Kann, A., Ivatek-Šahdan, S., Tian, W., Ma, X., Tascu, S., Bazile, E. (2011) The Central European limit‐ ed-area ensemble forecasting system: ALADIN-LAEF. *Quart. J. Roy. Meteor. Soc* 137,

[63] Warner, T., Peterson, R., Treadon, R. (1997) Atutorial on lateral boundary conditions as a basic and potentially serious limitation to regional numerical weather predic‐

for the primitive equations. *J. Meteor. Soc. Japan*, 60, 319-325.

models – A review. *Mon. Wea. Rev*. 119 2206–2223.

liminary results. *Cro. Met. J.* 44/45, 31-46.

*den*, pp 146, available at www.hirlam.org

tion. *Bull. Amer. Meteor. Soc*., 78, 2599–2617.

2007, *Meteorol. Z.,* 19, pp. 453-466.

*Rev.* 109, 758–766.

5246–5255.

1-14.

pp. 523–537.

483-502.


[51] Robert, A. (1982) A semi-Lagrangian and semi-implicit numerical integration scheme for the primitive equations. *J. Meteor. Soc. Japan*, 60, 319-325.

[36] Ivatek-Šahdan, S. & Tudor M. (2004) Use of high-resolution dynamical adaptation in

[37] Kann A., Seidl H., Wittmann C. & Haiden T. (2009) Advances in predicting continen‐

[38] Kessler E. (1969) On distribution and continuity of water substance in atmospheric

[39] Louis J.-F., Tiedke M. & Geleyn J.-F. (1982) A short history of PBL parameterization at ECMWF. In: *Proceedings from ECMWF workshop on planetary boundary layer parame‐*

[40] Lynch, P., X.-Y. Huang (1994) Diabatic Initialization using recursive filters. – *Tellus*

[41] Machenhauer, B., J.E. Haugen (1987) Test of a spectral limited area shallow water model with timedependent lateral boundary conditions and combined normal mode/ semi-lagrangian time integration schemes. *Techniques for Horizontal Discretization in*

[42] Mahfouf, J.,F., 1991: Analysis of soil moisture from near-surface variables: A feasibili‐

[43] Mlawer E.J., Taubman S.J., Brown P.D., Iacono M.J. & Clough S.A. (1997) Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for

[44] Morcrette J.-J. (1989) Description of the radiation scheme in the ECMWF Model. *Tech*

[45] Noilhan, J., Planton, S. (1989) A simple parapetrization of land surface processes for

[46] Pauluis O., Emanuel K. (2004) Numerical instability resulting from infrequent calcu‐

[48] Radnóti, G. (1995) Comments on ''A spectral limited-area formulation with time-de‐ pendent boundary conditions applied to the shallow-water equations". *Mon. Wea.*

[49] Redelsperger J.L., Mahé F., Carlotti P. (2001). A simple and general subgrid model suitable both for surface layer and free-stream turbulence. *Bound.-Layer Meteor.* 101,

[50] Ritter B. & Geleyn J.-F. (1992) A comprehensive radiation scheme for numerical weather prediction models with potential applications in climate simulations. *Mon*

[47] Pielke, R.A. (2002) Mesoscale meteorological modelling. *Academic press*, pp 676.

*Numerical Weather Prediction Models, 2–4 November 1987*, ECMWF, 361–377.

tal low stratus with a regional NWP model. *Wea. Forecasting*, 25, 290–302.

operational suite and research impact studies. *Meteorol Z* 13(2):1–10

circulations. *Meteorol Monogr Am Meteorol Soc* 10(32):84

*terization, 25–27 November 1981*, pp 59–79

ty study. *J. Appl. Meteor.*, 30, 1534-1547

*Memo 165, ECMWF*, 26 pp

*Rev*., 123, 3122– 3123.

*Wea. Rev* 120:303–325

pp.375–408.

the longwave. *J Geophys Res* 102D:16663–16682

meteorological models. *Mon. Wea. Rev.* 117, 536-549.

lation of radiative heating. Mon Weather Rev 132:673–686

46A, 583–597.

86 Climate Change and Regional/Local Responses


[64] Xu K.-M. & Randall D.A. (1996) A semi-empirical cloudiness parameterization for use in climate models. *J Atmos Sci* 53:3084–3102

**Chapter 3**

**Provisional chapter**

**Measurements and Observations of Meteorological**

**Measurements and Observations of Meteorological**

In the presence of dust, smoke, fog, haze or pollution, meteorological visibility is reduced. This reduction constitutes a common and vexing transportation problem for different public

First, low visibility is obviously a problem of traffic safety. Road crashes which occur in fog are generally more severe as the average crash [1]. According to NOAA [2], in the United States there are approximately 700 annual fog-related fatalities, defined as occurring

a smaller country, with over 100 annual fatalities attributed to low visibility, defined as

and significant problems on Northern America and French highways. The combination of fog and smoke presence on a motorway was the cause of dramatic pile-ups in France, e.g. on the A10 in 2002 near Coulombiers (58 vehicles involved, 40 injured, 8 deaths). Indeed, even if the origin of both phenomena differs, the combination of their mutual effect on the visibility is exponential, which leads to close to zero visibility areas. It should be stressed that the solution lies not necessarily in better low visibility detection but in drivers' response to fog that they encounter. Indeed, the behavior of drivers in fog is often inappropriate (e.g., reduced headways, altered reaction times) but to understand the origins of these dangerous behaviors is difficult [3]. Different countermeasures have been tested to mitigate the impact of critically reduced visibility [4]. The California San Joaquin and Sacramento Valley regions are particularly adequate test-beds for such measures, because of the well-known Tule fog phenomenon. In the Stockton area of Caltrans District 10, the Caltrans Automated Warning System (CAWS) employs roadside weather stations and visibility meters to provide automated detection [5]. In District 6, Caltrans has installed the "Fog Pilot" system, which

> ©2012 Hautière et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Hautière et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2013 Hautière et al. licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

<sup>4</sup> mile. Fog constitutes an equally important issue in France,

<sup>4</sup> mile). Indeed, fog causes similar

**Visibility at ITS Stations**

**Visibility at ITS Stations**

http://dx.doi.org/10.5772/55697

when visibility is less than a <sup>1</sup>

10.5772/55697

**1. Introduction**

Nicolas Hautière, Raouf Babari, Eric Dumont,

Nicolas Hautière, Raouf Babari, Eric Dumont, Jacques Parent Du Chatelet and Nicolas Paparoditis

Additional information is available at the end of the chapter

authorities in multiple countries throughout the world.

occurring when visibility is less than a 400 meters (≈ <sup>1</sup>

Additional information is available at the end of the chapter

Jacques Parent Du Chatelet and Nicolas Paparoditis


**Provisional chapter**

## **Measurements and Observations of Meteorological Visibility at ITS Stations Measurements and Observations of Meteorological Visibility at ITS Stations**

Nicolas Hautière, Raouf Babari, Eric Dumont, Jacques Parent Du Chatelet and Nicolas Paparoditis Nicolas Hautière, Raouf Babari, Eric Dumont, Jacques Parent Du Chatelet and Nicolas Paparoditis

Additional information is available at the end of the chapter

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/55697

## **1. Introduction**

10.5772/55697

[64] Xu K.-M. & Randall D.A. (1996) A semi-empirical cloudiness parameterization for

[65] Yessad, K (2011) Library architecture and history of the technical aspects in AR‐ PEGE/IFS, ALADIN and AROME in the cycle 37 of ARPEGE/IFS. *Meteo-France,* pp.

22, available at http://www.cnrm.meteo.fr/gmapdoc/IMG/pdf/ykarchi37t1.pdf

atia: 1961. - 1990. : 1971. - 2000. . Zagreb : CMHS *Monograph*, pp. 200.

[66] Zaninović, K., Gajić-Čapka, M., Perčec Tadić, M., Vučetić, M., Milković, J., Bajić, A., Cindrić, K., Cvitan, L., Katušin, Z., Kaučić, D., Likso, T., Lončar, E., Lončar, Ž., Mihaj‐ lović, D., Pandžić, K., Patarčić, M., Srnec, L., Vučetić, V., (2008) Climate atlas of Cro‐

[67] Žagar, M., Rakovec, J. (1999): Small-scale surface wind prediciton using dynamical

use in climate models. *J Atmos Sci* 53:3084–3102

88 Climate Change and Regional/Local Responses

adaptation. *Tellus,* Vol.51A, pp. 489–504.

In the presence of dust, smoke, fog, haze or pollution, meteorological visibility is reduced. This reduction constitutes a common and vexing transportation problem for different public authorities in multiple countries throughout the world.

First, low visibility is obviously a problem of traffic safety. Road crashes which occur in fog are generally more severe as the average crash [1]. According to NOAA [2], in the United States there are approximately 700 annual fog-related fatalities, defined as occurring when visibility is less than a <sup>1</sup> <sup>4</sup> mile. Fog constitutes an equally important issue in France, a smaller country, with over 100 annual fatalities attributed to low visibility, defined as occurring when visibility is less than a 400 meters (≈ <sup>1</sup> <sup>4</sup> mile). Indeed, fog causes similar and significant problems on Northern America and French highways. The combination of fog and smoke presence on a motorway was the cause of dramatic pile-ups in France, e.g. on the A10 in 2002 near Coulombiers (58 vehicles involved, 40 injured, 8 deaths). Indeed, even if the origin of both phenomena differs, the combination of their mutual effect on the visibility is exponential, which leads to close to zero visibility areas. It should be stressed that the solution lies not necessarily in better low visibility detection but in drivers' response to fog that they encounter. Indeed, the behavior of drivers in fog is often inappropriate (e.g., reduced headways, altered reaction times) but to understand the origins of these dangerous behaviors is difficult [3]. Different countermeasures have been tested to mitigate the impact of critically reduced visibility [4]. The California San Joaquin and Sacramento Valley regions are particularly adequate test-beds for such measures, because of the well-known Tule fog phenomenon. In the Stockton area of Caltrans District 10, the Caltrans Automated Warning System (CAWS) employs roadside weather stations and visibility meters to provide automated detection [5]. In District 6, Caltrans has installed the "Fog Pilot" system, which

©2012 Hautière et al., licensee InTech. This is an open access chapter distributed under the terms of the

provides a high-technology solution every <sup>1</sup> <sup>4</sup> mile along a 12-mile (20-km) portion of State Route 99.

10.5772/55697

91

http://dx.doi.org/10.5772/55697

meteorological visibility distance is lower than 1000 meters are obtained from 60 weather stations distributed over the entire territory. Today, about 160 meteorological stations with visibility measurements are available in real time. But as such, this data cannot be used for predicting fog events and warning road authorities and hence drivers. Indeed, the local nature of this phenomenon is not compatible with the current capacity of meteorological agencies to monitor it accurately. Some studies seek to exploit the photosensitive cells of

Measurements and Observations of Meteorological Visibility at ITS Stations

A survey has been conducted by Météo France on the French motorway networks to estimate the potential of existing CCTV networks to observe the visibility: In 2009, the French motorway network was 8,372 km long and was equipped with approximately 2,000 cameras. Accounting the fact that some are grouped together and some are dedicated to tunnel safety, a potential of 1,000 cameras available to monitor the weather was estimated. The French highway network is also equipped with cameras but they are less numerous. This whole network covers the territory quite uniformly. Consequently, a roadside sensor network constitutes a relevant mesh able to feed meteorological centers with geolocalized data.

The term Intelligent Transportation Systems (ITS) refers to information and communication technology (applied to transportation infrastructure and vehicles) that improves transportation outcomes such as transportation safety, transportation productivity, travel reliability, informed travel choices, social equity, environmental performance and network operation resilience. The recent development of real-time data exchange systems between vehicles and infrastructure allows linking operation centers, roadside sensors and vehicles by means of so-called "ITS Stations". Such technology fosters a new generation of Intelligent Transportation Systems. The objectives of this work is thus to design computer vision methods, which can be implemented into camera-based surveillance systems connected to ITS Stations, in order to detect and characterize reduced visibility conditions, so as to mitigate the risk of accidents by alerting the drivers or by computing adaptive speeds related to the

According to [11, 12], the "meteorological visibility distance" denotes the greatest distance at which a black object of a suitable dimension can be seen in the sky on the horizon, with a threshold contrast set at 5%. The meteorological visibility distance is thus a standard dimension which characterizes the opacity of the atmosphere. According to [13], the road visibility is defined as the horizontal visibility determined 1.2 m above the roadway. It may be reduced to less than 400 m by fog, precipitations or projections. Four visibility ranges are defined and are listed in Tab. 1. Based on these definitions, a visibility sensor should assign the visibility range to one of the four categories and detect the origin of the visibility reduction, i.e. it should detect fog, rain and projections. In this section, the focus is on daytime fog detection and visibility range estimation using two complementary systems.

fixed cameras to measure the visibility.

**2.2. The potential of CCTV networks**

**2.3. Intelligent Transportation System**

**3. Background on meteorological visibility**

**3.1. Visibility sensor requirements**

offered visibility distance.

In addition to the safety problem, reduced visibility is the cause of delays and disruption in air, sea and ground transportation for passengers and freight. On freeways, massive pile-ups create exceptional traffic congestions which sometimes force the operator to momentarily close the road. Fog-related road closures are not an uncommon subject for news headlines. For example, the Heathrow airport was blocked for three days during the 2006 Christmas time. Such events have of course important economic impacts [6]. According to [7], in 1974 fog was estimated to have cost over roughly £120 millions at 2010 prices on the roads of Great Britain. This number includes the cost of medical treatment, damage to vehicles and property, as well as the administrative costs of police, services and insurance, but they do not include the cost of delays to vehicle passengers not directly involved in the accident.

An ability to accurately monitor visibility helps resolve these problems. Important transportation facilities where safety is critical, such as airports, are generally instrumented for monitoring visibility with devices that are expensive and hence, scarce. Cost is precisely the reason why highway meteorological stations are seldom equipped with visibility metering devices. In this context, using already existing and ubiquitous highway cameras is of great interest, as these are low cost sensors already deployed for other purposes such as traffic monitoring [8]. Furthermore, introducing new functionalities into roadside cameras will make them multipurpose and thus more cost-effective, easing their deployment along the roads. In the United States, this potential has been identified by US DOT and was evaluated in the CLARUS Initiative [9], and these efforts may continue with the US DOT IntelliDrive program. In France, a similar initiative has been launched between Ifsttar (French institute of science and technology for transport, development and networks), Météo France (French National Weather Service) and IGN (French National Geographical Institute), three public research institutes dealing with road operation, weather monitoring and forecasting, and geography and cartography, respectively. The French initiative aims at assessing the potential of highway cameras to monitor visibility for different applications ranging from safety hazard detection to air quality monitoring. In the future, such initiatives might make it possible to monitor visibility reduction at the scale of a road itinerary. Prediction, which will soon be possible for airports [10], might even be envisioned.

## **2. Objectives**

## **2.1. Problematic**

Reduced visibility in the atmosphere is directly related to light scattering by air molecules and airborne particles. This tenet of physics is the basis of the operating principle of visibility meters. There are two types of instruments for measuring atmospheric visibility: transmissometers and scatter meters. The transmissometer extrapolates the attenuation of a light beam emitted from a source to a receiver at a known path length in order to estimate the distance for which the emitted light is attenuated by 95%. The transmissometer is also used to calibrate the scatter meter. A scatter meter assesses the dispersion of a light beam at a particular scattering angle, more often close to 40˚(forward-scatter meters). Visibility meters can measure the meteorological visibility distance up to a few tens of kilometers with an error of 20%. The annual statistics on fog occurrence in France, i.e. episodes where the meteorological visibility distance is lower than 1000 meters are obtained from 60 weather stations distributed over the entire territory. Today, about 160 meteorological stations with visibility measurements are available in real time. But as such, this data cannot be used for predicting fog events and warning road authorities and hence drivers. Indeed, the local nature of this phenomenon is not compatible with the current capacity of meteorological agencies to monitor it accurately. Some studies seek to exploit the photosensitive cells of fixed cameras to measure the visibility.

## **2.2. The potential of CCTV networks**

2 Climate Change and Regional/Local Responses

Route 99.

**2. Objectives**

**2.1. Problematic**

provides a high-technology solution every <sup>1</sup>

<sup>4</sup> mile along a 12-mile (20-km) portion of State

In addition to the safety problem, reduced visibility is the cause of delays and disruption in air, sea and ground transportation for passengers and freight. On freeways, massive pile-ups create exceptional traffic congestions which sometimes force the operator to momentarily close the road. Fog-related road closures are not an uncommon subject for news headlines. For example, the Heathrow airport was blocked for three days during the 2006 Christmas time. Such events have of course important economic impacts [6]. According to [7], in 1974 fog was estimated to have cost over roughly £120 millions at 2010 prices on the roads of Great Britain. This number includes the cost of medical treatment, damage to vehicles and property, as well as the administrative costs of police, services and insurance, but they do not include the cost of delays to vehicle passengers not directly involved in the accident.

An ability to accurately monitor visibility helps resolve these problems. Important transportation facilities where safety is critical, such as airports, are generally instrumented for monitoring visibility with devices that are expensive and hence, scarce. Cost is precisely the reason why highway meteorological stations are seldom equipped with visibility metering devices. In this context, using already existing and ubiquitous highway cameras is of great interest, as these are low cost sensors already deployed for other purposes such as traffic monitoring [8]. Furthermore, introducing new functionalities into roadside cameras will make them multipurpose and thus more cost-effective, easing their deployment along the roads. In the United States, this potential has been identified by US DOT and was evaluated in the CLARUS Initiative [9], and these efforts may continue with the US DOT IntelliDrive program. In France, a similar initiative has been launched between Ifsttar (French institute of science and technology for transport, development and networks), Météo France (French National Weather Service) and IGN (French National Geographical Institute), three public research institutes dealing with road operation, weather monitoring and forecasting, and geography and cartography, respectively. The French initiative aims at assessing the potential of highway cameras to monitor visibility for different applications ranging from safety hazard detection to air quality monitoring. In the future, such initiatives might make it possible to monitor visibility reduction at the scale of a road itinerary. Prediction, which

Reduced visibility in the atmosphere is directly related to light scattering by air molecules and airborne particles. This tenet of physics is the basis of the operating principle of visibility meters. There are two types of instruments for measuring atmospheric visibility: transmissometers and scatter meters. The transmissometer extrapolates the attenuation of a light beam emitted from a source to a receiver at a known path length in order to estimate the distance for which the emitted light is attenuated by 95%. The transmissometer is also used to calibrate the scatter meter. A scatter meter assesses the dispersion of a light beam at a particular scattering angle, more often close to 40˚(forward-scatter meters). Visibility meters can measure the meteorological visibility distance up to a few tens of kilometers with an error of 20%. The annual statistics on fog occurrence in France, i.e. episodes where the

will soon be possible for airports [10], might even be envisioned.

A survey has been conducted by Météo France on the French motorway networks to estimate the potential of existing CCTV networks to observe the visibility: In 2009, the French motorway network was 8,372 km long and was equipped with approximately 2,000 cameras. Accounting the fact that some are grouped together and some are dedicated to tunnel safety, a potential of 1,000 cameras available to monitor the weather was estimated. The French highway network is also equipped with cameras but they are less numerous. This whole network covers the territory quite uniformly. Consequently, a roadside sensor network constitutes a relevant mesh able to feed meteorological centers with geolocalized data.

## **2.3. Intelligent Transportation System**

The term Intelligent Transportation Systems (ITS) refers to information and communication technology (applied to transportation infrastructure and vehicles) that improves transportation outcomes such as transportation safety, transportation productivity, travel reliability, informed travel choices, social equity, environmental performance and network operation resilience. The recent development of real-time data exchange systems between vehicles and infrastructure allows linking operation centers, roadside sensors and vehicles by means of so-called "ITS Stations". Such technology fosters a new generation of Intelligent Transportation Systems. The objectives of this work is thus to design computer vision methods, which can be implemented into camera-based surveillance systems connected to ITS Stations, in order to detect and characterize reduced visibility conditions, so as to mitigate the risk of accidents by alerting the drivers or by computing adaptive speeds related to the offered visibility distance.

## **3. Background on meteorological visibility**

## **3.1. Visibility sensor requirements**

According to [11, 12], the "meteorological visibility distance" denotes the greatest distance at which a black object of a suitable dimension can be seen in the sky on the horizon, with a threshold contrast set at 5%. The meteorological visibility distance is thus a standard dimension which characterizes the opacity of the atmosphere. According to [13], the road visibility is defined as the horizontal visibility determined 1.2 m above the roadway. It may be reduced to less than 400 m by fog, precipitations or projections. Four visibility ranges are defined and are listed in Tab. 1. Based on these definitions, a visibility sensor should assign the visibility range to one of the four categories and detect the origin of the visibility reduction, i.e. it should detect fog, rain and projections. In this section, the focus is on daytime fog detection and visibility range estimation using two complementary systems.


**Table 1.** Ranges issued from the French standard NF P 99-320 on highway meteorology, in agreement with the international practice.

#### **3.2. Vision through the atmosphere**

The apparent luminance of the road pavement *L* is given by Koschmieder's law [14] which adds to Beer-Lambert's law a second term corresponding to the atmospheric veil:

$$L = L\_0 e^{-kd} + L\_f(1 - e^{-kd}) \tag{1}$$

10.5772/55697

93

http://dx.doi.org/10.5772/55697

f

z

θ

H

Z

*D* = 

*∂*2 *I*

(*v*−*vh* )<sup>3</sup> . The equation *<sup>∂</sup>*<sup>2</sup> *<sup>I</sup>*

physically plausible. The only useful solution is (7):

Following a change of *d* according to *v* based on (3), (2) then becomes:

*I*(*v*) = *R* − (*R* − *A*∞)(1 − *e*

− *<sup>k</sup><sup>λ</sup> <sup>v</sup>*−*vh*

*<sup>k</sup>* <sup>=</sup> <sup>2</sup>(*vi* <sup>−</sup> *vh*)

where *vi* denotes the position of the inflection point of *I*(*v*). In this manner, if *vi* > *vh*,

*<sup>V</sup>*met <sup>=</sup> <sup>3</sup>*<sup>λ</sup>*

daytime fog is detected and the parameter *k* is obtained. We deduce *V*met = <sup>3</sup>

 *kλ v* − *vh*

By taking the second derivative of *I* with respect to *v*, one obtains the following:

*<sup>∂</sup>v*<sup>2</sup> (*v*) = *<sup>k</sup>ϕ*(*v*)*<sup>e</sup>*

**C** x

y

**S** X Y

point *M* to the camera is as following:

**4.2. Daytime fog detection**

where *<sup>ϕ</sup>*(*v*) = *<sup>λ</sup>*(*R*−*A*∞)

d

**Figure 1.** Modelling of the camera within the road environment. *vh*: image line corresponding to the horizon line in the image. The relationship between the distance *d* on the ground and the distance *D* from the same

road plane

θ

D

v

image plane

M

*H*<sup>2</sup> cos2 *θ* + (*d* − *H* sin *θ*)<sup>2</sup> (4)

*<sup>v</sup>*−*vh* ) (5)

(6)

− *<sup>k</sup><sup>λ</sup>*

− 2 

*<sup>∂</sup>v*<sup>2</sup> = 0 has two solutions. The solution *k* = 0 is not

*<sup>λ</sup>* (7)

<sup>2</sup>(*vi* <sup>−</sup> *vh*) (8)

*<sup>k</sup>* [11]:

<sup>u</sup> vh

Measurements and Observations of Meteorological Visibility at ITS Stations

where *L*<sup>0</sup> denotes the intrinsic luminance of the pavement and *Lf* the atmospheric luminance. In a foggy image, the intensity *I* of a pixel is the result of the camera response function crf applied to (1). Assuming that crf is linear, (1) becomes:

$$I = \text{crf}(L) = \text{Re}^{-kd} + A\_{\infty}(1 - e^{-kd}) \tag{2}$$

where *R* is the intrinsic intensity of the pixel, i.e. the intensity corresponding to the intrinsic luminance value of the corresponding scene point and *A*∞ is the background sky intensity.

## **4. Detection and characterization of safety-critical visibility ranges**

The first system is able to detect daytime fog and estimate safety-critical visibility ranges. To be run, this system only needs an accurate geometrical calibration of the camera with respect to the road plane.

#### **4.1. Camera modelling**

Assuming that the road is locally planar, the distance of a point located at the range *d* on the roadway can be expressed in the image plane, assuming a pinhole camera model [15]:

$$d = \frac{\lambda}{(v - v\_h)}\tag{3}$$

where *<sup>λ</sup>* = *<sup>H</sup><sup>α</sup>* cos(*θ*) and *vh* <sup>=</sup> *<sup>v</sup>*<sup>0</sup> <sup>−</sup> *<sup>α</sup>* tan(*θ*). *<sup>H</sup>* and *<sup>θ</sup>* respectively denote the mounting height and the pitch angle of the camera. *α* = *<sup>f</sup> tp* is an intrinsic parameter of the camera based its focal length *f* and the size *tp* of a pixel. *v*<sup>0</sup> and *vh* respectively denote the vertical position in the image of the projection center and of the horizon line (see Fig. 1).

**Figure 1.** Modelling of the camera within the road environment. *vh*: image line corresponding to the horizon line in the image.

The relationship between the distance *d* on the ground and the distance *D* from the same point *M* to the camera is as following:

$$D = \sqrt{H^2 \cos^2 \theta + (d - H \sin \theta)^2} \tag{4}$$

#### **4.2. Daytime fog detection**

4 Climate Change and Regional/Local Responses

**3.2. Vision through the atmosphere**

practice.

to the road plane.

where *<sup>λ</sup>* = *<sup>H</sup><sup>α</sup>*

**4.1. Camera modelling**

and the pitch angle of the camera. *α* = *<sup>f</sup>*

**Visibility range index Horizontal visibility distance (m)** 1 200 to 400 2 100 to 200 3 50 to 100 4 < 50

**Table 1.** Ranges issued from the French standard NF P 99-320 on highway meteorology, in agreement with the international

The apparent luminance of the road pavement *L* is given by Koschmieder's law [14] which

<sup>−</sup>*kd* <sup>+</sup> *Lf*(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*

where *L*<sup>0</sup> denotes the intrinsic luminance of the pavement and *Lf* the atmospheric luminance. In a foggy image, the intensity *I* of a pixel is the result of the camera response function crf

*<sup>I</sup>* <sup>=</sup> crf(*L*) = *Re*−*kd* <sup>+</sup> *<sup>A</sup>*∞(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*

**4. Detection and characterization of safety-critical visibility ranges**

where *R* is the intrinsic intensity of the pixel, i.e. the intensity corresponding to the intrinsic luminance value of the corresponding scene point and *A*∞ is the background sky intensity.

The first system is able to detect daytime fog and estimate safety-critical visibility ranges. To be run, this system only needs an accurate geometrical calibration of the camera with respect

Assuming that the road is locally planar, the distance of a point located at the range *d* on the roadway can be expressed in the image plane, assuming a pinhole camera model [15]:

*<sup>d</sup>* <sup>=</sup> *<sup>λ</sup>*

focal length *f* and the size *tp* of a pixel. *v*<sup>0</sup> and *vh* respectively denote the vertical position in

the image of the projection center and of the horizon line (see Fig. 1).

cos(*θ*) and *vh* <sup>=</sup> *<sup>v</sup>*<sup>0</sup> <sup>−</sup> *<sup>α</sup>* tan(*θ*). *<sup>H</sup>* and *<sup>θ</sup>* respectively denote the mounting height

<sup>−</sup>*kd*) (1)

<sup>−</sup>*kd*) (2)

(*<sup>v</sup>* <sup>−</sup> *vh*) (3)

*tp* is an intrinsic parameter of the camera based its

adds to Beer-Lambert's law a second term corresponding to the atmospheric veil:

*L* = *L*0*e*

applied to (1). Assuming that crf is linear, (1) becomes:

Following a change of *d* according to *v* based on (3), (2) then becomes:

$$I(v) = R - (R - A\_{\infty})(1 - e^{-\frac{\hbar\lambda}{v - v\_{\hbar}}}) \tag{5}$$

By taking the second derivative of *I* with respect to *v*, one obtains the following:

$$\frac{\partial^2 I}{\partial v^2}(v) = k\varphi(v)e^{-\frac{k\lambda}{v-v\_h}} \left(\frac{k\lambda}{v-v\_h} - 2\right) \tag{6}$$

where *<sup>ϕ</sup>*(*v*) = *<sup>λ</sup>*(*R*−*A*∞) (*v*−*vh* )<sup>3</sup> . The equation *<sup>∂</sup>*<sup>2</sup> *<sup>I</sup> <sup>∂</sup>v*<sup>2</sup> = 0 has two solutions. The solution *k* = 0 is not physically plausible. The only useful solution is (7):

$$k = \frac{2(v\_i - v\_h)}{\lambda} \tag{7}$$

where *vi* denotes the position of the inflection point of *I*(*v*). In this manner, if *vi* > *vh*, daytime fog is detected and the parameter *k* is obtained. We deduce *V*met = <sup>3</sup> *<sup>k</sup>* [11]:

$$V\_{\text{met}} = \frac{3\lambda}{2(v\_i - v\_h)}\tag{8}$$

To implement this principle, a region within the image that displays minimal line-to-line gradient variation when browsed from bottom to top is identified by a segmentation algorithm. A vertical band is then selected in the detected area. Finally, taking the median intensity of each segment yields the vertical variation of the intensity of the image and the position of the inflection point. Details of the method are given in [15]. It has been applied to a sample image in Fig. 2(a). Even if there are many vehicles in the original image, the method is able to ignore them and to detect fog presence, as well as to estimate the meteorological visibility.

#### **4.3. Estimation of the visibility distance**

The previous method detects that the visibility is reduced by daytime fog and estimates its density. In the same way, methods dedicated to other meteorological phenomena quantification could be added. Nevertheless, to supervise these different methods, a generic method is needed to estimate the visibility. To achieve this aim, we proposed to compute the distance to the furthest visible object on the road surface. This distance is called the mobilized visibility distance *Vmob*, which is close to the definition of *Vmet* if a 5% contrast threshold is chosen [16]. Thus, a local contrast computation algorithm, based on Köhler's binarization technique and described in detail in [16], is applied to the image to compute local contrasts above or equal to 5%. The obtained contrast map contains objects of the road scene. A flat road may be assumed. As a matter of fact, along a top-bottom scanning line of the local contrast map starting from the horizon line, objects encountered get closer to the camera. Consequently, the algorithm consists of finding the highest point in the contrast map having a local contrast above 5%. *vc* denotes the corresponding image-line. The distance to this point can then be recovered using Eq. (3), which allows estimating *V*mob [17]:

$$V\_{\rm mob} = \frac{\lambda}{v\_{\rm c} - v\_{\rm h}} \tag{9}$$

10.5772/55697

95

http://dx.doi.org/10.5772/55697

(a) (b)

Measurements and Observations of Meteorological Visibility at ITS Stations

∆(*d*max) < 0.1*d*max (11)

*vh* > 0 (12)

< *v*<sup>0</sup> (13)

(14)

**Figure 2.** Detection algorithms applied to a fog highway image: (a) The vertical yellow curve represents the instantiation of (2); the horizontal red line represents the estimation of the visibility distance. The blue vertical segments represent the limits of

Second, the system must detect fog. Based on section 4.2, the horizon line must lie in the image. Third, the visibility system shall detect visibilities lower than *d*min (50 m in our case). To run correctly, the corresponding location of the inflection point must lie in the upper part of the image, i.e. *vi* must be lower than *v*0. Consequently, additional constraints on the

> 3*λ d*min

The admissible solutions of Eq. (14) can then be used to solve Eq. (11). To fulfill the requirements expressed in [13], we are thus able to specify relevant camera characteristics, which are partly detailed in Tab. 2. In this table, *b* denotes the diameter of the camera matrix. *H* denotes the sensor mounting height. *f* denotes the focal length of the camera optics (see Fig. 1). *tp* denotes the pixel size of the camera matrix and *dimy* the width of the matrix. *θ*

> *b* [inch] *H* [m] *f* [mm] *tp* [*µ*m] *dimy* [pix] *θ* [degree] 1/2 6 4.5 4.65 1360 28-29

<sup>&</sup>lt; *<sup>θ</sup>* <sup>&</sup>lt; tan−<sup>1</sup>

 *v*<sup>0</sup> *α* 

*vh* +

 *H* <sup>3</sup>*dmin*

From (12) and (13), the following inequation is obtained:

denotes the pitch angle of the camera.

**Table 2.** Technical solution for camera specification.

sin−<sup>1</sup>

the vertical band analyzed. (b) Map of local contrasts above 5%.

compromise between accuracy and cost [18]:

sensor are as following:

However, the image may also contain vertical objects, which do not respect the flat world assumption and alter the method. This scenario is the case in Fig. 2(b), where the vehicle lights are detected higher in the image than the road surface elements. Another step is thus needed to filter the vertical objects and correctly estimate the visibility distance. This task is achieved using a background modelling method [17].

#### **4.4. Camera specification**

First, according to the sensor requirements given in section 3.1, the visibility system shall detect visibility up to *d*max (400 m in our case). By using Eq. (3), the surface covered by a pixel at the distance *d* can be computed [15]:

$$\Delta(d) = \frac{\lambda}{\lfloor v\_h + \frac{\lambda}{d} \rfloor - v\_h} - \frac{\lambda}{\lceil v\_h + \frac{\lambda}{d} \rceil - v\_h} \tag{10}$$

where ⌊*x*⌋ designates the whole part of *x* and ⌈*x*⌉ the integer greater than or equal to *x*. We proposed this surface to be lower than 10% of *d*max (40 m in our case), which is a good

**Figure 2.** Detection algorithms applied to a fog highway image: (a) The vertical yellow curve represents the instantiation of (2); the horizontal red line represents the estimation of the visibility distance. The blue vertical segments represent the limits of the vertical band analyzed. (b) Map of local contrasts above 5%.

compromise between accuracy and cost [18]:

6 Climate Change and Regional/Local Responses

**4.3. Estimation of the visibility distance**

achieved using a background modelling method [17].

<sup>∆</sup>(*d*) = *<sup>λ</sup>* ⌊*vh* + *<sup>λ</sup>*

pixel at the distance *d* can be computed [15]:

**4.4. Camera specification**

visibility.

To implement this principle, a region within the image that displays minimal line-to-line gradient variation when browsed from bottom to top is identified by a segmentation algorithm. A vertical band is then selected in the detected area. Finally, taking the median intensity of each segment yields the vertical variation of the intensity of the image and the position of the inflection point. Details of the method are given in [15]. It has been applied to a sample image in Fig. 2(a). Even if there are many vehicles in the original image, the method is able to ignore them and to detect fog presence, as well as to estimate the meteorological

The previous method detects that the visibility is reduced by daytime fog and estimates its density. In the same way, methods dedicated to other meteorological phenomena quantification could be added. Nevertheless, to supervise these different methods, a generic method is needed to estimate the visibility. To achieve this aim, we proposed to compute the distance to the furthest visible object on the road surface. This distance is called the mobilized visibility distance *Vmob*, which is close to the definition of *Vmet* if a 5% contrast threshold is chosen [16]. Thus, a local contrast computation algorithm, based on Köhler's binarization technique and described in detail in [16], is applied to the image to compute local contrasts above or equal to 5%. The obtained contrast map contains objects of the road scene. A flat road may be assumed. As a matter of fact, along a top-bottom scanning line of the local contrast map starting from the horizon line, objects encountered get closer to the camera. Consequently, the algorithm consists of finding the highest point in the contrast map having a local contrast above 5%. *vc* denotes the corresponding image-line. The distance to

this point can then be recovered using Eq. (3), which allows estimating *V*mob [17]:

*<sup>V</sup>*mob <sup>=</sup> *<sup>λ</sup>*

However, the image may also contain vertical objects, which do not respect the flat world assumption and alter the method. This scenario is the case in Fig. 2(b), where the vehicle lights are detected higher in the image than the road surface elements. Another step is thus needed to filter the vertical objects and correctly estimate the visibility distance. This task is

First, according to the sensor requirements given in section 3.1, the visibility system shall detect visibility up to *d*max (400 m in our case). By using Eq. (3), the surface covered by a

*<sup>d</sup>* ⌋ − *vh*

where ⌊*x*⌋ designates the whole part of *x* and ⌈*x*⌉ the integer greater than or equal to *x*. We proposed this surface to be lower than 10% of *d*max (40 m in our case), which is a good

<sup>−</sup> *<sup>λ</sup>* ⌈*vh* + *<sup>λ</sup>*

*<sup>d</sup>* ⌉ − *vh*

*vc* − *vh*

(9)

(10)

$$
\Delta(d\_{\text{max}}) < 0.1 d\_{\text{max}} \tag{11}
$$

Second, the system must detect fog. Based on section 4.2, the horizon line must lie in the image. Third, the visibility system shall detect visibilities lower than *d*min (50 m in our case). To run correctly, the corresponding location of the inflection point must lie in the upper part of the image, i.e. *vi* must be lower than *v*0. Consequently, additional constraints on the sensor are as following:

$$v\_{\hbar} > 0\tag{12}$$

$$
\upsilon\_h + \frac{3\lambda}{d\_{\rm min}} < \upsilon\_0 \tag{13}
$$

From (12) and (13), the following inequation is obtained:

$$\sin^{-1}\left(\frac{H}{3d\_{\min}}\right) < \theta < \tan^{-1}\left(\frac{v\_0}{\alpha}\right) \tag{14}$$

The admissible solutions of Eq. (14) can then be used to solve Eq. (11). To fulfill the requirements expressed in [13], we are thus able to specify relevant camera characteristics, which are partly detailed in Tab. 2. In this table, *b* denotes the diameter of the camera matrix. *H* denotes the sensor mounting height. *f* denotes the focal length of the camera optics (see Fig. 1). *tp* denotes the pixel size of the camera matrix and *dimy* the width of the matrix. *θ* denotes the pitch angle of the camera.


**Table 2.** Technical solution for camera specification.

**Figure 3.** Experimental verification of camera specifications: (a) experimental setup; (b) quadratic error of calibration with respect to the number of considered pairs of points.

## **4.5. Experimental validation**

#### *4.5.1. Verification of Camera Specifications*

First, an experimental verification of the camera specifications has been carried out to check if we are able to reach the specifications. To achieve this aim, seven cones were set on a flat road section following the experimental setup of Fig. 3(a). Using the perspective projection model and the positions of the different cones, we are able to calibrate the camera:

$$\left(\lambda, v\_h\right) = \underset{n=1\dots\mathcal{T}}{\operatorname{argmin}} \sqrt{\sum\_{i=1}^{n} \left(d\_i - \frac{\lambda}{v\_i - v\_h}\right)^2} \tag{15}$$

10.5772/55697

97

http://dx.doi.org/10.5772/55697

(a)

0

00:13:00

00:16:00

00:19:00

00:22:00

00:25:00

00:28:00

00:31:00

00:34:00

00:37:00

00:40:00

00:43:00

00:46:00

00:49:00

(b)

**Figure 4.** Implementation of the designed camera-based visibility metering system: (a) sample images acquired during a fog

With the second system, we are able to monitor a whole spectrum of visibility ranges (from 0 to 10,000 m). In this system, we calibrated a response function of the contrast within the image with respect to reference visibility measurements obtained by external sensors. The camera needs not be a high resolution one. However, the calibration is more complex and

Let us consider an outdoor scene where targets are distributed continuously at increasing distances from the camera. When we assume that the surface of the targets is Lambertian,

> *L* = *ρ<sup>i</sup> E*

where *E* denotes the global illumination and *ρ<sup>i</sup>* denotes the albedo at *i*. Moreover, it is a

<sup>−</sup>*β<sup>d</sup>* <sup>≈</sup> (*ρ*<sup>2</sup> <sup>−</sup> *<sup>ρ</sup>*1)*<sup>e</sup>*

00:52:00

**Time** 

00:55:00

00:58:00

01:01:00

01:04:00

01:07:00

01:10:00

01:13:00

01:16:00

01:19:00

Measurements and Observations of Meteorological Visibility at ITS Stations

01:22:00

01:25:00

*<sup>π</sup>* (16)

*<sup>V</sup>* (17)

*π* so that the contrast of two Lambertian targets at distance

− <sup>3</sup>*<sup>d</sup>*

− <sup>3</sup>*<sup>d</sup> <sup>V</sup>* = ∆*ρe*

00:10:00

episode; (b) visibility distance estimated during this fog episode.

needs at least one fog episode.

**5.1. Contrast of a distant target**

classical assumption to set *<sup>L</sup>*<sup>∞</sup> = *<sup>E</sup>*

*d* becomes [19]:

**5. Monitoring of meteorological visibility**

the luminance *L* at each point *i* of the target is given by:

*<sup>C</sup>* = (*ρ*<sup>2</sup> − *<sup>ρ</sup>*1)*<sup>e</sup>*

100 200 300

400 500 600

**Visibility distance [m]**

where *<sup>λ</sup>* <sup>=</sup> *<sup>d</sup>*<sup>1</sup> <sup>−</sup> *dn* <sup>1</sup> *<sup>v</sup>*1−*vh* <sup>−</sup> <sup>1</sup> *vn*−*vh* . The quadratic error is plotted in Fig. 3(b) with respect to the number of pairs of points taken into account in the calibration process. Three pairs of points, i.e. four points on the ground, suffice to obtain a quadratic error which is smaller than 2% at *d* = 136 m. This error is in line with the theoretical error at the same distance and hence acceptable.

#### *4.5.2. Implementation of the System*

The complete image acquisition system has been installed in a van equipped with a pneumatic pole. Using this vehicle, we have grabbed a sequence of images during sunrise. Sample pictures of this fog episode are shown in Fig. 4(a). The visibility distance has been estimated and is plotted in Fig. 4(b) with respect to the time. As one can see, the visibility distance increases as well as the global illumination in the scene while fog dissipates. The behavior of the system is good except at time *t* ≈ 50 min, where the visibility is underestimated. This underestimation is due to the fact that the exposure time is momentarily too high so that the images are overexposed and contrasts deteriorated. Fortunately, the auto-exposure algorithm quickly solves the problem and the visibility increases again up to the maximum value, which is above 400 m as expected.

<sup>96</sup> Climate Change and Regional/Local Responses Measurements and Observations of Meteorological Visibility at ITS Stations 9 10.5772/55697 Measurements and Observations of Meteorological Visibility at ITS Stations http://dx.doi.org/10.5772/55697 97

**Figure 4.** Implementation of the designed camera-based visibility metering system: (a) sample images acquired during a fog episode; (b) visibility distance estimated during this fog episode.

#### **5. Monitoring of meteorological visibility**

With the second system, we are able to monitor a whole spectrum of visibility ranges (from 0 to 10,000 m). In this system, we calibrated a response function of the contrast within the image with respect to reference visibility measurements obtained by external sensors. The camera needs not be a high resolution one. However, the calibration is more complex and needs at least one fog episode.

#### **5.1. Contrast of a distant target**

8 Climate Change and Regional/Local Responses

d=16 m

respect to the number of considered pairs of points.

*4.5.1. Verification of Camera Specifications*

<sup>1</sup> *<sup>v</sup>*1−*vh* <sup>−</sup> <sup>1</sup> *vn*−*vh*

*4.5.2. Implementation of the System*

**4.5. Experimental validation**

where *<sup>λ</sup>* <sup>=</sup> *<sup>d</sup>*<sup>1</sup> <sup>−</sup> *dn*

acceptable.

20 m 20 m

d=136 m

(a) (b)

**Figure 3.** Experimental verification of camera specifications: (a) experimental setup; (b) quadratic error of calibration with

First, an experimental verification of the camera specifications has been carried out to check if we are able to reach the specifications. To achieve this aim, seven cones were set on a flat road section following the experimental setup of Fig. 3(a). Using the perspective projection

> *n* ∑ *i*=1

number of pairs of points taken into account in the calibration process. Three pairs of points, i.e. four points on the ground, suffice to obtain a quadratic error which is smaller than 2% at *d* = 136 m. This error is in line with the theoretical error at the same distance and hence

The complete image acquisition system has been installed in a van equipped with a pneumatic pole. Using this vehicle, we have grabbed a sequence of images during sunrise. Sample pictures of this fog episode are shown in Fig. 4(a). The visibility distance has been estimated and is plotted in Fig. 4(b) with respect to the time. As one can see, the visibility distance increases as well as the global illumination in the scene while fog dissipates. The behavior of the system is good except at time *t* ≈ 50 min, where the visibility is underestimated. This underestimation is due to the fact that the exposure time is momentarily too high so that the images are overexposed and contrasts deteriorated. Fortunately, the auto-exposure algorithm quickly solves the problem and the visibility

*di* <sup>−</sup> *<sup>λ</sup> vi* − *vh*

. The quadratic error is plotted in Fig. 3(b) with respect to the

2

(15)

model and the positions of the different cones, we are able to calibrate the camera:

*n*=1..7

increases again up to the maximum value, which is above 400 m as expected.

(*λ*, *vh*) = argmin

Quadratic error [%]

1 2 3 4 5 6

Number of pairs of points

Let us consider an outdoor scene where targets are distributed continuously at increasing distances from the camera. When we assume that the surface of the targets is Lambertian, the luminance *L* at each point *i* of the target is given by:

$$L = \rho\_i \frac{E}{\pi} \tag{16}$$

where *E* denotes the global illumination and *ρ<sup>i</sup>* denotes the albedo at *i*. Moreover, it is a classical assumption to set *<sup>L</sup>*<sup>∞</sup> = *<sup>E</sup> π* so that the contrast of two Lambertian targets at distance *d* becomes [19]:

$$\mathcal{C} = (\rho\_2 - \rho\_1)e^{-\beta d} \approx (\rho\_2 - \rho\_1)e^{-\frac{\mathcal{M}}{\mathcal{V}}} = \Delta \rho e^{-\frac{\mathcal{M}}{\mathcal{V}}} \tag{17}$$

Consequently, the contrast of a distant Lambertian target only depends on its physical properties and on its distance to the sensor and on the meteorological visibility distance, and no longer on the illumination. Such targets allow for computing contrasts in the scene in a way which is robust to strong variations in illumination [19].

## **5.2. Probabilistic modelling**

Let us denote *φ* the probability density function (p.d.f.) of observing a contrast *C* in the scene:

$$\mathbb{P}(\mathbb{C} < X \le \mathbb{C} + \text{d}\mathbb{C}) = \phi(\mathbb{C})\text{d}\mathbb{C} \tag{18}$$

10.5772/55697

99

http://dx.doi.org/10.5772/55697

Measurements and Observations of Meteorological Visibility at ITS Stations

(22)

0 5000 10000 15000

Meteorological visibility distance [m]

**Figure 5.** Blue: curves depicting the mean contrast in random scenes with respect to the meteorological visibility distance.

the faster the curves reach a 0.5 contrast. We thus define an indicator *τ* of the system quality which is the meteorological visibility distance at which two thirds of the "capacitance" is reached. A high value of *τ* also means a lower sensitivity of the model at low meteorological

In the previous section, we have computed an analytical expression of the mean contrast expectation *mu* with respect to the meteorological visibility distance *V*. Ultimately, we would like to compute *V* as a function of *mu*. To achieve this aim, we need to invert the mean contrast expectation function (21). The inversion of this model exists and is expressed by

*<sup>V</sup>*(*mu*) = <sup>3</sup>*mud*max

1 + *muW*

where the Lambert *W* function is a transcendental function defined by solutions of

In this section, we present an experimental evaluation of the proposed model for visibility

The observation field test we used is equipped with a reference transmissometer (Degreane Horizon TI8510). It serves to calibrate different scatterometers (Degreane Horizon DF320)

estimation. To achieve this aim, we have collected ground truth data.

 *e* −1 *mu mu* 

0

**5.4. Model inversion and error estimation**

0.1

0.2

0.3

0.4

Mean contrast

Red: expectation of the mean contrast.

visibility distances.

*W*(*x*)*eW*(*x*) = *x* [20].

*5.5.1. Methodology*

5.5.1.1. Instrumentation

**5.5. Experimental evaluation**

[19]:

0.5

0.6

0.7

The expectation of the contrast *m* in the image is expressed as [19]:

$$m = \mathbb{E}[\mathbb{C}] = \int\_0^1 \mathbb{C}\phi(\mathbb{C})d\mathbb{C} \tag{19}$$

Based on (17), *C* is a random variable which depends of the two random variables *d* and ∆*ρ*. These two variables are assumed to be independent, which allows expressing Eq. (19) as [19]:

$$m = \mathbb{E}\left[\Delta\rho\right]\mathbb{E}\left[e^{-\frac{\Delta t}{\Psi}}\right] = \overline{\Delta\rho}\int\_0^{+\infty} \psi(d)e^{-\frac{\Delta t}{\Psi}}\mathbf{d}d\tag{20}$$

where ∆*ρ* denotes the mean albedo difference between the objects in the scene and *ψ* denotes the p.d.f. of there being an object at the distance *d* in the scene. Choosing a suitable target distribution *ψ* allows us computing the expectation of the contrast using Eq. (20) with respect to the meteorological visibility distance *V*.

#### **5.3. Expectation of the mean contrast**

In this paragraph, we seek an analytical expression of Eq. (20). To achieve this aim, we assume a scene which contains *n* Lambertian targets with random albedos located at random distances between 0 and *d*max. For a given sample scene, we can compute the mean contrast of the targets with respect to the meteorological visibility distance and plot the corresponding curve. Sample curves are plotted in blue in Fig. 5 (*n* = 100 and *d*max = 1000 m). We can compute the mathematical expectation of the mean contrast and obtain the following analytical model:

$$m\_{\rm ll} = \frac{V \Delta \overline{\rho}}{6d\_{\rm max}} \left[ 1 - \exp\left( -\frac{\Im d\_{\rm max}}{V} \right) \right] \tag{21}$$

where ∆*ρ*¯ is the mean albedo difference of the targets in the scene. We plot this model in red in Fig. 5. When we do not have any a priori on the targets distribution in the scene, this analytical model is the most probable with which to fit the data [19]. This fact is experimentally assessed in section 5.5. At this stage, we can make a comparison with the charging/discharging of a capacitor. The capacitance of the system is determined by the distribution of Lambertian targets in the scene. The smaller the capacitance of the system is,

**Figure 5.** Blue: curves depicting the mean contrast in random scenes with respect to the meteorological visibility distance. Red: expectation of the mean contrast.

the faster the curves reach a 0.5 contrast. We thus define an indicator *τ* of the system quality which is the meteorological visibility distance at which two thirds of the "capacitance" is reached. A high value of *τ* also means a lower sensitivity of the model at low meteorological visibility distances.

## **5.4. Model inversion and error estimation**

In the previous section, we have computed an analytical expression of the mean contrast expectation *mu* with respect to the meteorological visibility distance *V*. Ultimately, we would like to compute *V* as a function of *mu*. To achieve this aim, we need to invert the mean contrast expectation function (21). The inversion of this model exists and is expressed by [19]:

$$V(m\_{\mu}) = \frac{3m\_{\mu}d\_{\text{max}}}{1 + m\_{\mu}\mathcal{W}\left(\frac{e^{\frac{-1}{m\_{\mu}}}}{m\_{\mu}}\right)}\tag{22}$$

where the Lambert *W* function is a transcendental function defined by solutions of *W*(*x*)*eW*(*x*) = *x* [20].

## **5.5. Experimental evaluation**

In this section, we present an experimental evaluation of the proposed model for visibility estimation. To achieve this aim, we have collected ground truth data.

#### *5.5.1. Methodology*

10 Climate Change and Regional/Local Responses

**5.2. Probabilistic modelling**

scene:

Consequently, the contrast of a distant Lambertian target only depends on its physical properties and on its distance to the sensor and on the meteorological visibility distance, and no longer on the illumination. Such targets allow for computing contrasts in the scene

Let us denote *φ* the probability density function (p.d.f.) of observing a contrast *C* in the

 <sup>1</sup> 0

Based on (17), *C* is a random variable which depends of the two random variables *d* and ∆*ρ*. These two variables are assumed to be independent, which allows expressing Eq. (19) as [19]:

where ∆*ρ* denotes the mean albedo difference between the objects in the scene and *ψ* denotes the p.d.f. of there being an object at the distance *d* in the scene. Choosing a suitable target distribution *ψ* allows us computing the expectation of the contrast using Eq. (20) with respect

In this paragraph, we seek an analytical expression of Eq. (20). To achieve this aim, we assume a scene which contains *n* Lambertian targets with random albedos located at random distances between 0 and *d*max. For a given sample scene, we can compute the mean contrast of the targets with respect to the meteorological visibility distance and plot the corresponding curve. Sample curves are plotted in blue in Fig. 5 (*n* = 100 and *d*max = 1000 m). We can compute the mathematical expectation of the mean contrast and obtain the following

> *mu* <sup>=</sup> *<sup>V</sup>*∆*ρ*¯ 6*d*max

 1 − exp

where ∆*ρ*¯ is the mean albedo difference of the targets in the scene. We plot this model in red in Fig. 5. When we do not have any a priori on the targets distribution in the scene, this analytical model is the most probable with which to fit the data [19]. This fact is experimentally assessed in section 5.5. At this stage, we can make a comparison with the charging/discharging of a capacitor. The capacitance of the system is determined by the distribution of Lambertian targets in the scene. The smaller the capacitance of the system is,

<sup>−</sup> <sup>3</sup>*d*max *V*

 <sup>+</sup><sup>∞</sup> 0

*ψ*(*d*)*e* − <sup>3</sup>*<sup>d</sup>*

**P**(*C* < *X* ≤ *C* + d*C*) = *φ*(*C*)d*C* (18)

*Cφ*(*C*)d*C* (19)

*<sup>V</sup>* d*d* (20)

(21)

in a way which is robust to strong variations in illumination [19].

The expectation of the contrast *m* in the image is expressed as [19]:

*m* = **E** ∆*ρ* **E** *e* − <sup>3</sup>*<sup>d</sup> V* = ∆*ρ*

to the meteorological visibility distance *V*.

**5.3. Expectation of the mean contrast**

analytical model:

*m* = **E**[*C*] =

#### 5.5.1.1. Instrumentation

The observation field test we used is equipped with a reference transmissometer (Degreane Horizon TI8510). It serves to calibrate different scatterometers (Degreane Horizon DF320)

10.5772/55697

101

http://dx.doi.org/10.5772/55697

0 1000 2000 3000 4000

Time [minutes]

Meteorological visibility versus time

Measurements and Observations of Meteorological Visibility at ITS Stations

(a) (b) (c)

0

0.5

1

1.5

Meteorological

(d) (e) **Figure 7.** Samples of data collected in winter 2008-2009: (a) images with strong illumination conditions and presence of shadows; (b) cloudy conditions; (c) foggy weather situation; (d) background luminance and (e) meteorological visibility distance

**Figure 8.** Mask of Lambertian surfaces on our field test: The redder the pixel, the higher the confidence that the surface is

the probability that the pixel belongs to a Lambertian surface. This technique provides an efficient way to locate some Lambertian surfaces in the scene. For our field test, the mask of Lambertian surfaces is shown in Fig.8. The redder the pixel, the more the surface is likely to

 visibility [m]

2

2.5

3

3.5 x 10 <sup>4</sup>

0 1000 2000 3000 4000

Time [minutes]

Luminance versus time

0

data collected in the field test during two days.

1000

2000

3000

Luminance

Lambertian.

be Lambertian.

 [cd.m-2]

4000

5000

6000

7000

**Figure 6.** Instrumentation of our observation field test: (a) the camera grabbing pictures of the field test;(b) the scatterometer along with the background luminancemeter.

used to monitor the meteorological visibility distance in France, one of which provided our data. They are coupled with a background luminance sensor (Degreane Horizon LU320) which monitors ambient light conditions. We have added a camera which grabs images of the field test every ten minutes. The camera is an 8-bit CCD camera (640 × 480 definition, *H*=8.3 m, *θ* = 9.8*o*, *fl* = 4 mm and *tpix* = 9 *µm*). Compared to the camera specified in section 4.4, it is thus a low cost camera which is representative of common video surveillance roadside cameras (cf. section 2.2). Fig. 6(a) shows the installed camera and its orientation with respect to the field test. Fig. 6(b) shows the scatterometer and the background luminancemeter.

#### 5.5.1.2. Data collection.

We have collected two fog events the 28th February and the 1st March 2009. The fog occurred early in the morning and lasted a few hours after sunrise. During the same days, there were strong sunny weather periods. Fig. 7 shows sample images of sunny weather and cloudy weather and foggy weather. The corresponding meteorological visibility distances and luminances are plotted in Figs. 7(d,e). Obviously, the meteorological visibility distance ranges from 100 m to 35.000 m and the luminance ranges from 0 to 6.000 cd.m<sup>−</sup>2.

We have thus collected exceptional experimental data. Indeed, we met rapidly changing weather conditions over a short period of time. The ranges of meteorological visibility distance and luminance were very large. In the literature, works are dedicated to limited ranges of visibility distances [17, 21]. For example, road safety applications are dealing with 0-400 m [17] whereas people working on environmental issues are dealing with meteorological visibility distances which are above 1000 m [21]. We are among the first to have collected data encompassing both ranges. Moreover, since the data was collected over a short period of time, we consider that the content of the scene did not change. For example, we assumed that the phenology of the trees did not change, so that the amount of texture in the scene without fog remains constant.

#### 5.5.1.3. Location of Lambertian surfaces

To estimate *mu* and thus *V*, we compute the normalized gradient only on the Lambertian surfaces of the scene as proposed in section 5.1. We thus need to locate Lambertian surfaces in the images. To achieve this aim, we compute the Pearson coefficient, denoted *P<sup>L</sup> i*,*j* , between the intensity of pixels in image series where the position of the sun changes and the value of the background luminance estimated by the luminancemeter. The closer *P<sup>L</sup> <sup>i</sup>*,*<sup>j</sup>* is to 1, the stronger <sup>100</sup> Climate Change and Regional/Local Responses Measurements and Observations of Meteorological Visibility at ITS Stations 13 10.5772/55697 Measurements and Observations of Meteorological Visibility at ITS Stations http://dx.doi.org/10.5772/55697 101

12 Climate Change and Regional/Local Responses

along with the background luminancemeter.

5.5.1.2. Data collection.

(a) (b)

**Figure 6.** Instrumentation of our observation field test: (a) the camera grabbing pictures of the field test;(b) the scatterometer

used to monitor the meteorological visibility distance in France, one of which provided our data. They are coupled with a background luminance sensor (Degreane Horizon LU320) which monitors ambient light conditions. We have added a camera which grabs images of the field test every ten minutes. The camera is an 8-bit CCD camera (640 × 480 definition, *H*=8.3 m, *θ* = 9.8*o*, *fl* = 4 mm and *tpix* = 9 *µm*). Compared to the camera specified in section 4.4, it is thus a low cost camera which is representative of common video surveillance roadside cameras (cf. section 2.2). Fig. 6(a) shows the installed camera and its orientation with respect to the field test. Fig. 6(b) shows the scatterometer and the background luminancemeter.

We have collected two fog events the 28th February and the 1st March 2009. The fog occurred early in the morning and lasted a few hours after sunrise. During the same days, there were strong sunny weather periods. Fig. 7 shows sample images of sunny weather and cloudy weather and foggy weather. The corresponding meteorological visibility distances and luminances are plotted in Figs. 7(d,e). Obviously, the meteorological visibility distance

We have thus collected exceptional experimental data. Indeed, we met rapidly changing weather conditions over a short period of time. The ranges of meteorological visibility distance and luminance were very large. In the literature, works are dedicated to limited ranges of visibility distances [17, 21]. For example, road safety applications are dealing with 0-400 m [17] whereas people working on environmental issues are dealing with meteorological visibility distances which are above 1000 m [21]. We are among the first to have collected data encompassing both ranges. Moreover, since the data was collected over a short period of time, we consider that the content of the scene did not change. For example, we assumed that the phenology of the trees did not change, so that the amount of

To estimate *mu* and thus *V*, we compute the normalized gradient only on the Lambertian surfaces of the scene as proposed in section 5.1. We thus need to locate Lambertian surfaces in

intensity of pixels in image series where the position of the sun changes and the value of the

*i*,*j*

*<sup>i</sup>*,*<sup>j</sup>* is to 1, the stronger

, between the

the images. To achieve this aim, we compute the Pearson coefficient, denoted *P<sup>L</sup>*

background luminance estimated by the luminancemeter. The closer *P<sup>L</sup>*

ranges from 100 m to 35.000 m and the luminance ranges from 0 to 6.000 cd.m<sup>−</sup>2.

texture in the scene without fog remains constant.

5.5.1.3. Location of Lambertian surfaces

**Figure 7.** Samples of data collected in winter 2008-2009: (a) images with strong illumination conditions and presence of shadows; (b) cloudy conditions; (c) foggy weather situation; (d) background luminance and (e) meteorological visibility distance data collected in the field test during two days.

**Figure 8.** Mask of Lambertian surfaces on our field test: The redder the pixel, the higher the confidence that the surface is Lambertian.

the probability that the pixel belongs to a Lambertian surface. This technique provides an efficient way to locate some Lambertian surfaces in the scene. For our field test, the mask of Lambertian surfaces is shown in Fig.8. The redder the pixel, the more the surface is likely to be Lambertian.

10.5772/55697

103

http://dx.doi.org/10.5772/55697

0 5000 10000 15000

Meteorological visibility distance [m]

We have to fit the mean contrast model (Eq. (21)) to the data shown in Fig. 9(a) using robust regression techniques. To ensure a mathematical solution, we have fitted the model (Eq. (24)), which is slightly different from the theoretical model. Three unknown variables *a*, *b* and *dmax*

<sup>1</sup> − exp

This model fits well with the data (*R*<sup>2</sup> = 0.91). In particular, we obtain *d*max = 307.2 m. The fitted curve is plotted in Fig. 10. We estimated a capacitance (as defined in section 5.3) of the

From the fitted model, we can now invert the model using (22) and estimate the

After having estimated the meteorological visibility distance, we can compute the error on this estimation. The results are given in Tab. 3. Since the applications are very different depending on the range of meteorological visibility distances, we have computed the error and the standard deviation for various applications: road safety, meteorological observation and air quality. One can see that the error remains low for critical safety applications. It

 *ae <sup>a</sup> b*−*mu b* − *mu*

 − *a*

*<sup>V</sup>*˜ <sup>=</sup> <sup>3</sup>*d*max(*<sup>b</sup>* <sup>−</sup> *mu*)

(*b* − *mu*)*W*

<sup>−</sup> <sup>3</sup>*d*max *V*

Measurements and Observations of Meteorological Visibility at ITS Stations

+ *b* (24)

(25)

have to be estimated, which can be easily done using classical curve fitting tools.

**Figure 10.** Data fitting with the mean contrast model. Dots: data. Red curve: fitted model.

*<sup>m</sup>*˜*<sup>u</sup>* <sup>=</sup> *aV <sup>d</sup>*max

meteorological visibility distance *V*˜ based on the mean contrast *mu*[19]:

Weighted sum of contrast

5.5.2.2. Model fitting

system *τ* ≈ 950 m. 5.5.2.3. Discussions

**Figure 9.** Visibility estimators: (a) the estimator is based on the contrast on Sobel gradient alone; (b) the estimator is based on Sobel gradient weighted by Lambertian surfaces.

Having located the Lambertian surfaces, we can compute the gradients in the scene by means of the module of the Sobel filter. For each pixel, we normalize the gradient *Gi*,*<sup>j</sup>* by the intensity of the background. Since our camera is equipped with an auto-iris, the background intensity *A*<sup>∞</sup> is most of the time equal to 28 − 1, so that this step can be avoided. Each gradient is then weighted by *P<sup>L</sup> i*,*j* , the probability of a pixel to belong to a Lambertian surface where no depth discontinuity exists (*P<sup>L</sup>* is mostly very small). Consequently, only relevant areas of the image are used, and the scene need not be totally Lambertian. Finally, the estimated contrast in the scene *m*˜*u* is given by [19]:

$$\mathfrak{M}\_{\mathcal{U}} = \sum\_{i=0}^{h} \sum\_{j=0}^{w} \Delta \rho\_{i,j} \exp\left(-\frac{\mathfrak{J}d\_{i,j}}{V}\right) \approx \sum\_{i=0}^{h} \sum\_{j=0}^{w} \frac{G\_{i,j}}{A\_{\infty}} P\_{i,j}^{L} \tag{23}$$

where <sup>∆</sup>*ρi*,*<sup>j</sup>* is the intrinsic contrast of a pixel in Eq. (17), and *<sup>h</sup>* and *<sup>w</sup>* are respectively the height and the width of the images.

#### *5.5.2. Results*

#### 5.5.2.1. Contrast estimators

We have computed Eq. (23) for our collection of 150 images with different meteorological visibility distances. For comparison purposes, we have also computed the simple sum of gradients in the image without weighting the Lambertian surfaces. The results are shown in Fig. 9. By using the Lambertian surfaces, we can see that the shape of the distribution in Fig. 9(a) looks like the curve proposed in Fig. 5, which is very satisfactory. Conversely, when all the pixels of the scene are used, the points are more scattered when the meteorological visibility distance is above 2500 m (see Fig. 9(b)): When the sky is clear and the visibility is high, the illumination from the sun strongly influences the gradients in the scene. Consequently, the estimation of the visibility is altered. These two distributions show the benefit of selecting the Lambertian surfaces to estimate the visibility distance.

10.5772/55697

**Figure 10.** Data fitting with the mean contrast model. Dots: data. Red curve: fitted model.

#### 5.5.2.2. Model fitting

14 Climate Change and Regional/Local Responses

Sum of gradients <sup>6</sup> x 104

0 5000 10000 15000

*i*,*j*

*h* ∑ *i*=0

*w* ∑ *j*=0

<sup>∆</sup>*ρi*,*<sup>j</sup>* exp

0 5000 10000 15000

Meteorological visibility distance (m)

, the probability of a pixel to belong to a Lambertian surface

*w* ∑ *j*=0 *Gi*,*<sup>j</sup> A*∞ *PL*

*<sup>i</sup>*,*<sup>j</sup>* (23)

(a)

1.4 1.6 1.8 2 2.2 2.4 2.6

Sum of contrasts taken on Lambertian

**Figure 9.** Visibility estimators: (a) the estimator is based on the contrast on Sobel gradient alone; (b) the estimator is based on

Having located the Lambertian surfaces, we can compute the gradients in the scene by means of the module of the Sobel filter. For each pixel, we normalize the gradient *Gi*,*<sup>j</sup>* by the intensity of the background. Since our camera is equipped with an auto-iris, the background intensity *A*<sup>∞</sup> is most of the time equal to 28 − 1, so that this step can be avoided. Each

where no depth discontinuity exists (*P<sup>L</sup>* is mostly very small). Consequently, only relevant areas of the image are used, and the scene need not be totally Lambertian. Finally, the

> <sup>−</sup>3*di*,*<sup>j</sup> V* ≈ *h* ∑ *i*=0

where <sup>∆</sup>*ρi*,*<sup>j</sup>* is the intrinsic contrast of a pixel in Eq. (17), and *<sup>h</sup>* and *<sup>w</sup>* are respectively the

We have computed Eq. (23) for our collection of 150 images with different meteorological visibility distances. For comparison purposes, we have also computed the simple sum of gradients in the image without weighting the Lambertian surfaces. The results are shown in Fig. 9. By using the Lambertian surfaces, we can see that the shape of the distribution in Fig. 9(a) looks like the curve proposed in Fig. 5, which is very satisfactory. Conversely, when all the pixels of the scene are used, the points are more scattered when the meteorological visibility distance is above 2500 m (see Fig. 9(b)): When the sky is clear and the visibility is high, the illumination from the sun strongly influences the gradients in the scene. Consequently, the estimation of the visibility is altered. These two distributions show

the benefit of selecting the Lambertian surfaces to estimate the visibility distance.

2.8 x 104

 surfaces

Meteorological visibility distance (m)

estimated contrast in the scene *m*˜*u* is given by [19]:

*m*˜*<sup>u</sup>* =

Sobel gradient weighted by Lambertian surfaces.

gradient is then weighted by *P<sup>L</sup>*

height and the width of the images.

5.5.2.1. Contrast estimators

*5.5.2. Results*

(b)

We have to fit the mean contrast model (Eq. (21)) to the data shown in Fig. 9(a) using robust regression techniques. To ensure a mathematical solution, we have fitted the model (Eq. (24)), which is slightly different from the theoretical model. Three unknown variables *a*, *b* and *dmax* have to be estimated, which can be easily done using classical curve fitting tools.

$$\tilde{m}\_{\rm u} = \frac{aV}{d\_{\rm max}} \left[ 1 - \exp\left( -\frac{\Im d\_{\rm max}}{V} \right) \right] + b \tag{24}$$

This model fits well with the data (*R*<sup>2</sup> = 0.91). In particular, we obtain *d*max = 307.2 m. The fitted curve is plotted in Fig. 10. We estimated a capacitance (as defined in section 5.3) of the system *τ* ≈ 950 m.

#### 5.5.2.3. Discussions

From the fitted model, we can now invert the model using (22) and estimate the meteorological visibility distance *V*˜ based on the mean contrast *mu*[19]:

$$\tilde{V} = \frac{3d\_{\text{max}}(b - m\_{\text{ul}})}{(b - m\_{\text{ul}})W\left(\frac{ae^{\frac{d}{b - m\_{\text{ul}}}}}{b - m\_{\text{ul}}}\right) - a} \tag{25}$$

After having estimated the meteorological visibility distance, we can compute the error on this estimation. The results are given in Tab. 3. Since the applications are very different depending on the range of meteorological visibility distances, we have computed the error and the standard deviation for various applications: road safety, meteorological observation and air quality. One can see that the error remains low for critical safety applications. It increases for higher ranges of visibility distances, and becomes huge for visibility distances above 7 km. Different issues may be discussed. First, the model presented in this chapter is relevant for uniform distributions of distances which happen in many environments, such as highway scenes. The scene in which the experimental data used in this paper were collected may be meet this assumption. Second, the Sobel operator is certainly not the best estimate for the gradient. Indeed, it is a simple high-pass filter which is problematic because of the impulse noise of camera sensors. Different filters may be used to enhance the images beforehand, or to compute the contrast more robustly.

10.5772/55697

105

http://dx.doi.org/10.5772/55697

**7.2. Fog nowcasting**

The forecasting of the weather within the next six hours is often referred to as nowcasting. In this time range, it is possible to forecast smaller features such as individual showers and thunderstorms with reasonable accuracy, as well as other features too small to be resolved by a computer model. [25] show that combining satellite-based forecasting of low clouds with terrestrial measurements of humidity allows computing a probability of fog occurrence. A camera-based visibility meter could easily substitute for a humidity sensor. Indeed, using camera-based visibility estimation and meteorological data, [26] showed that visibility can be predicted up to 15 minutes in advance with 1-km mesh meteorological data. Such camera-based nowcasting methods may be good solutions to allow the re-routing of vehicles

[27] has proposed a review of best practices in terms of mitigation of highway visibility problems, in particular fog related issues. In this paper, he describes existing installations in the USA dedicated to driver alert in case of low visibility on the highway. Apart from fog dispersal techniques, the best practices are related to the timely alert of the drivers which approach a foggy area. Then, depending on the fog density, different advisory speed limits may be posted. In the same time, the public lighting is adapted. The component of these systems are made of weather stations, CCTV cameras and Variable Message Signs (VMS). More recently, Caltrans has installed the "Fog Pilot" system in District 6, which provides a

This centralized solution relies on the use of infrastructure to vehicle communications to warn the drivers whose vehicles are equipped with a receiver, in case of sudden low speed area. Thanks to the proposed camera-based visibility monitoring techniques, we are able to build a decentralized fog-pilot, which makes use of CCTV cameras to monitor the visibility and allows optimizing the speed of drivers approaching a low visibility area, as well as the intensity of road studs. Based on best existing practices, its principle is to warn the drivers of a foggy area with enough time, so that they can adapt their speed to the prevailing visibility

Reliable solutions to accurately monitor the meteorological visibility along road networks at reasonable costs are still not available. The use of the cameras, which are multifunctional sensors and are already deployed along the roadsides, is a promising solution. However, progress in computer vision are still needed to obtain robust techniques, which are able to fulfill the needs of transportation safety, of meteorology and of environment in terms of

In this chapter, we have presented two different camera-based systems to estimate the visibility distance. These systems could be integrated in ITS stations to alert drivers as well public authorities in case of fog hazard. The first system is dedicated to safety-critical visibility ranges but needs a high-resolution camera along with a simple calibration process. The second technique is not restricted to low visibility ranges and needs a low cost camera but a more complex calibration procedure. In the future, it is obvious that the

<sup>4</sup> mile along a 12-mile (20-km) portion of State Route 99.

Measurements and Observations of Meteorological Visibility at ITS Stations

before they reach a low-visibility area in a timely manner.

**7.3. Pile-up prevention and mitigation**

high-technology solution every <sup>1</sup>

distance in the dangerous area.

observation.

**8. Conclusion and perspectives**

## **6. Outcome**

First, a camera-based system has been developed to detect safety-critical visibility conditions. Sample results based on data collected during a morning fog episode are shown in the Fig. 4. When visibility is below 400 m, the accuracy of the system is expected to be 90%. The second system has been assessed on the test site of Météo France in Trappes, where images have been grabbed with visibility data (see Fig. 7). Based on our experiments, we are able to obtained error estimates which are lower than 10% [19] as well.

We thus have at disposal two different techniques to estimate the visibility distance in case of fog or haze. The first technique is dedicated to low visibility ranges but needs a high-resolution camera along with a simple calibration process. The second technique is not restricted to low visibility ranges and works with a low cost camera but needs a more complex calibration procedure. It is clear that the complementarity between both approaches must be studied so as to build a single system, which is easy to set up and deploy.

## **7. Potential applications**

## **7.1. Winter maintenance**

The monitoring of meteorological visibility has different applications for winter maintenance [22]. First, the knowledge of a low visibility area is important for the safety of winter maintenance operations. Second, a sudden drop of the visibility can be due to heavy snow falls. The relationship between liquid equivalent snowfall rate and visibility has also been investigated [23], which means that a camera-based visibility meter is potentially a good snow sensor. Finally, meteorological models have been developed to forecast pavement temperatures as well as snow height [24]. Nebulosity and fog are phenomena which alter the prediction, since the radiative transfer between the pavement and the air is affected. The assimilation of visibility data in these prediction models may be useful to increase the accuracy of forecasts.


**Table 3.** Relative errors of meteorological visibility distance estimation with respect to the envisioned application.

## **7.2. Fog nowcasting**

16 Climate Change and Regional/Local Responses

**6. Outcome**

**7. Potential applications**

**7.1. Winter maintenance**

accuracy of forecasts.

beforehand, or to compute the contrast more robustly.

obtained error estimates which are lower than 10% [19] as well.

increases for higher ranges of visibility distances, and becomes huge for visibility distances above 7 km. Different issues may be discussed. First, the model presented in this chapter is relevant for uniform distributions of distances which happen in many environments, such as highway scenes. The scene in which the experimental data used in this paper were collected may be meet this assumption. Second, the Sobel operator is certainly not the best estimate for the gradient. Indeed, it is a simple high-pass filter which is problematic because of the impulse noise of camera sensors. Different filters may be used to enhance the images

First, a camera-based system has been developed to detect safety-critical visibility conditions. Sample results based on data collected during a morning fog episode are shown in the Fig. 4. When visibility is below 400 m, the accuracy of the system is expected to be 90%. The second system has been assessed on the test site of Météo France in Trappes, where images have been grabbed with visibility data (see Fig. 7). Based on our experiments, we are able to

We thus have at disposal two different techniques to estimate the visibility distance in case of fog or haze. The first technique is dedicated to low visibility ranges but needs a high-resolution camera along with a simple calibration process. The second technique is not restricted to low visibility ranges and works with a low cost camera but needs a more complex calibration procedure. It is clear that the complementarity between both approaches

The monitoring of meteorological visibility has different applications for winter maintenance [22]. First, the knowledge of a low visibility area is important for the safety of winter maintenance operations. Second, a sudden drop of the visibility can be due to heavy snow falls. The relationship between liquid equivalent snowfall rate and visibility has also been investigated [23], which means that a camera-based visibility meter is potentially a good snow sensor. Finally, meteorological models have been developed to forecast pavement temperatures as well as snow height [24]. Nebulosity and fog are phenomena which alter the prediction, since the radiative transfer between the pavement and the air is affected. The assimilation of visibility data in these prediction models may be useful to increase the

**Application Highway fog Meteorological fog Haze Air quality** Range [m] 0-400 0-1000 0-5000 0-15000 Number of data 13 19 45 150 Mean error [%] 12.6 18.1 29.7 - Std [%] 13.7 18.9 22 -

**Table 3.** Relative errors of meteorological visibility distance estimation with respect to the envisioned application.

must be studied so as to build a single system, which is easy to set up and deploy.

The forecasting of the weather within the next six hours is often referred to as nowcasting. In this time range, it is possible to forecast smaller features such as individual showers and thunderstorms with reasonable accuracy, as well as other features too small to be resolved by a computer model. [25] show that combining satellite-based forecasting of low clouds with terrestrial measurements of humidity allows computing a probability of fog occurrence. A camera-based visibility meter could easily substitute for a humidity sensor. Indeed, using camera-based visibility estimation and meteorological data, [26] showed that visibility can be predicted up to 15 minutes in advance with 1-km mesh meteorological data. Such camera-based nowcasting methods may be good solutions to allow the re-routing of vehicles before they reach a low-visibility area in a timely manner.

## **7.3. Pile-up prevention and mitigation**

[27] has proposed a review of best practices in terms of mitigation of highway visibility problems, in particular fog related issues. In this paper, he describes existing installations in the USA dedicated to driver alert in case of low visibility on the highway. Apart from fog dispersal techniques, the best practices are related to the timely alert of the drivers which approach a foggy area. Then, depending on the fog density, different advisory speed limits may be posted. In the same time, the public lighting is adapted. The component of these systems are made of weather stations, CCTV cameras and Variable Message Signs (VMS). More recently, Caltrans has installed the "Fog Pilot" system in District 6, which provides a high-technology solution every <sup>1</sup> <sup>4</sup> mile along a 12-mile (20-km) portion of State Route 99. This centralized solution relies on the use of infrastructure to vehicle communications to warn the drivers whose vehicles are equipped with a receiver, in case of sudden low speed area. Thanks to the proposed camera-based visibility monitoring techniques, we are able to build a decentralized fog-pilot, which makes use of CCTV cameras to monitor the visibility and allows optimizing the speed of drivers approaching a low visibility area, as well as the intensity of road studs. Based on best existing practices, its principle is to warn the drivers of a foggy area with enough time, so that they can adapt their speed to the prevailing visibility distance in the dangerous area.

## **8. Conclusion and perspectives**

Reliable solutions to accurately monitor the meteorological visibility along road networks at reasonable costs are still not available. The use of the cameras, which are multifunctional sensors and are already deployed along the roadsides, is a promising solution. However, progress in computer vision are still needed to obtain robust techniques, which are able to fulfill the needs of transportation safety, of meteorology and of environment in terms of observation.

In this chapter, we have presented two different camera-based systems to estimate the visibility distance. These systems could be integrated in ITS stations to alert drivers as well public authorities in case of fog hazard. The first system is dedicated to safety-critical visibility ranges but needs a high-resolution camera along with a simple calibration process. The second technique is not restricted to low visibility ranges and needs a low cost camera but a more complex calibration procedure. In the future, it is obvious that the complementarity between both approaches must be studied so as to build a single system, which will be easy to set up and to deploy. Second, the problem of night fog is still subject to research. In particular, we would like to adapt our in-vehicle techniques of night fog detection to visual surveillance. Third, we would like to develop an intelligent camera where fog detection is implemented along with existing traffic applications (virtual loops, AID, etc.). Finally, we would like to develop test beds for our sensing technologies and our decentralized fog-pilot.

10.5772/55697

107

http://dx.doi.org/10.5772/55697

[9] R. Hallowell, M. Matthews, and P. Pisano. An automated visibility detection algorithm utilizing camera imagery. In *23rd Conference on Interactive Information and Processing Systems for Meteorology, Oceanography, and Hydrology (IIPS), San Antonio, TX, Amer.*

Measurements and Observations of Meteorological Visibility at ITS Stations

[10] S. Roquelaure, R. Tardif, S. Remy, and T Bergot. Skill of a ceiling and visibility local ensemble prediction system (leps) according to fog-type prediction at paris-charles de

[12] WMO. *Guide to Meteorological Instruments and Methods of Observation*. Number 8. World

[13] AFNOR. Road meteorology - gathering of meteorological and road data - terminology.

[15] N. Hautière, J.-P. Tarel, J. Lavenant, and D. Aubert. Automatic fog detection and estimation of visibility distance through use of an onboard camera. *Machine Vision*

[16] N. Hautière, R. Labayrade, and D. Aubert. Real-time disparity contrast combination for onboard estimation of the visibility distance. *IEEE Transactions on Intelligent*

[17] N. Hautière, E. Bigorgne, and D. Aubert. Visibility range monitoring through use of a

[18] J.D. Crosby. Visibility sensor accuracy: what's realistic? In *12th Symposium on*

[19] N. Hautière, R. Babari, E. Dumont, R. Brémond, and N. Paparoditis. *Lecture Notes in Computer Science, Computer Vision - ACCV 2010*, volume 6495, chapter Estimating Meteorological Visibility using Cameras: A Probabilistic Model-Driven Approach, pages

[20] R. M Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth. On the Lambert

[21] C.-H. Luo, C.-Y. Wen, C.-S. Yuan, C.-C. Liaw, J.-L. ans Lo, and S.-H. Chiu. Investigation of urban atmospheric visibility by high-frequency extraction: Model development and

[22] Y. Nagata, T. Hagiwara, K. Araki, Y. Kaneda, and H. Sasaki. Application of road visibility information system to winter maintenance. *Transportation Research Records:*

W function. *Advances in Computational Mathematics*, 5:329–359, 1996.

[14] W.E.K. Middleton. *Vision through the atmosphere*. University of Toronto Press, 1952.

gaulle. *Airport. Weather and Forecasting*, 24:1511–1523, 2009.

[11] CIE. *International Lighting Vocabulary*. Number 17.4. 1987.

Meteorological Organization, 2008.

NF P 99-320, April 1998.

*Applications*, 17(1):8–20, 2006.

243–254. Springer, March 2011.

*Transportation Systems*, 7(2):201–212, June 2006.

*Meteorological Observations and Instrumentation*, 2003.

field test. *Atmospheric Environment*, 39:2545–2552, 2005.

*Journal of the Transportation Research Board*, 2055:128–138, 2008.

roadside camera. In *IEEE Intelligent Vehicles Symposium*, 2008.

*Meteor. Soc.*, 2007.

## **Author details**

Nicolas Hautière1, Raouf Babari1, Eric Dumont1, Jacques Parent Du Chatelet2 and Nicolas Paparoditis<sup>3</sup>

1 Université Paris-Est, Institut Français des Sciences et Technologies des Transports, de l'Aménagement et des Réseaux, France 2 Météo France, France

3 Université Paris-Est, Institut Géographique National, France

## **References**


10.5772/55697


18 Climate Change and Regional/Local Responses

Nicolas Hautière1, Raouf Babari1, Eric Dumont1, Jacques Parent Du Chatelet2 and Nicolas Paparoditis<sup>3</sup>

3 Université Paris-Est, Institut Géographique National, France

l'Aménagement et des Réseaux, France

fog-pilot.

**Author details**

2 Météo France, France

43:1730–1737, 2011.

(2069):9–15, 2008.

*Board*, (1937):87–95, 2005.

(2139):97–106., 2009.

Swansea, Wales, United Kingdom, 1991.

**References**

complementarity between both approaches must be studied so as to build a single system, which will be easy to set up and to deploy. Second, the problem of night fog is still subject to research. In particular, we would like to adapt our in-vehicle techniques of night fog detection to visual surveillance. Third, we would like to develop an intelligent camera where fog detection is implemented along with existing traffic applications (virtual loops, AID, etc.). Finally, we would like to develop test beds for our sensing technologies and our decentralized

1 Université Paris-Est, Institut Français des Sciences et Technologies des Transports, de

[1] M. Abdel-Aty, A.-A. Ekram, H. Huang, and K. Choi. A study on crashes related to visibility obstruction due to fog and smoke. *Accident Analysis and Prevention*,

[2] B. Whiffen, P. Delannoy, and S.. Siok. Fog: Impact on road transportation and mitigation options. In *National Highway Visibility Conference, Madison, Wisconsin, USA*, May 2004.

[3] J. Kang, R. Ni, and G. J. Andersen. Effects of reduced visibility from fog on car-following performance. *Transportation Research Record: Journal of the Transportation Research Board*,

[5] C. A. Mac Carley. Methods and metrics for evaluation of an automated real-time driver warning system. *Transportation Research Record: Journal of the Transportation Research*

[6] T. Pejovic, V. A. Williams, R. B. Noland, and R. Toumi. Factors affecting the frequency and severity of airport weather delays and the implications of climate change for future delays. *Transportation Research Record: Journal of the Transportation Research Board*,

[7] A. H. Perry and L. J. Symons. *Highway Meteorology*. University of Wales Swansea,

[8] N. Jacobs, Burgin W., N. Fridrich, A. Abrams, K. Miskell, B. Brswell, A. Richardson, and R. Pless. The global network of outdoor webcams: Properties and apllications. In *ACM*

*International Conference on Advances in Geographic Information Systems*, 2009.

[4] F.D. Shepard. *Reduced Visibility Due to Fog on the Highway*. Number 228. 1996.


[23] Jothiram Vivekanandan Jeffrey Cole Barry Myers Charles Masters Rasmussen, Roy M. The estimation of snowfall rate using visibility. *Journal of Applied Meteorology and Climatology*, 38(10):1542–1563, 1999.

**Section 2**

**Climate Modeling**


**Section 2**
