Doppler radar radial winds in HIRLAM. Part I: observation modelling and validation

Abstract An observation operator for Doppler radar radial wind measurements is developed further in this article, based on the earlier work and considerations of the measurement characteristic. The elementary observation operator treats radar observations as point measurements at pre-processed observation heights. Here, modelling of the radar pulse volume broadening in vertical and the radar pulse path bending due to refraction is included to improve the realism of the observation modelling. The operator is implemented into the High Resolution Limited Area Model (HIRLAM) limited area numerical weather prediction (NWP) system. A data set of circa 133 000 radial wind measurements is passively monitored against theHIRLAM six-hourly background values in a 1-month experiment.No data assimilation experiments are performed at this stage. A new finding is that the improved modelling reduces the mean observation minus background (OmB) vector wind difference at ranges below 55 km, and the standard deviation of the radial wind OmB difference at ranges over 25 km. In conclusion, a more accurate and still computationally feasible observation operator is developed. The companion paper (Part II) considers optimal super-observation processing of Doppler radar radial winds for HIRLAM, with general applicability in NWP.


Introduction
Weather radar is a remote sensing instrument which emits electromagnetic pulses and receives backscattered radiation from atmospheric targets. Properties of the backscattered radiation, such as intensity, phase shift and polarization, can be interpreted as atmospheric observations. Radar measurements do not directly represent atmospheric state variables, like most in situ measurements of the atmosphere, and further processing and analysis of the measurements is therefore necessary. This applies generally to all remote sensing data.
Meteorological interpretation of the radar measurements is usually obtained through geophysical inversion, where the unknown atmospheric variables are resolved from the known measurements. A well-known example of such an inversion in radar meteorology is the Z-R relationship, where reflectivity (Z) is related to rain rate at the surface (R) through regression (Marshal and Palmer, 1948). The velocity azimuth display (VAD; Lhermitte and Atlas, 1961;Browning and Wexler, 1968) and the volume velocity processing (VVP; Waldteufel and Corbin, 1979) techniques are examples of geophysical inversion methods for Doppler radar radial wind observations. These techniques interpret the Doppler radar radial wind velocity at different antenna azimuth directions as a horizontal vector wind, under certain assumptions on the true wind field in the vicinity of the radar.
The difficulty with the geophysical inversion of the radar measurements is that atmospheric phenomena which the radar measures are very complex. Even though the radar provides measurements in abundance, they may still not be sufficient alone for resolving the inversion problem. In other words, there are more degrees of freedom in the atmosphere than in the measurements, and thus the set of possible solutions can be infinite. One way to constrain the number of possible solutions is to introduce additional data to the geophysical inversion process. The inversion can be supplemented, for example, with other observations, or with atmospheric model data, to provide a complementary view to the atmospheric parameters. A systematic method to do this is data assimilation.
Data assimilation consists of modelling the measurement with the numerical weather prediction (NWP) model variables, constraining the solution by the dynamics of the atmosphere, and by using all other available observations and a priori information of the atmosphere simultaneously in the inversion process. These information sources put together are sufficient for obtaining a unique solution for the atmospheric variables with less severe assumptions on the properties of the solution. Data assimilation enables exploitation of a variety of indirect observations of the NWP model variables, such as satellite infrared radiances (e.g. Eyre, 1990), GPS zenith total delays and slant delays (e.g. Vedel and Huang, 2004;Eresmaa and Järvinen, 2006), radar reflectivities (e.g. Caumont et al., 2006) and Doppler radar radial winds (e.g. Sun and Crook, 1997;Salonen et al., 2003;Lindskog et al., 2004;Seko et al., 2004;Rihan et al., 2005). A common factor to the exploitation of these data types is observation modelling, where the observed quantities are expressed in terms of the model variables.
Observation modelling requires detailed understanding of the measurement. In case of Doppler radar radial wind measurements, the observation represents a weighted average radial wind speed in the radar pulse volume. There are several aspects which need to be considered in the modelling, such as measurement geometry, bending of the radar pulse path due to atmospheric refraction, broadening of the pulse volume, distribution of the atmospheric hydrometeors in the pulse volume and fall speed of the scattering hydrometeors. Also, as compared with many other atmospheric observation types, radar measurements are available in abundance. This enables either to use thinned raw measurements in NWP or to pre-process the raw measurements into spatial mean observations, hereafter called super-observations. The idea in super-observation processing is to average out random errors which appear as erratic noise from measurement pixel to pixel.
This article describes the modelling of the Doppler radar radial wind measurement in High Resolution Limited Area Model (HIRLAM; Undén et al., 2002) framework and introduces refinements to the basic version of the observation operator utilized in Lindskog et al. (2000) and Lindskog et al. (2004). The aim of this paper is to improve the modelling of the radial wind observations and to statistically intercompare different observation operator versions, with an objective to reduce the systematic and random observation minus background (OmB) model difference. The objective of the companion paper (Part II) is to reduce the random OmB difference by means of optimized super-observation processing of the raw measurements. Validation indicates that these objectives are achieved.
This article is organized as follows. Section 2 describes the Doppler radar radial wind measurement and Section 3 its error sources from the modelling point of view. Section 4 concentrates on the observation operator development. Test results are given in Section 5, followed by summarizing discussions in Sections 6 and 7.

