Assimilating ambiguous QuikScat scatterometer observations in HIRLAM 3-D-Var at the Norwegian Meteorological Institute

The Norwegian Meteorological Institute has implemented ambiguous QuikScat observations in the HIRLAM threedimensional variational data assimilation system, by approximating the method that will achieve the best mean squared error verification results according to the Bayesian decision theory. A four-month long impact study of this implementation during the 2003/2004 winter showed that assimilating QuikScat observations had a small positive effect on the average mean sea level pressure forecasting skill.


Introduction
Placed on board the QuikScat polar orbiting satellite, the Sea-Winds scatterometer instrument covers large ocean areas giving valuable wind observations for operational meteorology (Atlas et al., 2001;Liu, 2002). The scatterometer is a radar, emitting microwave pulses towards the sea surface and measuring the backscattered signal from different directions as the satellite passes by the target area. The backscatter properties of the sea surface depend on the small capillary sea waves. These waves in turn depend on the near surface wind. An empirical relationship can be found between the backscatter measurements and the 10-m wind speed and wind direction relative to the antenna-viewing angle. It turns out that several wind vector solutions are typically in agreement with a set of backscatter measurements. A scatterometer observation is therefore usually presented as a set of wind vector solutions, called ambiguities. The NOAA/NESDIS (National Oceanic and Atmospheric Administration/National Environmental Satellite, Data and Information Service) provides QuikScat observations operationally on BUFR (Binary Universal Form for the Representation of meteorological data) format to the user community as backscatter measurements and the corresponding ambiguous wind vectors on 25-km resolution. In the context of the EUMETSAT/NWP SAF (European Organisation for the Exploitation of Meteorological Satellites/Numerical Weather Prediction Satellite Application Facility), the KNMI (the Koninklijk Nederlands Meteo- Corresponding author. e-mail: f.t.tveter@met.no rologisch Instituut) provides software for averaging the 25-km resolution QuikScat backscatter measurements to 100 km, and retrieving the corresponding wind vector ambiguities (de Vries, 2004). By averaging the backscatter measurements, the retrieved winds become less ambiguous, and therefore more suitable for assimilation purposes.
Different techniques have been developed for assimilating scatterometer data in numerical weather prediction systems. Stoffelen and Anderson (Stoffelen and Anderson, 1997) developed a preprocessor that selected the ambiguities that gave the most spatially consistent wind field. These ambiguities were in turn assimilated in the ECMWF (European Center for Medium-Range Weather Forecasting) model as observations with Gaussian error characteristics. Portabella and Stoffelen (Portabella and Stoffelen, 2004) assimilated QuikScat observations by introducing a non-quadratic cost function to a (two-dimensional) variational analysis, similar to a method used, for instance, by Isaksen and Janssen (Isaksen and Janssen, 2004). They based their cost function on the principle that the resulting 'analysis' should be the 'most probable state of the atmosphere', as described by Lorenc (Lorenc, 1988). This principle had been suggested earlier by Lorenc (Lorenc, 1986) as one of two possible data assimilation principles, where the second principle was to use the 'expected state of the atmosphere' as the 'analysis'. Ingleby and Lorenc (Ingleby and Lorenc, 1993) used Bayesian decision theory in a discussion of the two assimilation principles, identifying the two different 'loss functions' for which they are optimal. They did not elaborate on this difference. It is known from the Bayesian decision theory that the loss function for which the 'expected state of the atmosphere' is optimal, is a representation of the root mean squared error (RMSE) verification method. These results are re-visited in Section 2. The paper further shows how an ambiguous QuikScat observation can be represented as a 'super-observation' if we use the 'expected state of the atmosphere' as the 'analysis'. The implementation of the 100-km resolution ambiguous QuikScat observations in HIRLAM 3-D-Var (3-dimensional variational analysis) (Gustafsson et al., 2001) at the Norwegian Meteorological Institute is presented in Section 3. Examples of super-observations are shown in Section 4, along with a comparison of the RMSE verification score from a forecast trial that uses QuikScat observations according to this implementation and scores from a trial that do not use these observations. Forecasts of two storms in the verification period are also compared.

