Deriving turbulence characteristics from the COSMO numerical weather prediction model for dispersion applications

At MeteoSwiss an integrated modelling system is used to simulate the dispersion of radioactive material in emergency situations. For the prediction of the atmospheric flow, the COSMO numerical weather prediction model is used. The model is run operationally at 6.6 and 2.2 km horizontal resolution, respectively and uses a 1.5 order turbulence closure with a prognostic equation for turbulent kinetic energy. Both versions of the COSMO model are coupled o ff-line with a Lagrangian particle dispersion model (LPDM). The aim of this study is to investigate the sensitivity of the dispersion model to di fferent interfacing approaches between LPDM and the COSMO model. The diagnosed turbulence variables are validated on an ideal convective case and two measurement campaigns. Simulations of hypothetical pollutant releases show that the di fferent interfacing approaches can lead to substantial changes in the forecasted concentrations.


Introduction
Lagrangian particle dispersion models are among the most sophisticated tools to simulate atmospheric dispersion of pollutants.For this type of model the pollutant cloud is simulated by a large number (more than 100 000) of individual particles.To be able to calculate the trajectories of each particle, information from the mean atmospheric variables as well as from the turbulence state of the atmosphere is required.In most operational systems, the mean meteorological variables can directly be extracted from a numerical weather prediction (NWP) model but turbulence characteristics have to be parameterized.This is done by using a meteorological preprocessor or interface.In the present study the sensitivity of a dispersion model to different interfacing approaches will be presented.First, the models used operationally at Me-teoSwiss and two different diagnostic methods for turbulence variables are described.Second, the validation of diagnosed turbulence characteristics is presented and finally, the impact on dispersion characteristics is studied on a hypothetical pollutant release.

Methodology
At MeteoSwiss an integrated modelling system is used to simulate the dispersion of radioactive material in emergency situations.In this system the COSMO numerical weather prediction model is used for the prediction of the atmospheric flow.The COSMO model is a limited-area numerical weather prediction model (Doms and Schaettler, 2002) which is being developed in the framework of the COSMO consortium (COnsortium for Small-scale MOdelling).At MeteoSwiss the COSMO model is run operationally at two horizontal resolutions.COSMO-7 has a horizontal resolution of 6.6 km and is integrated for 72 h twice a day on a European domain.COSMO-2 has a 2.2 km horizontal resolution and provides 24 h forecasts eight times a day for a smaller domain covering Switzerland.For the parameterization of atmospheric turbulence the COSMO model uses a one-and-a-half order closure (Buzzi et al., 2009), which corresponds to level 2.5 in the Mellor and Yamada notation (Mellor and Yamada, 1982).This closure type carries a prognostic equation for turbulent kinetic energy (TKE).The COSMO model predicts all the meteorological parameters (e.g.wind and temperature profiles) which are relevant for dispersion modelling with high accuracy.At MeteoSwiss the COSMO model is continuously verified against radio soundings (Arpagaus, 2005) and surface observations (Kaufmann, 2005).
Published by Copernicus Publications.

B. Szintai et al.: Deriving turbulence characteristics from the COSMO model
The model used for the calculation of pollutant dispersion in an emergency situation is the Lagrangian Particle Dispersion Model (LPDM), which was developed by Glaab et al. (1996) at the German Weather Service (DWD).In LPDM the trajectory of each particle is calculated using the actual wind speed at the position of the particle, which is decomposed into a mean and a turbulent component.The mean wind component is taken directly from the COSMO model, while the turbulent component is computed using the Langevin equation (Legg and Raupach, 1982).Evaluation of this equation requires the following turbulence variables at the particle's position: autocorrelation function and the Lagrangian timescale (parameterized according to Taylor (1921), and the standard deviation of the wind fluctuations (σ k ).The latter are derived from the turbulent kinetic energy (e), which is taken directly from the COSMO model: where m k is the portion of TKE for the given coordinate direction.In the standard application of LPDM the vertical portion of TKE is determined according to: R f is the flux Richardson number and L c =0.052 is the ratio of the vertical and horizontal diffusion length scales.In LPDM a horizontally isotropic turbulence is assumed, which leads to: Standard deviations of velocity fluctuations derived according to this approach will be further referred to as the "Direct" method.
A different approach to diagnose the turbulence variables for a dispersion model is to apply similarity theory considerations.In this case usually the surface fluxes and a diagnosed planetary boundary layer (PBL) height is needed from the NWP model.This approach is used e.g. in the Lagrangian dispersion model FLEXPART (Stohl et al., 2005).In the FLEXPART model the turbulence characteristics are parameterized according to Hanna (1982).In this approach for the diagnosis of turbulence characteristics the boundary layer parameters h, L, w * , z 0 and u * are used, i.e. the PBL height, the Obukhov length, the convective velocity scale, the roughness length and the friction velocity, respectively.During unstable conditions the standard deviations of wind fluctuations are computed as: For stable conditions standard deviations are assumed to decrease linearly with height: Standard deviations of the wind fluctuations derived according to this approach will be further referred to as the "Hanna" method.To be able to use similarity theory approaches for the determination of dispersion parameters, first the PBL height has to be diagnosed from COSMO model outputs.
Several methods have been tested for this purpose, using, e.g., the bulk and gradient Richardson number, a TKE criterion or criteria on heat and momentum fluxes from the model.Various theoretical approaches for the evolution of both the stable and convective PBL have been investigated over the years.For the stable case the diagnostic multi-limit formulation by Zilitinkevich et al. (2007) is used in the present study, while in the convective case the prognostic slab model of Batchvarova and Gryning (1991) is applied.Results were validated with profiles from five radio sounding stations over a one-month period (Fig. 1).Overall, methods based on the bulk Richardson number and momentum fluxes of the COMSO model yield good agreement with measurements (Szintai and Kaufmann, 2008).PBL heights diagnosed form TKE profiles are considerably overestimated.For the present purpose both approaches (the bulk Ri and the TKE methods) to determine the boundary layer heights are employed in Sect. 4 as a sensitivity test.