Doppler radar radial wind measurement
A Doppler radar emits compact electromagnetic pulses which backscatter from atmospheric hydrometeors. The hydrometeors are in three-dimensional motion and their radial velocity with respect to the radar can be determined from the phase shift between the backscattered returns of successive radar pulses (Doviak and Zrnić, 1993).
The received backscattered signal is a composite of echoes from a very large number of hydrometeors in a scattering volume. The scattering volume has the same cross-section as the radar beam at certain range and a depth roughly equal to half the length of the emitted pulse. Each hydrometeor in the scattering volume has its own speed and direction determined by its size, shape and motion of the surrounding air. The velocity distribution of the received signal is called Doppler spectrum. The width of the spectrum is a function of radar system parameters, such as pulse width and radar frequency, and meteorological parameters that describe the distribution of the hydrometeor density and velocity within the scattering volume. For instance, turbulence, wind shear and different hydrometeor fall speeds cause relative radial motion of the hydrometeors which tends to broaden the spectrum. In short, the radial wind velocity is neither an instantaneous value, nor it is a point measurement.
A major uncertainty related to the Doppler radar radial wind measurements is their velocity-range ambiguity. If the radial velocity is such that the phase difference between the successive pulses is greater than π radians, there is an ambiguity concerning the measurement of the phase and, consequently, of the radial wind. This phenomenon is characteristic to measurements with radar systems applying a constant pulse repetition frequency (PRF). Measurement techniques, such as, dual-PRF (Dazhang et al., 1984), staggered pulse repetition time (Sirmans et al., 1976) and simultaneous multiple PRF code (Pirttilä et al., 2005) can be used to alleviate the effects and limitations caused by the ambiguity problem. Measurements can be corrected also afterwards by using dealiazing algorithms (e.g. Ray and Ziegler, 1977;Haase and Landelius, 2004). In the following, it is assumed that the radial wind measurements are unambiguous.