Bayesian decision theory applied to numerical weather prediction
The analysis vector, X a , is in this context that the state of the atmosphere that is used as a starting point for the integration of the future state of the atmosphere in a forecast model. An analysis is determined using observations and a short forecast, also known as the first guess or background. If the first guess, X b , and the observations, Y, have Gaussian error characteristics and no bias, and if the observation operator H is linear, the 'best' analysis (Lorenc, 1986) is given by where R is the observation error covariance matrix and B is the background error covariance matrix. This analysis can be derived by seeking the 'most probable state of the atmosphere' according to the available model for the probability density of the true state of the atmosphere given the observations and first guess, P[X|Y ∩ X b ]. Note here that we are considering X, Y and X b to be random variables. MSE verification methods are commonly used to verify numerical weather prediction systems. These methods typically calculate where Y i is an independent reference observation (of a given type), H i (X f ) is the model forecast for the reference observation and m is the number of reference observations. The MSE penalty is calculated separately for different reference observation types, and a numerical weather prediction system is considered to be 'better' if it gives an overall lower MSE penalty for all reference observation types. Note that these results are not changed if we compare the square root of the MSE penalty function, also referred to as RMSE verification. For short forecasts, the MSE verification methods provide a measure of the quality of the analysis. In the Bayesian decision theory, the MSE verification method for short forecasts, linear reference observation operators and un-specified reference observation types, can be represented by the quadratic loss function, l[X, X a ] = (X − X a ) T (X − X a ), where X represents the unknown true state of the atmosphere. The corresponding Bayes risk is given by A small Bayes risk corresponds to a good MSE verification score. An important result from the Bayesian decision theory (Ingleby and Lorenc, 1993;Rice, 1995) is that the X * * a function that minimizes this particular Bayes risk, is given by which is 'the expected state of the atmosphere'.

Assimilating ambiguous observations in HIRLAM 3-D-Var
where the index i denotes a single observation, and r is the number of observations. Adding a new observation denoted by r + 1, is then simply done by adding a contribution J o,r +1 [X] = −log [P[Y r +1 |X ∩ X b ]] to the cost function. The cost function is minimized by a general mathematical algorithm, and the starting point for the minimization is the first guess, X b . The HIRLAM 3-D-Var minimization algorithm (conjugate gradients) works best when the cost function is quadratic and does not change the shape from one iteration to the next.

The Gaussian analysis
Let us assume that we have a set of observations with Gaussian error characteristics, Y, and another set of independent non-Gaussian observations,Ỹ. Using Bayes rule, we may write the probability density function for the state of the atmosphere as The probability density function for observations with linear observation operators and Gaussian error characteristics may be modelled according to In HIRLAM 3-D-Var, it is assumed that the background error is Gaussian, giving Next, let us define a 'Gaussian analysis' by By comparing with eq. (1), we observe thatX b is the analysis which we would have if we only assimilated the first guess, X b , and the Gaussian observations, Y. It can be shown, by using Bayes rule and some algebra, that Further, since our observation errors are independent of each other and the first guess error, we have that We may now rewrite eq. (4), using the Gaussian analysis, where we see that the normalization term now has been written as P [Ỹ|X b ]. We observe that the probability density function where the ef-fects of the Gaussian observations, Y, have been included in the Gaussian analysis,X b .

Assimilating a single ambiguous observation
An ambiguous scatterometer observation consists of n ambiguities,Ỹ = {Y 1 , Y 2 , . . . , Y n }. For QuikScat, n is 2, 3 or 4, but in principle it may be any positive integer. The probability density function associated with a single ambiguous observation can be modelled according to whereγ represents the a priori probability for 'gross error' (within a given range),q k represents the a priori probability that ambiguity k is 'correct', Y k is the value of the kth wind vector ambiguity, H is the (linearized) observation operator (which calculates the 10-m wind in the target area from the model state) andR is the observation error covariance matrix for the 'correct' ambiguity. The a priori probability that ambiguity k is correct is available from the retrieval algorithm that estimates the ambiguities from the backscatter measurements. These a priori probabilities are based on how well the simulated backscatter (derived for each ambiguity), agrees with the actual backscatter observation. We assume that the observation error is independent of the first guess error. Using the Gaussian analysis and our non-Gaussian scatterometer observation along with the Bayes rule, we may write the probability density function for the true atmospheric state as We may further write, This is the 'expected state of the atmosphere' analysis when we are assimilating one ambiguous QuikScat observation. Note how that the Gaussian analysis in eq. (8) can be considered to be another ambiguity with a probability given by gross error.

