In order to run a regional Numerical Weather Prediction (NWP) model, an initial condition and a host model providing lateral boundaries are required. In this work we have used the HIRLAM model in which the initial condition is determined by combining observations and Fa short forecast with variational data assimilation, and the host model is usually the operational ECMWF model.
The HIRLAM data-assimilation system is based on the variational technique described in Gustafsson et al. (2001) and Lindskog et al. (2001). In this work it is run in four-dimensional variational analysis (4D-Var) mode, described in Huang et al. (2002) and Gustafsson (2006), and the background error covariances include a statistical balance, Berre (2000).
The variational procedure finds the optimal model state constrained by the models internal balances and observations, i.e. the host model is not taken into consideration. The aim of this paper is to examine a way to include selected information from the host model in the data-assimilation procedure.
One reason for doing this is to make the initial model state more consistent with the host model that will be applied at the lateral boundaries during time-integration. Another reason is that the global models often have more advanced assimilation techniques that use more satellite data, for example, and therefore give a better description of the large-scale flow. Using observations inside a regional domain only, larger scales cannot be described in a proper way, and will be aliased on smaller scales, so by blending the large-scale information into the regional model we can also hope to improve the regional model forecasts.
In Guidard and Fischer (2008), information from the global Action de Recherche Petite Echelle et Grand Echelle (ARPEGE) model analysis was blended into the regional Aire Limitée Adaption Dynamique développement InterNational (ALADIN) model by adding the ARPEGE analysis as a constraint in the data-assimilation. All upper air variables and surface pressure from ARPEGE were applied with a diagonal error matrix, and the spectral information vector was truncated such that only the large scales were affected.
In this work we take a different approach by selecting only the vorticity field from the host model, ECMWF, and we let that represent the large-scale flow which we assume is better described due to its extensive use of satellite observations and due its use of a global domain. Furthermore, we use a short forecast from ECMWF, instead of an analysis, to reduce the possibility of error correlations between the host model field and the observations used in the regional model data-assimilation. Also, it is not even possible for us in practice to use the ECMWF analysis in an operational context since it is not available in real time. Finally, the error covariances of the ECMWF vorticity field are not assumed to be represented by a diagonal matrix.
The paper outline is briefly as follows:
Details of the background error constraint, Jb, in HIRLAM and how vorticity is treated is described in Section 2. This information is needed for Section 3, where the formulation of the host model constraint is explained in detail.
Thereafter, in Section 4, we explain how the errors of the host model constraint are modelled and compared with the background error covariances. Section 5 points out some details on the implementation. The impact of the new constraint on an analysis is also investigated with a simple three-dimensional variational analysis (3D-Var) single case test.
The impact on forecasts, using a multi incremental 4D-Var setting, is then shown in Section 6. In Section 7 we summarise and make concluding remarks.
The following notations will be used throughout the rest of this article:
The model state increment vector in the HIRLAM assimilation is
where B contains the error covariances of . All variables in are related to each other and therefore the B matrix becomes quite complex. To ensure a fast convergence of the minimisation, a pre-conditioning through a variable transform is introduced. The change of variable is designed as a series of operators that successively removes the correlations between the components in until we end up with a new variable, called , whose covariance matrix can be assumed to be the identity matrix.
The series of transforms can be written as
G=calculate vorticity and divergence.
σb=background error standard deviation. Horizontal average, i.e. one value per vertical level.
λb=horizontal spectral density of the background error correlation. Diagonal matrix.
Vb=eigenvectors of the vertical background error correlation matrix. Orthonormal matrix, i.e. .
γb=eigenvalues of the vertical background error correlation matrix. Diagonal matrix.
After the G operator is applied, the state vector contains (, , , , ). Then the F operator removes the cross-correlations between the variables and the state vector is written (, , , , ), where subscript u means that the variables are unbalanced. and are often combined into one variable, and the balance operators described as going from the unbalanced state to the model state. Following the notations in Berre (2000), and also omitting all δ and the ln for pressure to make the expression more compact, the balance operators acts as the following
where E is the identity matrix, M, N, P, Q, R and S are vertical balance operators derived by linear regression and H is a similar horizontal balance operator that relates a linearised geopotential with vorticity, including also a horizontally varying Coriolis parameter. It should also be noted that since F−1 is triangular, the inverse is also triangular, so bothF−1andFhave the same shape. This means that when we go from model state to the control vector, or from the control vector to the model state as in eq. (4), vorticity is not affected by the balance operators. If we only consider vorticity, the transform from model state to the control vector can be written as
We believe that the host model is better at describing the large-scale flow of which vorticity is a good descriptor. From the previous section it is also clear that vorticity has very appealing characteristics because it is also a data-assimilation control variable. Therefore we add an extra term to the cost-function that describes the difference between the regional- and host-model vorticity fields:
is represented by a short forecast from the host model and Bls contains the error covariances of . For simplicity we will ignore the cross-correlations between the errors of and as well as the cross-correlations between the errors of and the observation errors. The cross-correlations between the errors of and are certainly different from zero and will be examined in the next section. Using a host model forecast instead of a host model analysis eliminates the need to consider the cross-correlations between the errors in and the observation errors, since the observations used to determine the analysis from which originates are at least 6 h old. Also, most regional models do not have the global host model analysis available in real time.
Now we set in eq. (7):
where is the regional model background.
We now also denote which gives
Bls is not a diagonal matrix which makes this formulation difficult to use. Therefore a pre-conditioning is introduced using a similar transform as for the background error covariances B, but using the error characteristics of Bls:
Projecting both and through Lls will give a state vector whose error covariance matrix is the identity matrix.
This pre-conditioning hopefully also helps assuring a good minimisation towards the host model vorticity field.
As a last step, eq. (6) is used to express Jk as a function of the control vector.
The derivative of with respect to gives a gradient that can be directly added to the total gradient, that now becomes
To implement the host model constraint Jk we need the error covariances of the ECMWF model as seen in the HIRLAM model geometry. To achieve this, ECMWF forecasts were interpolated to the HIRLAM geometry and then processed with the software used to determine the HIRLAM background error statistics. Both HIRLAM and interpolated ECMWF fields were processed simultaneously which allowed for cross-covariances between the errors of two models to be calculated.
The statistics presented will focus on displaying contents of importance for this work. For example, the contents of the and Lls operators that will be used in the assimilation, eqs. (6) and (10).
All calculations are carried out on the HIRLAM reference area, called Regular Cycle with the Reference (RCR), Fig. 1. The horizontal grid is expressed in a rotated lat/lon geometry defined by shifting the South pole. Model top is at 10 hPa and there are 60 vertical levels. Model levels are numbered from the top and down, i.e. model level 1 is at the top and level 60 is closest to ground. The data-assimilation calculations are carried out in spectral space and therefore an area extension, following Haugen and Machenhauer (1993), is used to ensure periodic variations in both horizontal dimensions. The horizontal resolution and number of gridpoints are presented in Fig. 2.
On each vertical level, an arbitrary variable, a(x,y), is represented by a bi-Fourier series:
So the variable a is in spectral space defined by the complex coefficients. An elliptical truncation
is used to obtain an isotropic and homogeneous resolution over the whole area. Each wave-number pair (k, l) is linked to a total wave-number k* by
where Ns is a scaling factor that is equal to the length of the extended domain divided by a map factor.
The ECMWF data used in the forecast experiments in this work originate from a spectral T799 model with 91 vertical levels and the model top at 0.01 hPa. T799 corresponds to a gridpoint resolution of about 25 km in a reduced Gaussian grid.
It should be noted that the ECMWF forecasts used for the derivation of the host model error statistics in this study were retrieved from the ECMWF archives in a horizontal resolution of 0.5 latitude–longitude. These forecast fields were subject to a bi-linear horizontal interpolation followed by a vertical interpolation to the HIRLAM geometry. Due to the simple interpolation procedure, one should not trust statistical information concerning horizontal waves shorter than approximately 150 km.
Short-range forecast errors were simulated by applying an ad hoc assumption that differences between forecasts of different length valid at the same time can be taken as a proxy for such forecast errors. Thus, we applied what is referred to as the ‘NMC method’. Pairs of +36 h and +12 h forecasts, valid at the same time, were obtained from the ECMWF and HIRLAM RCR archives for 87 situations, all at 00UTC, during the period September–November 2008. Forecast differences were relaxed towards zero in the lateral boundary relaxation zone of the forecast model (Fig. 2), in order to avoid aliasing effects through the extrapolations to the area extension zone that are applied in order to get bi-periodic variations over the extended domain. All calculations of statistics were carried out in spectral space ignoring covariances between spectral component coefficients of different horizontal wave-numbers, thus implicitly assuming horizontal homogeneity with respect to the covariances in gridpoint space.
The horizontal variance power spectra for vorticity, Pz(k*), is an output entity from the NMC statistics software that shows the contribution to the total variance for different wavelengths. Pz(k*) is plotted on three different vertical levels in the left column of Fig. 3, while the cross correlation between the HIRLAM and ECMWF vorticity fields is plotted in the right column of the same figure.
We see that HIRLAM has larger values throughout the whole spectrum. Since the NMC statistics are based on forecast differences, one reason could be that HIRLAM is more active and has a larger spread due to its higher resolution. Another reason may be difference between the data-assimilation methods applied in the forecasting systems producing the forecast fields that enter into the statistical calculations. ECMWF applied 4D-Var, while the HIRLAM fields during autumn 2008 were produced with 3D-Var. It is known that introduction of 4D-Var reduces the ‘jumpiness’ in the time-series of forecast evolutions, and the difference between +12 h and +36 h forecasts valid at the same time is essentially a measure of jumpiness. Also note that the horizontal resolution used in the retrieval of the ECMWF forecast fields from the archives makes any comparison of covariance spectra meaningless for horizontal wave-numbers greater than 64 (wave length 160 km).
The cross-correlation between ECMWF and HIRLAM vorticity forecast differences is strongest for the longest scales, right column in Fig. 3. Since longer scales are more predictable, and since these larger scales are to a large extent controlled by the host model through the lateral boundary relaxation, it is likely that both models produce somewhat similar large-scale structures which make the shape of the statistics reasonable. It seems like it is safe to ignore the cross-correlations for shorter scales while it may be more questionable for longer scales.
The horizontally averaged standard deviations, Fig. 4, also show that the HIRLAM values exceed the ECMWF ones on all levels. They both have rather similar shapes from the top and down to level 30 with a maximum around model level 15 (≈350 hPa) which corresponds to the jet-stream level. Below, HIRLAM has another maximum around level 43 (≈850 hPa) which is not present in the ECMWF statistics. These horizontally averaged standard deviations are applied in the assimilation weights and Lls but are adjusted with an empirical scaling factor. The scaling factor for σls is used to adjust the weight of the host model term, see Section 5.2.
Figure 5 shows the horizontal spectral densities, , as calculated by the NMC statistics software. and Ps(k*) are related via a simple expression:
is the input to the HIRLAM data-assimilation, but before applied in assimilation they are scaled such that their sum, divided with the number of gridpoints, is equal to 1. It should be noted that the result of that scaling is resolution dependent. So what we see in Fig. 5 is not exactly what is used in the Lls and operators, i.e. and . In Fig. 5 we see that the ECMWF spectrum is sharper than HIRLAM on all levels. At the same time, the ECMWF values are larger for large scales (small wave-numbers) and vice versa for short scales. For ECMWF in relation to HIRLAM, this means that large scales contribute more than short scales to the horizontal correlation which gives broader horizontal structure functions for ECMWF (not shown). This is natural due to the difference in resolution.
The eigenvector and eigenvalue for mode n and wave-number k* can be denoted . z=1, Number of vertical LEVels (NLEV) and γn,k. The sum of eigenvalues over all modes, for a specific wave-number k*, is constant and equal to the number of modes, which is the same as the number of vertical levels, thus:and the ratiotells us how much of the total variance that is explained by mode n. How the variance is distributed vertically for eigenmode n is shown by the eigenvector .
Figure 6 shows the eigenvalues for the first six leading modes. At the bottom, the combined effect of these modes is shown. We see that these modes include 90% of the variance for the longest scales and then gradually decrease as the scales get shorter. This is important to remember for the ECMWF data because we choose to use only the six first, or leading, modes of Vis and in the operator Lls, eq. (10). Lls then acts as a vertical filter that removes part of the small scale variance from the host model, including also vertically small scale noise introduced by the vertical interpolation procedure. In the same figure we also see how the higher order modes contain less large-scale information. Fig. 7 shows the structures of some eigenvectors. For each mode, the wave-number that is associated with the largest eigenvalue is chosen. The shapes are quite similar with a slight shift of phase for modes 5–8. It should however be mentioned that the shape of the eigenvectors for a specific mode varies a little between the wave-numbers.
A short ECMWF forecast, with the geometry described in the first paragraph of Section 4.2, is interpolated to the HIRLAM geometry described in Section 4.1. The interpolated ECMWF wind-field is thereafter transformed to spectral space. It is important to mention this because it means that a 25 km grid is interpolated to a 16 km grid, and the ECMWF field is thus oversampled and contains false high-resolution information.
In practice, however, this is never a problem for us due to two reasons:
First, the actual minimisation is a carried out on a grid mesh with lower resolution than the original grid (16 km), this incremental formulation is also described in Gustafsson et al. (2001). One of the reasons for this is to reduce the computational cost. In this work, all analyzes are carried out in a resolution either three or six times lower than the original resolution, i.e. a grid mesh of 48 km or 96 km. This means that the oversampled part of the ECMWF vorticity field is never used.
Second, The Jk term itself is truncated to an even lower resolution as described in Section 5.2.
In Section 4.3, a scaling of the standard deviation σls in the operator Lls [eq. (10)] was mentioned. The purpose of this scaling is to adjust the weight for the Jk-constraint in relation to the other terms in the cost-function. In the HIRLAM variational analysis code, a similar scaling is applied to the background error standard deviation to adjust the weight of the Jb term in relation the observation term Jo. The tuning factor, here called α, is applied as such:and the actual tuning procedure in this study is rather ad hoc. A few different values of scaling-factors, α=1, 5, 10, 100, were applied on a 3D-Var test case with J=Jb+Jo+Jk and the behaviour of the cost-function, and the individual terms separately, during minimisation, were studied. α=100 gave Jk no weight at all, while α=1 made Jk too dominating. With α=5 and α=10 the Jk and Jo terms both minimised smoothly. α=5 however gave a more distinct minimisation of the Jk term without losing the fit towards observations. Therefore α=5 was used throughout all experiments in this study.
The Jk term is also truncated at k*=45, wave length around 230 km, to make sure only large scales are penalised. In terms of low-resolution increments, described in Section 5.1, k*=45 corresponds to a resolution about seven times lower than the original resolution.
As a first step, we perform a 3D-Var analysis with only the background and Jk as constraints, i.e. J=Jb+Jk. The date of the single case test was 22 October 2008 at 06UTC, and both xb and are +06 h forecasts from the previous cycle.
The analysis converges in four iterations, not shown, and the resolution was 48 km, i.e. three times lower than the original resolution but the Jk term was truncated at k*=45 as mentioned in the previous section. Both the first guess xb and analysis xa are then compared with the host model field xls. This is done by computing Root Mean Square (RMS) values in gridpoint space for (xb–xls) and (xa–xls) for the upper air variables б, u, ⊼, T, q. The differences are computed for every vertical level and are presented as profiles in Fig. 8.
Because Jk contains vorticity only, with no cross-correlations to any other variables, the analysed vorticity field must be closer to ECMWF than the first guess, otherwise, there would be something wrong with the implementation.
We now see in the upper left plot in Fig. 8 that the analysed vorticity field is clearly closer to ECMWF compared with the first guess. So the Jk implementation seems to do what it is supposed to. Differences in the other variables are produced by HIRLAM's internal balances via the B matrix and do therefore not have to be closer to ECMWF.
It is also interesting to investigate how the host model constraint acts on different scales. For simplicity, we only look at the vorticity field and compare the complex spectral coefficients of the differences and . As explained in Section 4.1, these differences are represented as bi-Fourier series as shown in eq. (14), and each wave-number pair (k, l) is linked to the total wave-number k* [eq. (16)].
To be able to visualise these difference fields they have been expressed as 1-D arrays over k*. Each element is the average of the absolute values of all spectral coefficients with wave-number pairs (k, l) that are linked to the same k*.
In Fig. 9, the top figure shows the averaged 1-D spectrum of the difference as a full line and as a dotted line. It can now be seen that the analysed vorticity is drawn towards ECMWF on all scales where Jk is active.
The size of the increments is also shown in the middle of Fig. 9. This plot shows the average absolute value of the vorticity increments at the last iteration which corresponds to the difference between the full and dotted line in the top plot in Fig. 9.
The largest increments are introduced at k*=2–8, peaking at 5 and thereafter slowly decreases until k*=45 where it becomes zero.
However, the relative size of the increments is rather uniformly distributed over the spectrum as shown in the bottom plot of Fig. 9.
Now we have verified that the implementation of the Jk term draws the analysis closer to the host model (ECMWF) vorticity field, and that the largest increments are introduced on the longest scales.
A 4D-Var setup with two outer loops was chosen for the forecast impact experiment in order to be close to what is used operationally at SMHI.
Two outer loops in HIRLAM means that the inner loop begins in a coarse resolution and performs a number of iterations. Then, as an intermediate step, the control vector is saved and a temporary analysis file written from which a new, non-linear trajectory is calculated with the forecast model to make a re-linearization of the observation operators and for the tangent-linear model. The second inner loop then continues in a higher resolution to finish the analysis.
In HIRLAM the 4D-Var time-window is ±3 h around the nominal analysis time and the first outer loop uses the +03 h forecast from the previous cycle as first guess. Since the control vector is updated and the cost-function gradient is calculated at the beginning of the time-window it is appropriate to also use a +03 h forecast for the host model constraint.
The first loop is carried out with 96 km resolution (six times lower than original) for a maximum of 60 iterations, and the second loop is carried out in 48 km resolution for a maximum of 30 iterations.
The model was run up to +48 h every 6th hour from 24 January 00UTC to 25 February 2009 at 18UTC, i.e. 33 d. For simplicity, only conventional observations were used in the data-assimilation. Conventional observation types include SYNOP, SHIP, BOUY, TEMP, PILOT and AIRCRAFT.
Two sets of assimilation experiments were carried out throughout the 33-d period:
As in the 3D-Var singe case test, it is interesting to see the cost-function behaviour for a 4D-Var case with observations. One date in the impact experiment was chosen arbitrarily, 24 January 2009 at 12UTC, and studied. The cost-function contains three terms Jb, Jo and Jk and each part is plotted separately, together with the total value Jtot in Fig. 10.
It is important to note that the cost-function values of the Jb and Jk terms are resolution dependent. Jb is the scalar product of the control vector, , and similarly the Jk value is the scalar product described in eq. (12). These vectors are thus longer in the 2nd outer loop. At the same time, the horizontal spectral density, λb and λis, undergoes a resolution-dependent scaling described in Section 4.3, as well as the coefficients of the spectral fields. These three factors combined, cause the cost-function values of Jb and Jk to make a discontinuous jump between the 1st and 2nd outer loops. Because of this discontinuity, the decrease of the cost-function value should be compared with the initial value of each loop.
In Fig. 10 we see that the total cost-function value, as well as the observation term Jo, decreases and converges in both loops, while Jk decreases and converges, already in the 1st loop, and is then almost constant in the 2nd. It is possible that the reason for this is that Jk only acts on large scales which converges faster. The main conclusion from Fig. 10 is that the analysis draws closer to both observations and Jk simultaneously.
Root Mean Square differences for this particular case, as described in Section 5.3 and Fig. 8, also show that the analysed vorticity is clearly closer to ECMWF (not shown).
Forecasts were verified by computing mean and RMS errors against all available SYNOP and radiosonde stations in the area. In this case it adds up to around 2000 SYNOP and 175 radiosonde stations.
For surface parameters like 2 m temperature, humidity and 10 m wind the impact was neutral. Mean sea level pressure showed a clear positive impact, Fig. 11. Upper air parameters, computed against radiosondes, also had a clearly reduced RMS error for the jk_active experiment for temperature, wind, geopotential and specific humidity. Fig. 12 shows vertical profiles of mean and RMS errors of temperature calculated at midnight, 00UTC and noon, 12UTC.
Before making the concluding remarks, we also compare one case in the forecast experiment period and see what the differences looks like. The date was chosen from a time-series of RMS errors of +24 h forecasts of surface pressure, Fig. 13. Around 30 January 2009, the RMS errors of +24 h forecasts differed with approximately 1.0 hPa between the experiments, with advantage to the jk_active experiment. Therefore the 12UTC run from 29 January was chosen for the study since the +24 h forecast from that time is valid on 30 January. We are thus 5 d into the experiment period and the two experiments have evolved separately through a number of data-assimilation cycles and forecast integrations. One cannot expect large-scale differences only just because Jk acts on large scales. The non-linear forecast integration can generate differences on all scales.
Figure 14 shows the analyzes of the mean sea level pressure (MSLP) fields, from the experiment (jk_active) and reference (ref) runs, and the difference between them. The overall large-scale synoptic features are quite similar with an intense low at the coast of Newfoundland, and a high pressure system to the south-east, over the Atlantic. A large trough stretches from Greenland to the south-east with a low-pressure system just to the West of Ireland. There is also another low-pressure minima near Iceland. Then there is a massive high pressure over northern Russia with a ridge that stretches to the south-west into Scandinavia. The difference map reveals that the most noticeable differences are associated with the low-pressure system to the West of Ireland. A closer look at the MSLP fields shows that it is the shape of the low-pressure system that is different, and that the low is bit deeper in the experiment (jk_active).
When this is integrated forward in time we see that it is a blocking high pressure over Russia/Scandinavia, and the low to the West of Ireland moves straight to the North while the low outside of Newfoundland moves to the East. After 24 h, Fig. 15, the low originally to the West of Ireland is located near Iceland. The differences have grown in both spatial scale and magnitude and are still associated with the low-pressure system near Iceland. If we compare the MSLP fields we see that it is the shape and the depth of the low that have changed. The low in the reference run is now deeper and more intense than the low in the (jk_active) experiment. This low has thus developed differently in the experiments. In the reference run, the low moved northwards and intensified while the low in (jk_active) experiment moved to the North and weakened a bit.
In this paper, we have described and tested a method to assimilate the host model information (ECMWF) into a regional NWP model (HIRLAM). The vorticity field from the host model was chosen to act as a constraint in the regional model variational analysis. For simplicity, the host model vorticity constraint was assumed to be univariate, i.e. no cross-correlations to any regional model variables or observations. This, along with the characteristics of vorticity in the statistical balance, led to an error covariance matrix that could be made diagonal by using the horizontal spectral densities and vertical eigenmodes.
A forecast experiment of 33 d has given very promising results with reduced RMS errors of upper air forecast fields compared with radiosonde observations as well as mean sea level pressure forecast fields compared with synop observations.
The study clearly indicates that the host model information is important to include when determining the initial condition.
Constructive comments from Xiang-Yu (Hans) Huang, NCAR and one anonymous reviewer helped us to improve this manuscript and these contributions are acknowledged.
Gustafsson, N, Berre, L, Hörnquist, S, Huang, X.-Y, Lindskog, M and co-authors. 2001. Three-dimensional variational assimilation for a limited area model. Part I: general formulation and the background error constraint. Tellus. 53A, 425–446.
Lindskog, M, Gustafsson, N, Navascués, B, Mogensen, K. S, Huang, X.-Y and co-authors. 2001. Three-dimensional variational assimilation for a limited area model. Part II: observation handling and assimilation experiments. Tellus. 53A, 447–468.