Observation errors
Atmospheric observations are typically irregular in their spatial and temporal distribution, whereas model variables are regularly discretized, for instance, in a model grid at coherent time levels. Furthermore, observations and model variables usually represent different physical quantities. For example, the Doppler radar measures the radial wind component only, whereas the model variable is vector wind.
In data assimilation, comparison of observations with model data is arranged in the following way. Observations y are modelled with an observation operator H from model state variables x within errors such that The observation operator models the observations with the right dimension through interpolation of the model variables to the observation locations, and with the right physical quantity through transformation of the model variables to the observed quantity. Assuming that x represents the atmospheric true state (i.e. omitting analysis and forecast errors), the observation operator is accurate within observation errors o , which consist of instrumental, modelling and representativeness errors. In case of a Doppler radar, systematic instrumental errors are typically due to inaccurate radar calibration or antenna pointing. These errors cannot be modelled, and they should be minimized prior to the use of the measurement data in NWP. A fraction of these biases may remain in the data and they can be estimated and partly removed in a specific bias removal in operational NWP, as data assimilation algorithms assume random gaussian observation errors with zero expectation (e.g. Kalnay, 2003). Random instrumental errors are due to fluctuations in the functioning and data processing of the instrument, and their impact can be reduced through optimizing the super-observation processing (cf. objectives of Part II).
Modelling errors are caused by deficiencies in the observation operator design. The operator interpolates the model variables to the observation locations and transforms them to the observed quantities. The modelling error thus consists of interpolation and numerical errors, and of errors in the physical modelling of the measurement. These errors contribute both to the systematic and random observation errors. The objective here is to reduce the systematic modelling error by careful observation operator design. Computation of the model radial wind from more than one model gridpoint value tends to average out random modelling errors which appear as erratic noise from one model gridpoint to another. This is in analogy to the smoothing of random errors in the super-observation processing of raw measurements. Thus, the observation operator development can also reduce the random part of the modelling error.
Representativeness errors are due to atmospheric phenomena which are detected by the instrument but are not represented in the finite vector of the discretized model variables. Representativeness errors are considered here to be mostly random, and their impact cannot be reduced through observation modelling. The viewpoint is that the impact of random representativeness errors may be partly alleviated by optimizing the super-observation processing (cf. objectives of Part II).
In summary, two strategies to minimize the observation errors are explored here. First, careful design of the observation operator is applied to minimize the systematic and random modelling errors. Second, super-observation processing (Part II) is applied and optimized to reduce the random part of the instrumental and representativeness errors.

Observation operator
This section describes in detail the observation operator developed for the Doppler radar radial wind measurements. In the following, it is assumed that the motion of air and the fall ve-locity of hydrometeors are inseparable. The separation would require knowing the microphysical composition of the hydrometeor population in the finest details. The inseparability assumption implies that it is not possible to determine the vertical wind component from the vertical velocity of the hydrometeors alone, and that the horizontal motion of the hydrometeors represents the horizontal motion of the air. This assumption restricts the following treatment to Doppler radar radial winds measured at low radar antenna elevation angles, low meaning elevation angles where the projection of the vertical wind component on the radar pulse path can be assumed negligible.
A pre-processing step prior to the observation modelling of radar observations expresses the observation location in terms of geographical position and height. A raw radar measurement contains information on the measurement range, and antenna azimuth and elevation angles. The observation location and height can thus be obtained with trigonometrical considerations and applying a 4 3 r-law for the radar pulse path, where r is the radius of the Earth (Doviak and Zrnić, 1993). The 4 3 r-law assumes standard atmospheric refraction conditions of the so-called International Civil Aviation Organization (ICAO) atmosphere. Actual bending of the radar pulse path depends on the atmospheric refractivity profile along the pulse path, and therefore some uncertainty remains in the pre-processed measurement height at a given range. This aspect is considered more in Section 4.3.

Radial wind component
The observation modelling task is to express the observed Doppler radar radial wind component in terms of appropriate model variables. The starting point here is the algorithm of Lindskog et al. (2000) and Lindskog et al. (2004), which considers observations as point measurements.
The algorithm is such that, first, NWP model profiles of the horizontal wind components u and v are interpolated to the observation location with a bi-linear interpolation. A further linear interpolation is performed in vertical (linear in ln p, i.e. linear in mass) to the observation height. The vertical interpolation is essentially a convolution between a full model wind profile and a 'hat' function which has exactly two model levels width and which peaks at the observation height. Second, the interpolated wind components are projected on a horizontal plane towards the radar using where θ is the radar antenna azimuth angle. Finally, v h is projected on a vertical plane to the slanted direction towards the radar using where v r is the model counterpart of the measured Doppler radar radial wind component and φ is the radar antenna elevation Tellus 61A (2009) angle. The argument α approximately accounts for the curvature of the Earth through the formula where h is the height of the radar antenna above the mean sea level, d is the measurement range and r is the radius of the Earth (Doviak and Zrnić, 1993). Observing geometry of this projection is illustrated in Fig. 1 [Note, that in Lindskog et al. (2000) and Lindskog et al. (2004) the term 4 3 r in the denominator of (4) was erroneously r. This causes an error of the order of 10 −5 ms −1 to the modelled radial wind component].
This modelling of the radial wind observation implicitly assumes that a radar measurement essentially represents a point in space, not a volume. As was discussed in Section 2, the radar pulse volume is spacious and further broadening takes place as the pulse propagates away from the radar antenna. This modelling aspect is considered next.