The super-observation approach
Let us assume that we know the a posteriori probabilities that each ambiguity is correct, q k [X b ]. We may then define a superobservation according tõ If we assimilate the superobservation as a conventional wind observation with Gaussian error characteristics given byR and in accordance to eq. (1), we get which is identical to eq. (8), i.e. the 'expected state of the atmosphere' analysis when we are assimilating one ambiguous QuikScat observation.
We could say that our ambiguous QuikScat observation contains less information than a conventional wind observation, since additional information (X b ) is needed to estimate the superobservation, which can be assimilated in the same way as a conventional wind observation. As a consequence, the error in the superobservation is correlated with the error in the additional information used. However, assimilating the superobservation (Ỹ s ) together with the additional information (X b ) will still give an optimal analysis, according to the theory above. If we instead had calculated the superobservations (Ỹ s ) using the first guess (X b ), and then assimilated the superobservation (Ỹ s ) and first guess (X b ) together with Gaussian observations (Y), the first guess would have been given a greater weight relative to the Gaussian observations due to the correlation between the superobservation error and the first guess error. By using the Gaussian analysis to calculate the superobservations, we do not upset the balance between the first guess and the Gaussian observations. Also, the quality of the superobservation depends on the quality of the additional information used. The Gaussian analysis should in theory have a smaller error than the first guess, and should therefore yield higher quality superobservations.
The a priori probabilities,q k , may peak in a range of values, while the a posteriori probabilities, q k [X b ], calculated according to eq. (7) usually are ∼0% or ∼100%. This is illustrated in Fig. 1, which shows an example of the probability distribution for a priori probabilities, and the corresponding a posteriori probability distribution. It follows that the superobservation in most cases is almost identical to one of the ambiguities. Note that when the wind speed is very low, the ambiguities are relatively close to each other compared to the Gaussian analysis error, and the a posteriori probabilities will in these cases remain close to the a priori probabilities. We could say that the Gaussian analysis contains little additional information when the wind speed is very low. It should also be noted that the superobservation combines two sources of information, namely the QuikScat ambiguous observation and a Gaussian analysis based on a weather prediction model. If the Gaussian analysis wind field is changed, by for instance changing the forecast model, a superobservation in an area where there is a large change in the Gaussian analysis may align itself with a completely different ambiguity. Obviously, the most accurate Gaussian analysis field will yield the 'best' superobservation.
As other ambiguous observations are added to our problem, the analysis becomes increasingly complex since the QuikScat observations themselves have correlated errors. A correct and complex analysis would have ensured that observations with correlated errors received lesser weight than similar observations with uncorrelated errors. We approximate this error correlation effect by reducing the weight of the superobservations (by increasingR), and otherwise treating the superobservations as having independent errors. In HIRLAM 3-D-Var, the following cost function contribution is added for each QuikScat observation, where the index k refers to superobservation number k. The (reduced weight) observation error covariance used to calculate the cost function contribution in eq. (10) is (2.5 m s) −2 on the diagonal, and zero otherwise. The superobservation is initially calculated using the first guess and an a priori gross error ofγ = 0.1. If the resulting superobservation is closest to the first guess (in wind vector space), the observation is rejected, otherwise the superobservation is re-calculated, this time using zero gross error. Figure 2 shows the distribution of the difference between the resulting superobservation and the first guess wind speed for a four-month long period (using 400 000 QuikScat observations). The average bias between the superobservation and the first guess wind speed is 0.8 m s −1 . Note the large 'tail' where the superobservations have higher wind speed than the first guess. The large 'tail' is interpreted as the effects of rain contamination (which gives a stronger radar backscatter resulting in an overestimation of the wind speed). The superobservations with wind speeds more than 2.0 m s −1 greater than the first guess wind speed are discarded, so that the resulting bias comes down to 0.07 m s −1 . Note that this is a very crude and restrictive quality control method that will cause rejection of some 'correct' observations with potentially high information content. About 25% of the observations are discarded on average due to the gross error check and the biascorrection procedure. About 1% of the QuikScat observations are discarded in the 100-km thinning procedure.
There are always 80 iterations in the minimization algorithm. The QuikScat observations are ignored for the first 20 iterations. At iteration number 20, it is assumed that the control vector has approached the Gaussian analysis (i.e. the analysis we would have if we ignored all QuikScat observations). The superobservations are then calculated according to eq. (10) (using the control vector) and their cost contributions are from now on added to the global cost function (no other changes are made to the assimilation system).
The Norwegian Meteorological Institute uses a fixed Gaussian innovation covariance matrix (HBH T +R), with (2.0 m s) −2 on the diagonal and zero otherwise, to calculate the a posteriori probabilities (in eqs 6 and 7). Note thatB (which is never calculated explicitly) is implicitly defined by the global cost function (before the contributions from the superobservations are added). The superobservations are re-calculated every five iterations using the current state of the control vector. The matrix HBH T +R should ideally be re-estimated every time the superobservations are calculated so that it reflects the covariance in the differences between the control vector (Gaussian analysis) and the 'correct ambiguity'. It turns out that the a posteriori probabilities are not too sensitive to HBH T +R. There are for instance almost no visible changes in plots of superobservations, if the diagonal elements in HBH T +R are increased from (2.0 m s) −2 to (2.5 m s) −2 . Since the a posteriori probabilities are not sensitive to HBH T +R, this matrix is not changed during the minimization.
Variational quality control (Andersson and Järvinen, 1999) is activated from iteration number 40. In this setup, the superobservations will influence on how the variational quality control rejects observations from other data sources. However, variational quality control has not been applied to the Gaussian analysis before it is used to calculate the QuikScat superobservations in the HIRLAM 3-D-Var implementation, and it is therefore not clear if the quality of these superobservations is better than those calculated from the first guess. In the end, it turns out that there is no significant difference in the impact from superobservations calculated using the Gaussian analysis compared to using the first guess (which is how superobservations are calculated when they are assimilated operationally). With the introduction of 'outer loops' (Courtier et al., 1994) in future releases of HIRLAM, a Gaussian analysis that has undergone variational quality control will be available. However, if one uses this Gaussian analysis to calculate the superobservation, the variational quality control for the conventional data in the first loop will not take the QuikScat observations into consideration. It is not clear at this point how QuikScat observations best should be combined with variational quality control in an outer loop implementation of HIRLAM 3-D-Var. Note that the current HIRLAM 3-D-Var implementation of QuikScat observations violates the assumption in the mathematical minimization algorithm that the cost function is quadratic, but apparently not enough to cause any serious problems.