Validation of turbulence characteristics
The diagnosed turbulence variables are validated on several measurement data sets, before applying them to the dispersion model.First, an ideally convective case is investigated, which is described in Mironov et al. (2000).The setting for this simulation was a horizontally homogeneous and flat terrain with constant heating rate at the bottom.In the simulation no phase changes were considered (dry case) and wind shear was neglected.For this case the Large Eddy Simulation (LES) data set is available that contains all the necessary turbulence characteristics, which can be compared to single column runs of the COSMO model.Figure 2 shows profiles of the standard deviations of wind fluctuations after the steady state was achieved in the simulation.In the horizontal direction the standard deviation is considerably overestimated by the "Hanna" approach throughout the whole PBL, while the "Direct" approach gives good results, especially in the middle of the PBL.As wind shear was neglected in the simulation, the along-wind and cross-wind standard deviations are identical.The vertical standard deviation is simulated well by both methods: However, in the upper PBL the "Hanna"   method overestimates while the "Direct" method underestimates the LES values.It has to be noted, that in the case of the "Hanna" approach, the bulk Richardson number was used to determine the PBL height, which show good agreement with the PBL height evaluated subjectively from the heat flux profile of the LES.
Diagnosed turbulence variables were furthermore evaluated on the LITFASS-2003 measurement campaign (Beyrich and Mengelkamp, 2006).In this case both the single column and the full three dimensional version of the COSMO model were used.In the case of the single column model, simulations were initialized with the measured soil temperature and moisture profiles and the radio sounding.In the case of the three dimensional model a long term (1 month) run was performed with the standalone version of the soil module of COSMO to produce a correct soil analysis.The single column and the 3-D runs gave similar results considering turbulence characteristics.The model results were compared to surface micrometeorological measurements and to turbulence data from a 100 m high tower.Figure 3 presents verification results of the single column model for the standard deviations of horizontal and vertical wind fluctuations at 90 m height.It can be noted, that the "Hanna" approach always predicts higher turbulence values than the "Direct" approach.In the case of horizontal wind fluctuations the measured turbulence intensity lies between the two predicted values.For the vertical fluctuations the "Direct" gives accurate values while the "Hanna" approach overestimates the measured turbulence intensity, especially during daytime.
The diagnosed turbulence characteristics were also evaluated on a MeteoSwiss measurement campaign using four sonic anemometers near Swiss nuclear power plants (CN-Met campaign).Operational forecasts of the COSMO-2 model were evaluated on a three month period between 1 August 2008 and 31 October 2008.The verification results show an overall good performance of the COSMO-2 model (Fig. 4), with all the selected turbulence parameters being in an acceptable range (20-30% relative bias).Turbulent kinetic energy, which is the only turbulence related model variable in COSMO, is generally underestimated by the model, except for the Beznau site, where the model grid points are characterized by significantly higher roughness lengths, compared to other sites.Very good performance was observed in the case of vertical turbulence, which is the most important turbulence variable with respect to mesoscale  Verified parameters: wind speed, standard deviation of horizontal and vertical wind fluctuations ("Direct" approach: sigma x, sigma y and sigma z; "Hanna" approach: sigma u, sigma v, sigma w), Turbulent Kinetic Energy (TKE).dispersion modelling.The standard deviations of horizontal wind speed are not so well predicted, as that for the vertical component.The "Hanna" method shows slightly better performance than that based on the direct use of TKE from the COSMO model.