Radar pulse volume and broadening
The shape of the emitted radar pulse in the azimuth and elevation directions is determined by the radar antenna properties. A typical half-power width of a pulse, often called the radar pulse main lobe, is ∼1 • . The power distribution in the radar pulse main lobe is approximately gaussian (Probert-Jones, 1962). The directional gain of the antenna has to be taken into account for transmission and reception when calculating the shape of the scattering volume. As the same antenna is used for both, the shape of the scattering volume is obtained by squaring the antenna pattern. Thus, if the pulse is of the width of ∼1 • , the half-power width of the scattering volume is ∼0.7 • . This is often referred to as a two-way beamwidth. As the angular width of the beam remains constant, the scattering volume broadens in the transverse direction as the pulse propagates away from the antenna.
The power distribution in the radial direction is determined by the pulse shape, the receiver impulse response and the processing performed in the radar signal processor. A full theory of radar ambiguity functions is presented in Lehtinen and Huuskonen (1996). Here, it suffices to note that the radial weighting function is a square of the convolution of the pulse form and the receiver impulse response, modified by any averaging performed. The scattering volume shape remains unchanged in the radial direction when the pulse propagates away from the antenna. Figure 2 shows the weighting function calculated for the measurements used in this study. The pulse form and the impulse response are based on theoretical functions calculated with the real pulse length of 0.6 μs and the receiver bandwidth of 800 kHz. The resulting range weighting function has a halfpower width of 110 m, and an extent of 360 m, which corresponds to the sum of the lengths of the pulse and the impulse response. In the final averaging stage, 12 consecutive estimates, separated by 83 m, are averaged. The resulting weighting function has a half-power width of 1000 m, and a total extent slightly below 1300 m. The size of the scattering volume varies greatly depending on a selection of available pulse length and processing parameters. It should be noted that the length of 1000 m is valid for this case only. In practice, widths from 200 m upwards can be encountered in different radar installations.   Next, an approach to model the vertical broadening of the radar pulse volume is described. In this approach, the 'hat' function is replaced with a gaussian averaging kernel in the vertical interpolation. The gaussian averaging kernel reads where z is the model level height and z 0 is the observation height.
The argument κ defines the kernel width. Here z k is the height of the upper limit of the pulse volume half-power width, calculated with the 4 3 r-law at measurement range d. Figure 3 illustrates the pulse volume broadening and the vertical averaging kernel at ranges of 50 and 150 km. At the range of 50 km (150 km), the half-power width of a purely gaussian volume is ∼600 m (∼1800 m). Clearly, the pulse volume extends over several model levels (denoted by dotted lines in Fig. 3) in the vertical at these ranges from the antenna.
Horizontal broadening of the pulse volume also takes place. At a range of 50 km (150 km), the half-power width of the gaussian volume is ∼600 m (∼1800 m) whereas the length of the pulse volume is independent of range (1000 m in our case). Depending on the grid spacing of an NWP model, this horizontal extent of the pulse volume needs or needs not to be accounted for. As long as the pulse volume is small compared to an NWP model grid cell, a bi-linear horizontal interpolation to the radar pulse volume centre point is expected to be accurate enough.
An unwanted feature of the gaussian averaging kernel is that it is non-zero from the Earth's surface up to the top of the atmo- sphere. Obviously, only the wind information which the radar is able to measure should be included into the model counterpart.
Lower and upper bounds are therefore introduced to the vertical averaging kernel which is set to zero outside these bounds. First, the radar is unable to see below radar horizon. This obscuring effect is taken into account by assuming a radar horizon of 0 • elevation angle, below which the model information is not used. The lower dashed line in Fig. 3 denotes the lower bound of the averaging kernel.
Second, the observation operator with a gaussian averaging kernel implicitly assumes a homogeneous field of scatterers between the radar horizon and the top of the model grid. To avoid this unrealistic modelling, an empirical upper bound is set to the height of 1.5 times the pulse volume half-power width (Jarmo Koistinen, personal communication, 2003). This is based on an empirical fact that the upper part of the radar pulse volume is located, at longer ranges and with higher elevation angles, above the scattering hydrometeors or in a region where the radar reflectivity factor is several orders of magnitude smaller than in the lower parts of the pulse volume. Hence, the upper part of the pulse volume will not significantly contribute to the measured volume-weighted reflectivity or the Doppler wind. The upper limit for the averaging kernel is denoted by a dashed line in Fig. 3.