Examples
The superobservation approach has been applied to the 25-km resolution ambiguous QuikScat observations from NOAA/NESDIS in a now-casting product at the Norwegian Meteorological Institute. The 25-km resolution QuikScat observations in the NOAA/NESDIS BUFR product come with a NCEP (National Center for Environmental Prediction) model wind and an ambiguity selection. The QuikScat superobservations are calculated independently of each other, using the NCEP model wind as the first guess and no other Gaussian observations. In this setup, the Gaussian analysis,X b , is equal to the first guess, X b . The gross error term is set to zero. Observations flagged as rain-contaminated by the NOAA/NESDIS QuikScat processing software or by a MLE rain-contamination check designed by KNMI, are rejected. Superobservations that are closer to the first guess than any of the ambiguities are also removed. Figure 3 shows the NCEP model first guess wind for a weak low-pressure centre in the Atlantic Ocean, the ambiguities selected by NOAA/NESDIS and the corresponding 25-km resolution QuikScat superobservations. The QuikScat superobservations and the NOAA/NESDIS selected ambiguities are in close agreement, and suggest a more easterly location of the low-pressure centre. Figure 4 shows the NOAA/NESDIS selected ambiguities and the 25-km resolution superobservations for another case in the North Atlantic Ocean. Note that the NOAA/NESDIS selected ambiguities are more spatially consistent than the superobservations. This is what we would expect since the superobservation analysis is done separately for each observation, while the procedure behind the NOAA/NESDIS selection takes spatial consistency into consideration in the analysis so that the selected neighbouring winds have similar wind direction. Figure 5 shows another case where the superobservations provide a wind field in agreement with independent ship observations. The corresponding NOAA/NESDIS selected  ambiguities have here supported each other in such a way that the resulting selection is spatially consistent, but completely wrong.