Impact on dispersion
In order to study the impact on dispersion characteristics, the two different interfacing approaches are introduced to the emergency system applied at MeteoSwiss, and the impact is evaluated on hypothetical case studies.In the following, the case study of 8 September 2008 will be presented.As the simulated concentrations cannot be compared to measurements, we will compare the different simulations to each other and investigate the relative difference.The synoptic situation was characterized by an anticyclone over Central-Europe, which caused weak southwesterly flow over Switzerland.Due to the calm winds and the lack of precipitation, turbulence potentially plays a significant role in the dispersion of pollutants.Three different simulations were made with the COSMO-7 -LPDM system for this case.First, the "Hanna" approach was applied with PBL heights determined with the bulk Richardson number method.In the second case, the "Hanna" approach was used again, but with PBL heights determined from TKE profiles.In the third case, the "Direct" approach was used, which does not use the PBL height explicitly as input variable for the turbulence calculations.The hypothetical release of Cs-137 was made in northern Switzerland on 8 September 2008 between 00:00 UTC and 06:00 UTC, with an emission rate of 46290 MBq/s, and the pollutant transport was calculated for 18 h.Figure 5 shows the forecasted nearsurface (below 500 m a.g.l.) mean concentration fields for the three simulations on 8 September 2008 18:00 UTC.Highest concentrations (maximum: 258 Bq/m 3 ) occur in the case of the first simulation ("Hanna" with Ri bulk).If PBL heights are derived from TKE profiles, the resulting cloud is more dispersed, and smaller concentrations (max: 163 Bq/m 3 ) are predicted.The lowest concentrations (max: 81 Bq/m 3 ) are simulated with the "Direct" approach.Comparing only the first two simulations, it is observed that changing the diagnostic approach for the PBL height leads to a difference in maximum concentrations by a factor of about 1.6.This is due to the fact, that PBL heights diagnosed from TKE profiles are higher than those calculated with the bulk approach with PBL heights determined with the bulk Richardson number method; (b) "Hanna" approach with PBL heights determined from TKE profiles; (c) "Direct" approach.Richardson number method (Fig. 6).Since the "Direct" approach implicitly contains a PBL height definition based on TKE, the third simulation should be compared to the second one.The difference between these two runs is due to the different interfacing approach.For the case investigated, the "Direct" approach results in lower maximum concentrations by a factor of 2. This can mainly be traced back to the stronger near-surface turbulence diagnosed in this approach.

Conclusions
The impact of two different interfacing approaches on dispersion has been studied with the COSMO-LPDM system operationally used at MeteoSwiss.One of the interfacing approaches estimates velocity variances directly from the TKE profile of the COSMO model ("Direct"), while the other employs similarity relations ("Hanna") and thus requires the boundary layer height as an input.An extensive validation exercise for methods to diagnose the PBL height had rendered methods based on the bulk Richardson number or momentum flux profiles of the NWP model as the two best approaches.The diagnosed turbulence variables are evaluated on three turbulence measurement campaigns, with the "Hanna" approach performing slightly better.Simulations of hypothetical pollutant releases showed that the different interfacing approaches can lead to substantial changes in the forecasted concentrations.It is planned to evaluate the COSMO-LPDM system on real dispersion experiments to further investigate the relative performance of the two interfacing methods.

Figure 1 .
Figure 1.Relative root mean square error (RMSE) of the diagnosed PBL height from the COSMO-7 (blue) and COSMO-2 (red) model for a one-month period in March 2008.Model results are validated with 5 radio sounding stations in Europe.Left: unstable situations, right: stable situations.Methods: gradient Ri number method ("Ri"), bulk Ri number method ("Bulk Ri"), TKE profile (TKE rel), Momentum flux profile ("Mom.Flux"), heat flux profile ("Heat flux"), Slab model, Zlitinkevich method ("Zil.").Applied thresholds for the different methods are indicated in parentheses.

Figure 2 .
Figure 2. Profiles of the (a) horizontal and (b) vertical standard deviations of wind fluctuations for the ideal convective case.Turbulence variables are diagnosed with two different approaches ("Direct" and "Hanna") from single column simulations of the COSMO model.The reference is output from LES simulations according to Mironov et al. (2000).

Figure 3 .
Figure 3.Time series of the standard deviations of cross-wind (left) and vertical (right) wind fluctuations at 90 m height for 30 May 2003 in Falkenberg, Germany (time in UTC).Red and yellow lines: COSMO single column simulations ("Hanna" and "Direct" approach).Black line: tower turbulence measurements.

Figure 4 .
Figure 4. Relative bias (upper panel) and relative standard deviation (lower panel) scores for the four measurement sites (in different colours).Verified parameters: wind speed, standard deviation of horizontal and vertical wind fluctuations ("Direct" approach: sigma x, sigma y and sigma z; "Hanna" approach: sigma u, sigma v, sigma w), Turbulent Kinetic Energy (TKE).

Figure 5 .
Figure5.Forecasted near-surface (below 500 m a.g.l.) mean concentration fields by the COSMO-7 -LPDM system for 8 September 2008 18:00 UTC (18 h forecasts, hypothetical case).Simulations differ from each other in the interfacing approach applied: (a) "Hanna" approach with PBL heights determined with the bulk Richardson number method; (b) "Hanna" approach with PBL heights determined from TKE profiles; (c) "Direct" approach.

Figure 6 .
Figure 6.Mixing heights (m AGL) diagnosed from outputs of the COSMO-7 model for three times indicated on 8 September 2008.Left: bulk Richardson method; right: mixing heights diagnosed from TKE profiles.

Figure 6 .
Figure 6.Mixing heights (m a.g.l.) diagnosed from outputs of the COSMO-7 model for three times indicated on 8 September 2008.Left: bulk Richardson method; right: mixing heights diagnosed from TKE profiles.