Bending of the radar pulse path
As discussed earlier, the radar pulse path bends due to atmospheric refraction. This has two effects on the observation modelling. First, the observation height at range d is different from the one obtained by applying the 4 3 r-law, and second, the elevation angle of the radar pulse path gradually differs from the radar antenna elevation angle. The actual elevation angle of the radar pulse path is hereafter referred as the effective elevation angle.
The radar pulse path bending can be calculated using the Snell's law where n 1 (n 2 ) is the refraction index in the first (second) layer and φ 1 (φ 2 ) is the radar pulse path elevation angle in the first (second) layer. Local refraction index n can be calculated from the NWP model temperature T, pressure p and water vapour partial pressure e for microwave frequencies using n = 1 + 77.6 × 10 −6 p T + 0.373e T 2 (8) (Doviak and Zrnić, 1993). In this observation operator implementation, the total bending of the radar pulse path across the model levels is accumulated until the pulse reaches the observation location. The last calculated elevation angle is used in (3) as the effective elevation angle when projecting the horizontal wind v h on the slanted Tellus 61A (2009), 2 direction towards radar. Also, the pre-processed observation height obtained with the 4 3 r-law at range d is replaced with the observation height obtained through the refraction calculation. For simplicity, it is assumed that the refraction index profile calculated for all NWP model full levels at the radar location is representative also at the measurement location.

Numerical tests
This section presents results from a 1-month validation experiment for January 2002. In the experiment, Doppler radar radial wind measurements are passively monitored against the modelled radial wind values using the HIRLAM data assimilation and forecasting system. Three different observation operator versions are used: (i) the basic version given by (3) and (4) (hereafter referred as the BASIC-version); (ii) the basic version plus modelling of the radar pulse volume broadening (VBROAD-version) and (iii) a version further improved by modelling of the pulse path bending (VBEND-version). The passive monitoring allows to evaluate the performance and accuracy of these three different observation operator versions. No data assimilation experiments are performed in this validation test.

The data
The measurement data set consists of circa 133 000 Doppler radar radial wind observations and their HIRLAM model counterparts. Quality of the radar wind observations varies from radar to radar and depends among other things on the radar siting. To minimize the effects of radar dependent error sources, and to make a reliable comparison of the observation operator versions, observations from one radar at Karlskrona (56.30 • N, 15.61 • E), Sweden, are considered here.
The measurements are made with the dual-PRF technique at four to three ratio (1200 s −1 /900 s −1 ). The resulting maximum measurement range is 120 km, and the unambiguous velocity interval is ±48 ms −1 . Elevation angles of 0.5 • , 1.1 • , 2.3 • and 3.2 • are included in the data set. Resolution of the raw data is 1 km in range and 0.9 • in azimuth (420 azimuth gates per 360 • scan). For this study, raw observations are spatially averaged so that the nominal resolution of the super-observations is 10 km in range and 1.7 • in azimuth. The spatial averaging has a tendency to smoothen the observed wind field and average out random errors. Details of the averaging procedure can be found in Lindskog et al. (2004) and a justification for this choice of averaging parameters from the Part II of this article.
The model counterparts for the radial wind observations are calculated from the 6-h forecasts obtained with the hydrostatic HIRLAM model version 7.1alpha1 (Yang, 2007). The model is run with 11 km horizontal grid length and with 40 vertical levels. The integration area covers Northern Europe including the northern part of Central Europe. The HIRLAM three-dimensional variational data assimilation system Lindskog et al., 2001) is applied. In the experiment, only conventional observations of temperature, specific humidity, wind and geopotential at the pressure level of the observation site are assimilated with six-hourly cycling.

