1. Introduction, context
The invitation to tender for the C3S regional reanalysis projects requested that "consistent estimates of uncertainty shall be generated based on explicit statistical assumptions on model errors and observation errors". We have proposed a task as a part of building the CARRA (Copernicus Arctic Regional Reanalysis) reanalysis system to devise a method to provide such uncertainty estimates.
In the project proposal, we noted that to get perfect figures for reanalysis uncertainty statistics, it would be needed to know its deviations from the true atmospheric state. As all potential information sources to assess this contain errors, this true state is not directly available, and uncertainty is generally provided as statistical estimates with varying degree of reliability. Usually uncertainty estimates are based on some reference, which the reanalysis can be compared to or based on some proxy dataset with statistical properties assumed to be similar to that of the deviations from the true state. It should be noted that uncertainty estimation is an area where all methods have their weaknesses, and there are continuous developments in methods and discussions on how to present it to users, so the area of deriving and presenting uncertainties will for sure be subject to future evolution.
The project proposal contained a sketch for deriving reanalysis uncertainty information from the dataset also used to construct background error covariance statistics for the assimilation, which we will follow with more details in this report. The spread of the ensemble serves as proxy for background errors when used for generation of background error covariance, and will also be used as a basis for our uncertainty estimates. This is not a perfect proxy for the real uncertainties, but will describe uncertainties in a consistent way for the prognostic variables of the model at any location. It will provide information on the variation of uncertainty with height as well as physically consistent relations between the uncertainties of the different physical quantities. Since this approach has no guarantee that the ensemble perturbations have the correct magnitude, some scaling of this statistics will be necessary to go from ensemble statistics to reanalysis uncertainty. More details for the proposed method are outlined in the following, in particular how such a scaling factor, called α below, will be estimated.
The purpose of the reanalysis applications is to reconstruct the past evolution of the atmospheric flow striving for as high as possible level of spatial-temporal consistency of the data. In order to achieve this a latest state-of-the-art numerical weather prediction model is initialized with consistent records of observations collected over (a long) time period. The CARRA reanalysis data set aims to provide a high resolution reconstruction of the atmospheric flow over Arctic over the last 24 years. The CARRA reanalysis system is based on the HARMONIE-AROME convective-scale non-hydrostatic numerical weather prediction system, see details on the system in Yang et al (2020a) and how it has been verified and tested in Yang et al (2020b). The HARMONIE-AROME system is used for operational weather forecasting at several National Weather Centers in Europe, including Danish Meteorological Institute and MetCoOp, the operational collaboration between Norway, Sweden and Finland. The CARRA reanalysis system employs a 3D-Variational data assimilation scheme with three hour update frequency to constrain the model simulations with observations. The observations are combined optimally with the short-range high resolution background field taking into account uncertainty about the background fields and the observations. The analysis is released in a so-called cycled environment, where a short- range HARMONIE AROME forecast initiated from the produced analysis is used as a background field for the next analysis.
Because of the relatively sparse observation coverage over the Arctic, the background error statistics play an important role in constructing the reanalysis fields. The background error covariance model provides information on spatial scales of the short-range forecast error fields and on the balance relationships between different model state components. This information, together with observation error statistics, is used to filter out observation noise and to spread out the observed information spatially and from observed quantities to unobserved model state variables, in a way that physically consistent gridded field with realistic spatial variations would be obtained.
2. Approach for deriving uncertainty related information
The CARRA reanalysis is a deterministic numerical weather prediction system which is unable to incorporate an actual error-of-the-day estimate of background error statistics constructing the analysis, even though the true background error statistics are orography- and weather-situation dependent. Nevertheless, designing CARRA reanalysis system significant efforts were spent on the estimation of the climatological background error covariance. A model of the climatological background error covariance was formulated by Berre (2000). This model is based on spectral representation of spatial covariances and imposes severe assumptions of homogeneity and isotropy in order to make computations efficient and feasible. In other words, this model is able to provide an estimate of the background error covariance, which is valid on average, but cannot resolve spatial variation within the domain due to orography or temporal variations due to different weather regimes.
The parameters of the climatological background error covariance are derived from time-series of forecast differences which are used as a proxy of background errors. In order to obtain forecast differences an ensemble approach was employed, the forecast differences are constructed as differences between ensemble members and control. The ensemble members are produced using the so-called BRAND scheme in deterministic mode (randomization of the B-matrix) for generation of initial condition perturbations. The ensemble of short- range forecasts was aggregated for a period of three weeks. The first 10 days of the simulations was excluded in order to minimize impact of the spin-up issues on the estimate of the background error statistics. The background error covariance is estimated separately for winter and summer periods over the CARRA-East and CARRA-West domains in order to account for strongly pronounced seasonal variability of the forecast error growth over Arctic. 800 ensemble members seemed to be sufficient to obtain a statistically stable estimate of climatological background error covariance. The generation of the CARRA reanalysis ensemble was a time consuming procedure with high computational and manpower costs. To run a continuous ensemble CARRA system was unaffordable under present time constraints for the production of CARRA dates sets. As a compromise the CARRA ensemble was computed for several short time slots on both domains. We speculate that the climatological uncertainty of the reanalysis is mostly determined by the climatological uncertainty of the background field used for reanalysis. Thus for the estimation of the re analysis uncertainty we plan to reuse the same ensemble data sets that are used for the estimation of the background error statistics, namely aggregated summer, July 2012, and winter, January 2017, ensembles over CARRA- East and CARRA- West domains.
2.1. Framework for a method
The estimation of the reanalysis uncertainty may be deduced from comparison of the reanalysis data sets to observations taking into account statistical information on observation errors and short- range forecast errors. The data assimilation provides a framework for such comparison. Innovations, the observations minus background equivalents in observation space, naturally combines these two sources of information. In well-tuned systems, provided that background and observation errors are uncorrelated, the variance of innovations equals to the sum of observation error variance and variance of the background equivalents in observation space.
\[ E(d_{b}d_{b}^T) = R + P, \quad (1) \]where R is the observation error covariance, P is the covariance of model equivalents in observation space, \( d_{b} = y - h(x^b) \) denotes innovations, y are observations, \( h(∘) \) is observation operator, \( x^b \) is the background of the control run and \( E(∘) \) denotes mathematical expectation. Here we assume that innovations are unbiased, namely \( E(d)=0 \) Note that even if the true background error covariance is clearly flow-dependent and varies with time, the 3D Variational data assimilation assumes climatological background error covariance.
For the well-tuned ensemble a sample estimate can be used to compute P, the background error covariance in observation space. \( P^S ≈ Z Z^T \) , where \( Z = (Z_{1}, ..., Z_{i},..., Z_{K}) \) is the matrix of the normalized ensemble perturbations in observation space \( Z_{i}= \frac{1}{\sqrt{K}}(h(X_{i}^b)-h(x^b)) \) . Here, \( X_{i}^b \) denotes the background of the i-th ensemble members. K is the ensemble size.
The CARRA reanalysis ensemble used for estimating the background error statistics, however, is driven by the low-resolution ERA5 EDA on the lateral boundaries and contains only 9 members per time instance. It is unrealistic to expect a proper spread from this small sized ensemble. Therefore we assumed that \( P= \alpha P^S \) , where α is a tunable factor that can be estimated from equation (1) provided that we have a large enough sample to obtain a stable estimate of the innovation variance E( dbdbT) and realistic assumptions on observation error statistics R.
\[ E(d_{b}d_{b}^T) ≈ \alpha P^S + R= \alpha ZZ^T + R, \quad (2) \]Note that the screening and data assimilation procedures implemented in the HARMONIE-AROME forecasting system project the +3h forecast of each ensemble member and of the control run into observation space and store the observation-minus-forecast-differences in the corresponding ODB (Observation Database) CCMA database. Both screening and data assimilation were activated for every cycle in the BRAND scheme for each ensemble member and for the control. For example, \( Z_{i}= \frac{1}{\sqrt{K}}(d_{b} - d_{b}^i) \) , where \( d_{b}^i \) is the innovation corresponding to the i-th ensemble member
There are simple off-line procedures available in the ODB model code module to extract the needed quantities from the ODB.
As it is obvious from Eq. (2) the estimate of the tuning factor α depends on how realistic the assumed observation error statistics are. Observation error statistics are often tuned in order to obtain an optimal performance of the system for certain forecast ranges. The observation error statistics might be adjusted to guarantee certain weights for a particular component of the observing system. The HARMONIE-AROME system assumes uncorrelated observation error. For example, the observation error variance might be increased artificially in order to force the system to put less weight on certain components of observation system or make certain components be less dominating. This can be done as a compensation for actually correlated observation errors or if the certain observed quantities are overrepresented.
Therefore, we are here planning to estimate observation error variances for each observation type through posterior statistics diagnostics,
\[ R^⋆ = E(d_{a} d_{b}^T), \quad (3) \]where \( d_{a} = y - x_{a} \text{ and } x^a \) is the analysis.
Under assumptions of Gaussianity for observation and background errors the analysis can be expressed as
\[ x^a = x^b+ Kd_{b}, \text{ where } K= BH^T(HBH^T+ R)^{-1} = BH^T (P+R)^{-1} \]Here B denotes the background error covariance in the model space and H is the linearization of the observation operator. Thus we use a tangent linear approximation to compute the background error covariance in observation space, P = HBHT
Indeed
\[ E(d_{a} d_{b}^T)=E((y-h(x^a))(y-h(x^b))^T) = E(y-h(x^b+K(y-h(x^b)))(y-h(x^b))^T)≈ E((I- HK)(y-h(x_{b})) (y - h(x_{b}))^T) = (I-HK) E(d_{b} d_{b}^T) ≈(HBH^T + R -HBH^T) (HBH^T+R)^{-1}(HBH^T + R)=R \]
We propose a final estimate of the observation error variance to be an average between the assumed and the diagnosed observation error variance
The observation error variance will be estimated for each observation type separately including the dependency on vertical coordinate where this is applicable. Note that the 3D Variational data assimilation applies a climatological background error covariance. This means that the diagnostic equation (3) holds only on average. Diagnosing R one has to pool together innovations over a long enough period of time. We will use innovations pooled together from the whole data set used to derive the background error statistics and we will use only those observations which were accepted for assimilation by all ensemble members.
Let N be a total number of cases, the total number of aggregated data assimilation cycles. Lets pn to be a number of observations of a particular type assimilated during cycle = 1,...,N. Lets \( d_{b,n}^* = \frac{1}{\sqrt{R^*}}d_{b} \) to be a vector of standardized innovations of the particular observation type assimilated during cycle n. The standardization is performed in order to allow for dependency of the observation error variance on vertical coordinate. Then for this particular observation type, the estimate of the standardized innovation variance can be computed as follows.
\[ E(d_{b}R^{*-1}d_{b}^T) = \frac{1}{N} \sum_{N}^1 \frac{d_{b,n}^{*T}d_{b,n}^*}{p_{n}}, \quad (4) \]Note that this quantity should always exceed 1. Otherwise the observation error variance seems to be heavily overestimated and such observation type should be excluded from the computation of the tuning factor α . The ensemble estimate standardized background error covariance in observation space for a particular observation type P*s can be computed as follows
\[ P^{*s} = \frac{1}{N} \sum_{N}^1 \frac{1}{p_{n}} Tr(Z_{i,n}^*Z_{i,n}^{*T}), \quad (5) \]Here \( Z_{i,n}^* \) is the standardized normalized perturbation in the observation space of the particular observation type computed during cycle n
The inflation factor αp can be then computed from Equations (4) and (5)
\[ α_{p}=\frac{E(d_{b}R^{*-1}d_{b}^T)-1}{P^{*S}} \]The final inflation factor α will be obtained by optimally weighing together αp corresponding to different observation types. The weighting will reflect the importance of different observation types and the number of assimilated observations.
The final inflation factor α will be used to compute the representative ensemble of the main reanalysis variables
This representative ensemble will be used to deduce the time- and domain- averaged statistical estimate of the reanalysis uncertainty for the main reanalysis variables. The uncertainty can be presented in form of usual statistical quantities, such as for instance standard deviations, inter-quantiles ranges or box-and-whiskers diagrams.
In addition to statistical estimates of the uncertainty, the presented diagnostic analysis will outline level of consistency of the data assimilation scheme used for the reanalysis. Such an analysis will provide valuable information on how trustable the analysis system is.
3. Verification information as a complement
The analysis will be complemented with verification results that also give an impression on the uncertainty, i.e. higher errors give on average less trust in the reanalysis. To support the uncertainty estimates described above, we will prepare a verification summary to be included in an updated version of the user guide/documentation. The summary will focus on how errors for near-surface weather parameters (e.g. Mean Sea Level Pressure, 2m air temperature, 10m wind speed, relative humidity in 2m, precipitation) varies between parameters and with season, region and lead time. This will give additional useful information for many users of the re analysis.
4. Practical implementation and output to users
The CARRA reanalysis has a quite comprehensive list of output fields to users, see the Data User Guide (Nielsen et al, 2019). We note from the statistical approach for deriving uncertainty information described above, that we have the restriction that we can only assess the uncertainty of the physical variables available and stored within the 3D-Var system, which is a small subset of the full output list from the model system. Also, to obtain robust statistics of uncertainty related parameters from a limited size ensemble, we need a certain aggregation of statistics, so we expect for instance that local gridpoint-to-gridpoint variations will be noisy with variations not reflecting the real variations of uncertainty. In summary, we propose to provide uncertainty estimates for a subset of the full output list at a certain aggregation level to be described in the following.
As for what uncertainty measures to provide to users, we propose to provide estimates of standard deviations of errors. The framework we have will for instance not be suited for providing biases.
We will not provide uncertainty estimates for most of the single level parameters and soil level parameters listed in the output tables of Sections 3.1 and 3.2 of the User Guide (Nielsen et al, 2019), but will provide uncertainty information uniquely for:
- Mean sea level pressure (GRIB code 151)
Several of the remaining single level parameters in these tables will be assessed by verification information at observing stations, as outlined in Section 3 above, which we propose to include in the next update of the Data User Guide. The variables provided in this way and the breakdown on sub-regions/areas will be as in Tables 3.1 and 3.2 of the full system documentation (Yang et al, 2020a).
For model level variables (see Section 3.3 of the Data User Guide), we will provide uncertainty information as outlined in Section 2 above for:
- Specific humidity (GRIB code 133)
- Temperature (GRIB code 130)
- u component of wind (GRIB code 131)
- v component of wind (GRIB code 132)
We will also provide uncertainty information for these physical variables on pressure levels (a subset of the complete list of pressure level variables in Section 3.4 of the Data User Guide).
We will provide uncertainties as horizontally averaged vertical profiles for each physical parameter, not as full 3D fields due to the issues of statistical robustness mentioned above. These estimates will possibly be provided to the MARS archive (and the Climate Data Store) in the next update of the CARRA output for MARS storage. (MARS archiving is already running for our main output variables and there will anyway be a second step with updating of the MARS archiving and later data release, for the parameters which depend on the GRIB approval process.) The feasibility and effectiveness of storing on MARS will be further discussed with the ECMWF MARS team, and this discussion will form the basis for a proposed decision on in which form(s) the uncertainties will be presented to users.
In the derivation of the uncertainty, standard deviation of errors will be computed for both domains and both for winter and summer. We will assess the differences between domains and seasons to decide on whether there are significant and robust differences in the uncertainties. We will probably provide average values between the domains unless we identify large differences, and most likely provide average values over the year and not seasonal values.
The provided uncertainty statistics is available from Copernicus Arctic Regional Reanalysis (CARRA): known issues and uncertainty information#Uncertaintyinformation
5. References
Nielsen, K. P. Yang, X., Agersten, S., Støylen, E. and Schyberg, H., 2019: C3S Arctic regional reanalysis – Data User Guide. Copernicus Climate Change Service contract 2018/C3S_D322_Lot2_METNorway/SC2 Deliverable Report no C3S_D322_Lot2.1.1.1. Available at Copernicus Arctic Regional Reanalysis (CARRA): Data User Guide
Yang, X., Nielsen K.P., Dahlbom, M., Amstrup, B., Peralta, C., Høyer, J., Englyst, P.N., Schyberg, H., Homleid, M., Køltzow, M.A.Ø., Randriamampianina, R., Dahlgren, P, Støylen, E., Valkonen, T., Palmason, B., Thorsteinsson, S., Bojarova, J., Körnich, H., Lindskog, M., Box, J., Mankoff, K., 2020a: C3S Arctic regional reanalysis – Full system documentation. Copernicus Climate Change Service contract 2018/C3S_D322_Lot2_METNorway/SC2 Deliverable Report no C3S_D322_Lot2.1.1.1. Available at Copernicus Arctic Regional Reanalysis (CARRA): Full system documentation
Yang, X., Nielsen K.P., Dahlbom, M., Peralta, C., Schyberg, H., Homleid, M., Køltzow, M.A.Ø., Randriamampianina, R., Dahlgren, P., Vignes, O., Støylen, E., Valkonen, T., Bojarova, J., Lindskog, M., Hagelin, S., Körnich, H., Palmason, B., Thorsteinsson, S., 2020b: CARRA report on test and verification. Copernicus Climate Change Service contract 2018/C3S_D322_Lot2_METNorway/SC2 Deliverable Report no C3S_D322_Lot2.1.2.1. Available at Copernicus Arctic Regional Reanalysis (CARRA): Complete test and verification report on fully configured reanalysis and monitoring system
This document has been produced in the context of the Copernicus Climate Change Service (C3S).
The activities leading to these results have been contracted by the European Centre for Medium-Range Weather Forecasts, operator of C3S on behalf of the European Union (Delegation agreement signed on 11/11/2014). All information in this document is provided "as is" and no guarantee or warranty is given that the information is fit for any particular purpose.
The users thereof use the information at their sole risk and liability. For the avoidance of all doubt , the European Commission and the European Centre for Medium - Range Weather Forecasts have no liability in respect of this document, which is merely representing the author's view.