Verification results
An impact study was performed over four months from November 2003 until February 2004. Version 6.2.1 of the HIRLAM forecast model was used. The model used digitalfiltering initialization, semi-Lagrangian advection scheme and physical parametrizations according to 'Savijärvi radiation', 'STRACO condensation', 'CBR turbulence scheme' and 'ISBA surface scheme' (Undén et al., 2002). Frames were taken from the ECMWF global model as lateral boundaries. The resolution was 20 km. Version 6.2.1 of the HIRVDA (HIRLAM Variational Data Assimilation) was used to assimilate observations in a 3-hr assimilation cycle with first guess at appropriate time (FGAT) enabled. The reference trial assimilated QuikScat observations along with ATOVS (Advanced TIROS Operational Vertical Sounder data) and conventional observations. The experiment trial was identical to the reference trial, except that it did not use the QuikScat observations. There was no observation cut-off (QuikScat data typically has a 2-hr time delay), and typically 500 QuikScat observations (100-km resolution) were available for an assimilation cycle. Note that the geometry of the QuikScat polar orbit causes the observation density to increase with latitude. Figure 6 shows the geographical distribution of the standard deviation in the MSLP (mean sea level pressure) analysis between the experiment and reference trials, for the model area. Note that there are no visible differences over central Europe, where there are many conventional observations (and no QuikScat observations). Figure 7 shows the corresponding   standard deviation in the MSLP 48-hr forecast. The differences in central Europe originate from perturbations in weather systems that are advected inland from the ocean during the forecast period. Figure 8 shows an example of a storm that moved in the Barents Sea along the coast of Norway, before turning inland at the Kola Peninsula. The figure shows 24-hr forecast from the reference and experiment (no QuikScat) trials, and the experiment analysis at the relevant time, just before the storm made landfall. We observe that the low in the experiment analysis is located in between the low in the experiment forecast and the low in the reference forecast. The relevant reference analysis was not too different from the experiment analysis (and in better agreement with the reference forecast, as we could expect).  Figure 9 shows an example of a storm that passed over Denmark on 14 January 2004. We observe that the experiment analysis low in this case is located closer to the reference trial than the experiment (no QuikScat) trial. Note that the pressure in a storm centre usually increases when it makes landfall.
The RMSE verification method used in this paper is based on the following penalty function, where m is the number of new independent verification observations (as suggested by the European working group on limited area modelling, for details on observations in the EWGLAM list, see for instance Vedel and Huang, 2004), i denotes each new independent observation (which has not yet been used in the assimilation), H r,i is the forward operator for each new independent observation and X[t] is a +t-hr forecast for the state of the atmosphere which is valid for the observation time. Only forecasts that start at 00 hr and observations at 00, 06, 12 and 18 hr are used in the verification. The different parameters (pressure, wind speed, etc.) are verified separately. The Rms is calculated and compared for the two trials. Note that the trial that gives lowest Rms penalty (wrt exactly the same, independent observations), is better according to this verification method. Figure 10 shows the standard deviation (Std) and bias (Bias) for MSLP for both the reference (with QuikScat) and experiment (no QuikScat) for the whole trial period (four months). Note that Rms[t] = √ Std[t] 2 + Bias[t] 2 . We observe that the trial where QuikScat observations were assimilated has a lower error (RMSE) than the trial without QuikScat observations. Figure 11 gives the same impression for the 10-m wind (FF10). The QuikScat observations made no significant impact on   the other surface parameters or parameters higher up in the atmosphere.
By studying variations in verification statistics for sub-periods of the four-month long verification period, we get an impression of how much we may expect the impact of assimilating QuikScat observations to vary. Figure 12 shows the daily contribution to the Rms ( t Rms[t]) for the MSLP, and it is clear that the positive verification results accumulated during a period from the middle of December until the middle of February.
Sinceq k is the a priori probability for each ambiguity, and q k is the corresponding a posteriori probability (where we have taken the Gaussian analysis into consideration), we should have thatq k E[q k ] = E[E[q k ]] q k . We may use this relationship to check if we are using an innovation model which is consistent with the data. Table 1 gives the average a priori and a posteriori values for the probability that each ambiguity is correct, for a typical assimilation cycle. There is a reasonable agreement between the average a priori and a posteriori probabilities, indicating that the innovation model is acceptable.

Conclusions
When the purpose of the analysis is to find the 'expected state of the atmosphere', the ambiguous QuikScat observations may be presented to the analysis system as superobservations with Gaussian error characteristics. This superobservation approach has been implemented at the Norwegian Meteorological Institute in HIRLAM 3-D-Var. Assimilating QuikScat observations according to this approach gave an overall small positive impact in MSLP and FF10 for a four-month period during the 2003/2004 Winter.

Acknowledgments
This work was sponsored by the Norwegian Space Centre. The author is also grateful for the support from Dr. Ole Vignes, John Bjørnar Bremnes and suggestions contributed by Lars-Anders Breivik. The perspectives and comments from the two anonymous reviewers are also greatly appreciated.