Vertical pulse volume broadening
First an illustration of the effects of applying linear interpolation and gaussian averaging kernel in the vertical interpolation is given in a selected case. In the left-hand panel of Fig. 4 In the right-hand panel of Fig. 4, the observation is made at a range of 100 km from the radar at a pre-processed observation height of 1460 m. There is a local minimum in the model wind profile near the observation height. Now, the model v wind component obtained with the BASIC-version is −2.3 ms −1 , and with the VBROAD-version it is −3.2 ms −1 . There is thus a considerable difference of 0.9 ms −1 in the interpolated wind components due to the different interpolation methods. In conclusion, the impact of the gaussian averaging kernel in the vertical interpolation is large in case there is a local minimum or maximum in the model wind profile near the observation height. The gaussian averaging kernel accounts for model winds at several model levels and is in that sense faithful to the radar measurement, under the assumption of a homogeneous field of scatterers.
Next, the whole 1-month data set is considered. The statistics calculated over the data set reveal that the mean difference between the wind components obtained with the BASIC versus the VBROAD-version is very close to zero. In extreme cases, covering less than 0.8% of the observations, the difference is over 2 ms −1 . These cases originate mainly from long measurement ranges where the gaussian averaging kernel is quite wide.

Radar pulse path bending
Modelling the radar pulse path bending affects both the observation height and the effective elevation angle in (3) and (4). Impacts of these effects are considered next.
First, the impact of the radar pulse path bending on the observation height is considered. Differences between the observation heights obtained with the 4 3 r-law versus the Snell's law are computed for the 1-month data set. The mean difference of the observation heights increases with measurement range (Fig. 5), as does the standard deviation around the mean value. The mean and the standard deviations are less than 20 and 36 m, respectively, for all measurement ranges. On average, the effect is thus small. However, minimum and maximum differences reach values as high as −105 and 215 m, respectively. In such cases, the peaking height of the averaging kernel is shifted significantly from the pre-processed height. Depending on the shape of the model wind profile, the impact on the modelled radial wind can occasionally be large.
Second, the impact of the radar pulse path bending on the effective elevation angle is estimated. Let us consider normal atmospheric refraction conditions of the ICAO atmosphere for a hypothetical measurement at a radar antenna elevation angle of 1.5 • and a measurement range of 66 km, implying a measurement height of 2 km. The effective elevation angle derived from (7) is 1.25 • . As a result, the term cos(φ + α) in (3) changes by a factor of 10 −4 . In extreme super-refraction conditions for low elevation angle measurements, changes in the term cos(φ + α) can be of the order of 10 −2 . More generally, in the 1-month data set the mean difference between the radar antenna elevation angle and the effective elevation angle is very close to zero with a standard deviation of 0.07 • . In 0.5% of all cases the difference is over 0.3 • with a maximum difference of 0.48 • , which changes the term cos(φ + α) by a factor of 10 −4 . Magnitudes of these changes are very small and hardly impact the modelling results.
In conclusion, modelling of the radar pulse path bending mainly impacts through changes in the observation height and to a minor extent through changes in the effective elevation angle.

Intercomparison of the observation operator versions
Accuracy of the observation modelling is evaluated next. Observation minus background (OmB) statistics are studied for the three observation operator versions in the 1-month data set.
Obtaining an estimate of the mean OmB (bias) for the radial wind components is a non-trivial issue. The conventional way of calculating the radial wind bias by summing up OmB values at different azimuth directions results in a near zero bias even when there are systematic differences in the modelled and observed wind speed and/or direction. A bias estimation method which accounts for this aspect (Salonen et al., 2007) is thus applied. This method provides a bias estimate for vector wind (i.e. wind speed and direction bias) instead of radial wind. Calculation of the OmB standard deviation for radial winds is straightforward as there is no such azimuthal dependence of the random errors. In the following, observations up to a measurement range of 105 km are considered since the amount of the data at longer measurement ranges is too small for a reliable bias estimation.
The 95% confidence intervals for the OmB statistics are calculated by using a bootstrap method (Efron, 1982;Efron and Gong, 1983). The bootstrap method is chosen as it is the most suitable method for determining the confidence intervals in the case of radial wind observations. The vector wind bias and the standard deviation of the radial wind component represent different quantities, thus the use of conventional t-test or χ 2 -test is not convenient. The sample size used in bootstrapping is 10 000. This value is considered to be high enough to obtain reliable estimates for the confidence intervals. The wind speed bias (not shown) is negative for all observation operator versions at all measurement ranges, such that the observed wind speed is on average weaker than the modelled one. The bias is about −0.12 ms −1 at a range of 15 km, and reaches −0.4 ms −1 (−0.9 ms −1 ) at a range of about 45 km (85 km). At measurement ranges of 15-85 km, the VBROADand the VBEND-versions perform equally well, while for the BASIC-version the wind speed bias is slightly (by up to about 0.06 ms −1 ) larger.
The wind direction bias (not shown) varies between −1.4 • and 1.7 • . The bias is positive at measurement ranges shorter than about 55 km and negative at longer ranges. VBROADand VBEND-versions perform equally well at all measurement ranges. Up to a range of about 50 km, the wind direction bias of the BASIC-version is the largest (by up to about 0.3 • ). At ranges 55, 65 and 95 km, the bias of the BASIC-version is the smallest (by up to about 0.2 • ).
The vector wind bias (Fig. 6) generally increases with the measurement range from 0.5 ms −1 , and below, to about 1.1 ms −1 . It appears that, apart from the ranges of 65 and 95 km, the smallest vector wind bias is obtained with the VBROADversion. The bias of the VBEND-version is equal or marginally larger than the bias of the VBROAD-version. The 95% confidence intervals of the bias estimates are denoted by vertical bars (for the BASIC, VBROAD and VBEND on the left-hand side, middle and right-hand side bar, respectively). At longer ranges, the confidence intervals become wider due to smaller data amounts. At a range of 45 km, for instance, the bias estimate of the BASIC-version is on the edge of the confidence intervals for the VBROAD-and VBEND-versions. This can be interpreted such that at 95% confidence level, vector wind bias of the VBROAD-and the VBEND-versions is smaller than for the BASIC-version. Generally, most bias estimates are within the 95% confidence intervals and statistically significant conclusions about the differences between the observation operator versions cannot be drawn.
The OmB standard deviation of the radial wind speed (Fig. 7) generally increases with the measurement range from about 3.4 to about 4.9 ms −1 . Apart from the range of 15 km, the standard deviation is the smallest for the VBROAD-version. Virtually equally small standard deviation is obtained for the VBENDversion. However, the 95% confidence interval bars overlap, and the differences between the observation operator versions are not statistically significant. At a range of 95 km, the 95% confidence threshold is almost reached.

Discussion and conclusions
Three different observation operator versions for Doppler radar radial wind measurements have been validated based on a 1-month data set of OmB differences. Two OmB statistics are used: the mean OmB difference for vector wind (bias) and the standard deviation of the OmB difference for the radial wind component (random error).
The validation reveals that the elementary observation operator (point measurement, pre-processed observation height) is a suitable option due to its relatively good validation result and conceptual simplicity. The validation indicates also that it is beneficial to improve the elementary observation operator by modelling the Doppler radar pulse volume broadening. This modelling results in reduced vector wind bias at ranges below 55 km, and reduced random error at ranges above 25 km. Finally, the validation reveals that an additional modelling of the Doppler radar pulse path bending has a neutral, or very slightly detrimental, effect on the validation statistics.
In the 1-month validation data set, bias and random error tend to increase with the measurement range for all versions of the observation operator, and the modelled radial wind tends to be stronger than the measured one. This might occur -bearing in mind that winds generally increase with height in the troposphere -if the observed wind in fact originates from lower levels than is assumed in the observation operator. One possible mechanism is as follows: the radar pulse volume broadens with increasing range and the pulse volume may be partly, or even mainly, above the scatterers. Clearly, the probability for this overshooting increases with the measurement range. In such cases, the measurement originates from the lower part of the pulse volume and the actual measurement height is lower than the modelled one [assuming either the two model levels closest to the 4 3 r-law observation height (the BASIC-version) or assuming a homogeneous field of scatterers between the lower and upper limits defined for the gaussian averaging kernel in the VBROAD-and VBEND-versions].
Modelling of the radar pulse path bending has a neutral, or even slightly detrimental, impact on the validation statistics. The neutral impact is easy to understand since in the 1-month data set the mean difference between the observation heights obtained with the 4 3 r-law and the Snell's law is less than 20 m for all measurement ranges. Thus, on average, the modelled winds are computed from very similar heights regardless of this modelling detail. In individual super-refraction cases, the effect of accounting for the pulse path bending can be much higher. Furthermore, the model winds originate from a fairly coarse resolution NWP model (11 km horizontal resolution, 40 vertical levels). As model resolutions are increasing, one can expect the refractive profiles to be better represented in the future models, and justify modelling of the radar pulse path bending.
Our hypothesis is that, first, physically more faithful observation modelling reduces the systematic modelling error and that, second, use of more than one gridpoint value in the computation of the model radial wind with the observation operator reduces the random modelling error. The statistical comparison of the different observation operator versions indicates that this is indeed the case.
The reduction of the random modelling error occurs as follows. The gaussian averaging kernel smoothens the random errors contained in the model u and vprofiles more effectively than the linear interpolation which only uses the two nearest model levels. This mechanism has been demonstrated by generating 10 000 synthetic wind profiles by adding random errors with zero expectation and standard deviation of 2 ms −1 on a 'true' wind profile. Modelled radial wind values calculated from the synthetic and 'true' wind profiles, respectively, have been compared. The mean difference is zero both for linear interpolation and for gaussian averaging kernel, as it should be for random errors with zero expectation. However, the standard deviation decreases significantly when the gaussian averaging kernel is applied.
At the European Centre for Medium-Range Weather Forecasts (ECMWF) IBM p690+ computer, the consumed CPU times per observation are for the BASIC-, VBROAD-and VBENDversions 0.1, 0.3 and 0.4 ms, respectively. The VBEND-version requires thus four times more CPU time than the BASIC-version, but in general, the CPU times are small for all versions.
In summary, the recommendation is to use the elementary observation operator if a quick implementation (and quick testing of correct functionality) is emphasized. For obtaining the best possible validation results and physical realism, especially in fine resolution modelling, the observation operator should also include the Doppler radar pulse volume broadening and the pulse path bending.

Summary
Modelling of the Doppler radar radial wind measurements has been discussed in detail in this article. The observation modelling effort builds on earlier work (Lindskog et al., 2000(Lindskog et al., , 2004 which represents the radar pulse volume by a point, and projects the corresponding model horizontal wind on the slanted path towards the radar antenna.
Two refinements of the observation operator are introduced here: first, modelling of the vertical broadening of the radar pulse volume, and, second, modelling of the pulse path bending due to atmospheric refraction. The vertical pulse volume broadening affects the modelled radial wind the most in cases where there is a local minimum/maximum in the NWP model wind profile. The pulse path bending affects the value of the modelled radial wind mainly through modifying the height assigned to the observation.
The observation operator is implemented into the HIRLAM limited area NWP system. Doppler radar radial wind measurements have been passively monitored against the modelled radial wind values in a 1-month validation experiment (January 2002). Statistics of the OmB values improve such that the modelling of the pulse volume broadening reduces the mean vector wind OmB difference at ranges below 55 km and the standard deviation of the radial wind speed OmB difference at ranges above 25 km. Modelling of the pulse path bending has a neutral, or very slightly detrimental, effect on the validation statistics.
In conclusion, an accurate and computationally feasible observation operator has been developed in Part I of this article. Part II considers optimizing the pre-processing of Doppler radar radial wind observations prior to exploiting them in data assimilation.
The observation operator computer codes are written in FORTRAN90 programming language and are available on request from the corresponding author.