Data assimilation of GNSS zenith total delays from a Nordic processing centre

Abstract

Atmospheric moisture-related information estimated from Global Navigation Satellite System (GNSS) ground-based receiver stations by the Nordic GNSS Analysis Centre (NGAA) have been used within a state-of-the-art kilometre-scale numerical weather prediction system. Different processing techniques have been implemented to derive the moisture-related GNSS information in the form of zenith total delays (ZTDs) and these are described and compared. In addition full-scale data assimilation and modelling experiments have been carried out to investigate the impact of utilizing moisture-related GNSS data from the NGAA processing centre on a numerical weather prediction (NWP) model initial state and on the ensuing forecast quality. The sensitivity of results to aspects of the data processing, station density, bias-correction and data assimilation have been investigated. Results show benefits to forecast quality when using GNSS ZTD as an additional observation type. The results also show a sensitivity to thinning distance applied for GNSS ZTD observations but not to modifications to the number of predictors used in the variational bias correction applied. In addition, it is demonstrated that the assimilation of GNSS ZTD can benefit from more general data assimilation enhancements and that there is an interaction of GNSS ZTD with other types of observations used in the data assimilation. Future plans include further investigation of optimal thinning distances and application of more advanced data assimilation techniques.

Introduction

Data assimilation in numerical weather prediction (NWP) optimally blends observations with an atmospheric model in order to obtain the spatial distribution of atmospheric variables and to produce the best possible model initial state. It was realized early that the forecast quality is strongly dependent on an accurate description of the initial state . There are strong requirements on the infrastructure and computing power for today's state-of-the-art high-resolution modelling systems. As model resolutions increase it is increasingly important to utilize observations with high spatial and temporal resolution to initialize mesoscale phenomena, such as convective storms and sea breezes.

The meteorological weather services of Sweden, Norway and Finland recently joined forces around a common operational kilometre-scale forecasting system named MetCoOp . The forecast model used within MetCoOp is developed in the framework of the shared Aire Limitée Adaptation dynamique Developpement InterNational (ALADIN) High Resolution Limited Area Model (HIRLAM) NWP system. This system can be run with different configurations and in MetCoOp the HIRLAM-ALADIN Regional Meso-scale Operational NWP In the Europe Application of Research to Operations at Mesoscale (HARMONIE-AROME) is used . The main components of the ALADIN-HIRLAM NWP system are surface data assimilation, upper-air data assimilation and the forecast model for the forward time integration.

The only direct humidity measurements used in the MetCoOp upper-air analysis are provided by vertical profile measurements from radiosondes. In addition, humidity-related information is provided by radar measurements , by satellite-based information and by moisture-related observations from the Global Navigation Satellite System (GNSS) zenith total delay (ZTD). Satellite observations are coupled to the moisture through the dependence of the radiative transfer at the top of the atmosphere on the atmospheric moisture distribution. The disadvantage of all of these humidity-related observation types, except GNSS ZTD, is that they are only available at particular times of the day (radiosonde and satellite measurements) or their availability is dependent on weather conditions (radar measurements). GNSS ZTD estimates, on the other hand, are available at all times with high temporal resolution (15 min), for all weather conditions. The ZTD is in fact an estimation, but for simplicity we hereafter refer to it as an observation. Moisture-related observations in the form of GNSS ZTD are a relatively new source of mesoscale atmospheric humidity information. ZTD observations obtained from the network of ground-based GNSS receivers contain horizontally dense information on the total columnar amount of water vapour (TCWV). A number of assimilation studies have shown a positive impact of GNSS ZTD observations on NWP systems at a horizontal model grid resolution on the order of 10 km . The importance of combining the GNSS data with other types of observations has been highlighted in several studies . Some encouraging results from assimilation of these observations at a kilometre-scale horizontal resolution have been obtained , and GNSS ZTD from 28 receiver stations are assimilated operationally in MetCoOp. These 28 receiver stations have been selected from the rather few receiver stations over the MetCoOp domain. Most of these so-called supersites, processed by several centres for comparison purposes. MetCoOp operationally uses data processed by the Met Office processing centre in the United Kingdom (METO) and by the Royal Observatory processing centre of Belgium (ROBH).

The EUMETNET GPS Water Vapour Program (E-GVAP) is a collaborative effort between the European geodetic community and several European national meteorological institutes. The purpose of E-GVAP is to provide atmospheric water vapour observations for use in operational meteorology. ZTD observations obtained from the E-GVAP network of ground-based GNSS receivers contain horizontally dense information and are available with a temporal resolution of up to 5 min and therefore have the potential to provide humidity-related data for kilometre-scale short-range weather forecasting. To stimulate further enhancements in the preprocessing and use of GNSS ZTD observations in NWP and nowcasting applications, in particular when forecasting severe weather, a European COST Action (ES1206) is ongoing between 2013 and 2017. One outcome of the action was a recent review of the current state of the art and future prospects of the ground-based GNSS meteorology in Europe . The action resulted furthermore in revitalization of the Nordic GNSS Analysis Centre (NGAA), now located at Lantmäteriet in Sweden, where GNSS data are processed for a large number of receiver stations, mainly from the Nordic countries. The dense network of GNSS ZTD observations provide an attractive source of supplementary humidity information to the MetCoOp modelling system.

Like all other types of measurements, the GNSS ZTD observations are associated with errors that need to be properly characterized. It has earlier been demonstrated that adaption of variational bias correction to be used together with GNSS ZTD data was successful for handling systematic observation errors . The sources of bias in the ZTD observation data with respect to the ZTD model data may be for several reasons, such as GNSS data-processing algorithms (use of mapping functions, formulation of hydrostatic delay, errors in the conversion of ray path to zenith delay) and systematic errors in both the model fields and the ZTD observation operator. In particular, a low model top will generally result in a systematically too-low model equivalent of the GNSS ZTD observations. In Sánchez-Arriola et al. only one predictor was used in the variational bias correction. Earlier, studied the behaviour of a non-adaptive multilinear bias correction scheme inspired by the one proposed by and found a benefit in using more than one predictor. The question is whether an adoptive bias correction scheme like the one used by Sánchez-Arriola et al. would also benefit from using more predictors.

Due to the measurement and processing techniques, GNSS observations are very likely to have correlated errors. The difficulties of spatially and temporally correlated observation errors have generally been circumvented in data assimilation by applying thinning of data, or through observation processing algorithms that are assumed to remove the observation error correlations . Methods have been developed to account for serially correlated errors , but there is certainly room for improvement regarding spatially correlated errors, although some general research within this area has been carried out . Some studies have focused on GNSS ZTD observations , but the handling of the correlated observation errors is still an active area of research.

Example of validation of ZTD (upper panel) and integrated water vapour (lower panel) from a week in March 2017 for station Onsala, Sweden. Statistics for different processing centres compared to a NWP run are shown in the table. The data are taken from the E-GVAP website.

[Figure omitted. See PDF]

GNSS ZTD observations processed by the NGAA centre have been used within the MetCoOp forecasting system, aiming at improving short-range forecasts of, in particular, moisture, clouds and precipitation. Two different GNSS ZTD processing techniques applied at NGAA are described, compared and evaluated. The sensitivity of the results to various aspects of the GNSS ZTD observation handling and data assimilation was investigated. The evaluation includes both an individual case study and statistics based on extended parallel experiments.

The paper is organized as follows. The GNSS data processing is the topic of Sect. . In Sect.  the NWP modelling system is described. Section  deals with the design of parallel data assimilation experiments, and their corresponding results are presented in Sect. . Finally, conclusions are presented in Sect.  together with some future plans.

GNSS data processing

In June 2016 Lantmäteriet (Swedish Mapping, Cadastre and Land Registration Authority) became NGAA, one of the analysis centres in E-GVAP, and it is in charge of the data processing for the GNSS stations in Sweden, Finland, Norway, Denmark and some IGS stations in order to provide near-real-time (NRT) ZTDs. For E-GVAP data processing the NRT product means that the estimated ZTD for the previous hour needs to be ready within 45 min. NGAA includes in total approximately 700 stations and currently provides two NRT ZTD products (NGA1 and NGA2). The NGA1 product is obtained from the Bernese v5.2 network solution, while NGA2 is given by the GIPSY/OASIS II v.6.2 data processing using the precise point positioning (PPP) strategy .

In a network solution there is no need for the precise clock product for the GNSS satellites due to the differential observables. However, the computing time will be exponentially increased as the number of GNSS stations in the data processing increases while the station-related errors are correlated to each other. During PPP processing, the data from only one GNSS station is processed at a time, meaning that station-related errors are independent from others. However, high-quality satellite clock product is critical for the accuracy of PPP data processing. More details about the two types of data processing can be found in Sect.  and .

Post-data processing

In order to obtain the best accuracy for the estimated hourly ZTD, the coordinates of the stations need to be fixed in NRT data processing. The fixed coordinates are provided by post-data processing, which is carried out once per day. However, due to the latent time of the final orbit products, the data used for the processing were actually obtained 14 days before the day being examined. The estimated coordinates will be averaged together with the coordinates estimated for the previous 6 days. The weekly averaged coordinates will be used as the fixed coordinates for hourly NRT data processing. Although for each station the fixed coordinates are the ones estimated 14 days prior, the maximum difference in the height component is less than 1 mm if no significant movements happened at the station, e.g. an earthquake, in the previous 14 days. Such a small difference will only have an insignificant impact (smaller than 0.3 mm) on the estimated ZTD.

In the post-data processing the acquired GPS phase-delay measurements are used to form ionospheric free linear combinations (LC) that are analysed by Bernese v5.2 using a network solution to estimate station coordinates together with tropospheric parameters. We used the final GPS orbit products provided by CODE (ftp.unibe.ch) and included an ocean tide loading correction using the FES2004 model . The absolute calibration of the phase centre variations (PCVs) for all antennas (IGS14.atx) was implemented . The tropospheric estimates are updated every 2 h, while 1 h estimates are given for the station coordinates. A 10${}^{\circ }$ elevation cut-off angle is used, and the slant delays are mapped to the zenith using the Vienna Mapping Function 1 (VMF1; ).

Time series of ZTD for June 2016 for the NGA1 (blue) and NGA2 (red) solutions together with the ones obtained from post-processing (green circles). The data are from the Ballerup (Copenhagen) station in Denmark. The $x$ axis shows the days in June, whereas the $y$ axis shows the ZTD in mm.

[Figure omitted. See PDF]

NGA1 data set

The NGA1 product is obtained from a Bernese hourly data processing running in near-real time and using the fixed station coordinates. We use the ultra-rapid GPS orbit products provided by CODE (ftp.unibe.ch). The ocean tide loading correction (FES2004) and the antenna PCV absolute calibration are implemented. The tropospheric estimates are updated every 15 min and a 10${}^{\circ }$ elevation cut-off angle is used with a global mapping function (GMF; ). The NGA1 product is currently under the operational status with a time delay of 45 min.

NGA2 data set

The NGA2 product is obtained from GIPSY NRT data processing where the GPS data were analysed by GIPSY-OASIS v6.2 using the PPP strategy with the fixed station coordinates. Currently we use the ultra-rapid GPS orbit and clock products provided by JPL ftp://sideshow.jpl.nasa.gov/pub/JPL_GPS_Products/Ultra. The same set-ups are used for the GIPSY data processing, i.e. FES2004 model, antenna PCV absolute calibration, a 10${}^{\circ }$ elevation cut-off angle and a GMF. The tropospheric estimates are updated every 5 min. In addition the single receiver phase ambiguity resolution is also implemented . The NGA2 product is now under a test mode due to a longer time delay of about 1.5 h for fetching the JPL ultra-rapid orbit and clock products.

Comparing the data sets

Due to the long time delay in the NGA2 data, the NGA1 data set is sent to E-GVAP for redistribution between member countries. At E-GVAP it is still in test mode, but this will change to operational in the near future. At the E-GVAP website all stations are validated against an NWP model run carried out at the Royal Netherlands Meteorological Institute (KNMI). This comparison against the model should not be taken as a validation of truth, but it makes it possible to compare the results from different processing centres calculating ZTD for the same stations. A number of stations around Europe have been selected for comparison to be supersites, which are processed by all centres that are part of E-GVAP. In Fig.  an example of such a comparison, taken from http://egvap.dmi.dk/, is shown for the station Onsala in southern Sweden. It can be seen that NGA1 compares well with most of the other centres.

In Fig.  the two solutions from NGAA are also assessed with respect to the ZTDs estimated by post-processing using the IGS final satellite orbits and clock products. The data are from the receiver in Ballerup (BUDP) just outside Copenhagen. The figure shows that the Bernese network solution (NGA1) and the GIPSY PPP solution (NGA2) have similar results with mean differences of $-$0.5 and $-$0.1 mm, respectively, while the corresponding standard deviations are 4.2 and 4.9 mm. The increased ZTD values around 22–25 June (see Fig. ) when a major convective storm passed over south-western Denmark and southern Sweden are interesting to note. This will be discussed further in the case study in Sect. .

Time series of ZTD for 22 to 24 June 2016 for the NGA1 (blue) and NGA2 (red) solutions together with the ones obtained from a post-processing (green circles). The data are from the Ballerup (Copenhagen) station in Denmark. The $x$ axis shows the days in June, while the $y$ axis shows the ZTD in mm.

[Figure omitted. See PDF]

The NWP modelling system

The main components of the MetCoOp ALADIN-HIRLAM NWP system are surface data assimilation, upper-air data assimilation and the forecast model.

The forecast model configuration, e.g. dynamical core and physical parameterizations, are described in detail in and . It has a spectral representation with a non-hydrostatic formulation. Stratiform and deep convective clouds are explicitly represented, while for shallow convection a sub-grid parameterization is applied using the EDMF (eddy-diffusivity mass-flux) scheme. The representation of the turbulence in the planetary boundary layer is based on a prognostic turbulent kinetic energy (TKE) equation combined with a diagnostic mixing length . The radiative transfer of the shortwave spectrum is described with six spectral bands , and the long-wave radiation is modelled in accordance with . Surface processes are modelled using SURFEX . Lateral boundary conditions were provided by 6-hourly European Centre for Medium-Range Weather Forecasts (ECMWF) operational forecasts with a 1 h time resolution. In addition, a spectral large-scale mixing of the background state, the 3 h HARMONIE forecast, fields with the lateral boundary ECMWF fields was applied. In this way we hoped to benefit from the high-quality large-scale information from the ECMWF global forecasts in the regional MetCoOp data assimilation.

In the MetCoOp setup there are 750 $\times$ 960 horizontal grid-points at each of the 65 vertical levels extending up to 10 hPa, which approximately corresponds to a height of 32 km in the atmosphere. The horizontal grid distance is 2.5 km. The model domain is illustrated by the black frames in Fig. . In the surface data assimilation SYNOP observations of temperature and relative humidity at the vertical level of 2 m were utilized. In addition sea surface temperature and sea ice concentration from an oceanographic model were used. In the MetCoOp upper-air data assimilation conventional types of in situ measurements were used and these include radiosonde, pilot-balloon wind, SYNOP, ship and aircraft measurements. In addition radiances from the AMSU-A, AMSU-B/MHS and IASI instruments onboard polar-orbiting satellites, as well as surface winds from the Advanced Scatterometer (ASCAT) instrument, were used. Furthermore, humidity observations from networks of ground-based weather radars and GNSS receiver stations were used. The radar reflectivity is not directly assimilated into the model since there is a complicated, non-linear relationship between the model variables and reflectivity. This includes parameterizations of microphysical processes and non-Gaussian error distributions. Instead a vertical moisture profile is retrieved through a one-dimensional (1-D) Bayesian retrieval based on a comparison between observed and simulated reflectivities . Observations used were obtained from the Global Telecommunications System (GTS), the EUMETSAT Advanced Retransmission Service (EARS), the advanced weather radar network for the Baltic Sea Region (BALTRAD) data hub and the E-GVAP retransmission service.

The surface data assimilation uses an optimal interpolation scheme . In the current study a 3-D variational upper-air data assimilation (3D-Var) scheme was applied within a 3 h data assimilation cycle. The climatological background error statistics used in the current study were derived from an ensemble of MetCoOp forecast differences obtained through downscaling of the forecast fields based on the ECMWF ensemble data assimilation (EDA). The ECMWF EDA-based forecast fields were horizontally and vertically interpolated to the HARMONIE AROME 2.5 configuration geometry and used as initial conditions for high-resolution non-hydrostatic model runs. The ECMWF EDA uses T399 horizontal resolution and 91 vertical levels. Then the evolved high-resolution ensemble was scaled to be consistent with the amplitude of the 3 h forecast error for HARMONIE-AROME. The values applied correspond roughly to a GNSS ZTD background error standard deviation. Recently ECMWF have increased the horizontal resolution of the EDA system to T639 and demonstrated clear improvements from this change of resolution . One could expect that re-derivation of the MetCoOp forecast differences utilizing the enhanced ECMWF EDA system would lead to improved MetCoOp background error statistics and thus an improved data assimilation system. We have therefore re-calculated the background error statistics utilizing the enhanced ECMWF EDA system and carried out sensitivity experiments with the new background error statistics. Results are presented in Sect.  as an example of how GNSS ZTD data assimilation can gain from more general data assimilation improvements. The background error statistics are specified for assimilation control variables. These are vorticity, unbalanced divergence, unbalanced temperature, unbalanced surface pressure and unbalanced specific humidity . Important upper-air data assimilation observation-handling components are the modelling of observation counterparts with an observation operator, quality control mechanism, thinning, bias correction and error specification. The observation operator projects the model state onto the GNSS ZTD observation. Since a variational framework is used, non-linear and the corresponding tangent-linear and adjoint versions of the observation operator are needed. The ZTD observation operator ($H$), given a station location (including altitude), calculates the model equivalent of the GNSS ZTD by integrating the model-calculated refractivity vertically from the station height to the model top, as described in . In the MetCoOp system, following the ideas of , we have extended the observation operator with the possibility of accounting for the contribution to the ZTD by the part of the atmosphere above the model top. Details of the observation handling within the data assimilation with emphasis on GNSS ZTD are given in .

The GNSS ZTD observation errors of the observations accepted for the data assimilation were assumed to have a Gaussian error distribution with an observation error standard deviation of 12 mm. This observation error standard deviation was derived from observation minus background and observation minus analysis departures, and it was empirically adjusted so that roughly the same weight was given to the observation and to the background. Objective methods such as the one proposed by could in future be tried instead to tune the observation error variance.

There is also an additional quality control mechanism within the assimilation. The purpose of this quality control mechanism is to remove observations affected by gross errors and a central part is the background check. The observation, $y_i$, is rejected if it does not satisfy the following inequality:

$${\left(\left[,H,\right(,{\mathit{x}}^b,),{]}_i,-,y_i\right)}^{\mathrm{2}}/{\mathit{\sigma }}_{b,i}^{\mathrm{2}}>L\times \mathit{\lambda },$$ where $\mathit{\lambda }=\mathrm{1}+{\mathit{\sigma }}_{o,i}^{\mathrm{2}}/{\mathit{\sigma }}_{b,i}^{\mathrm{2}}$, $L$ is the rejection limit and $\left[H\right({\mathit{x}}^b){]}_i$ denotes the projection of the model state on observation $i$. In the background guess check, the background and observation error standard deviation were set to 10 and 12 mm, consistent with the values used in 3D-Var. The rejection limit for GNSS ZTD observations in the HARMONIE system was set to 4. This value resulted in a relatively strict background quality control mechanism of GNSS ZTD observations.

To alleviate the effects on the initial state of spatially correlated observation errors caused by, for example, orographic effects, we applied spatial thinning to GNSS ZTD observations. The default thinning distance was on the order of 80–100 km. The thinning distance was used when selection receiver stations so that receiver stations closer than 80–100 km to each other were not used. This thinning distance was also rather empirically determined. The thinning is applied in the form of a static whitelist based on studies of data availability and statistics of observation minus background equivalent statistics from a spin-up period. Thus the thinning is static so that each data assimilation cycle of observations from the same set of GNSS ZTD receiver stations is used. A study of the sensitivity of reducing the thinning distance can be found in Sect. . A next step would be to apply objective methods such as the one proposed by ; instead of tuning the thinning distance, the observation error covariance could be modelled.

Systematic errors in the GNSS ZTD data were handled by an adaptive variational bias correction (VarBC). Within VarBC the bias was represented by coefficients for the selected predictors. These predictors were estimated within the variational data assimilation process simultaneously as deriving the assimilation control vector for the model state while minimizing the cost function . The bias correction was carried out individually for each receiver station and in the default version only one predictor, in the form of an offset value, was used. However, there is also the possibility of introducing more predictors, like 1000–300 hPa thickness and TCWV. The sensitivity to introducing extra predictors is investigated in Sect. .

MetCoOp model domain (black frame) and GNSS ZTD observation usage for operational MetCoOp (a), NGAA usage with 80–100 km thinning distance (b) and NGAA usage with 40 km thinning distance (c).

[Figure omitted. See PDF]

Experimental design

In order to investigate the potential benefit in the MetCoOp system of utilizing NGAA GNSS ZTD, a number of parallel data assimilation and forecast experiments have been carried out. Furthermore the parallel experiments were aimed at investigating the sensitivity of the GNSS ZTD data assimilation to various aspects of the data assimilation. A copy of the MetCoOp operational configuration was run with a 3 h data assimilation cycle for the period 1–30 June 2016 and with a 1-month spin-up period before that. This particular month was chosen because it was characterized by several heavy precipitation events. We expect that additional moisture-related observations should be particularly beneficial for prediction of such weather conditions. We ran short-range forecasts every 3 h to provide the background for the next analysis, and we launched forecasts of up to 36 h, four times per day, from 00:00, 06:00, 12:00 and 18:00 UTC. In total there were four data assimilation studies (A–D), each involving two or more parallel data assimilation experiments. These parallel experiments are explained as follows:

• A 

  Assessing the impact of assimilating GNSS ZTD from the NGAA processing centre.

• Observation usage as in MetCoOp operational, including GNSS ZTD from ROBH and METO processing centres.

• Observation usage as in MetCoOp operational, except that GNSS ZTD observation usage was extended to also include receiver stations from the NGAA processing centre, processed with the Bernese approach.

• Observation usage as in MetCoOp operational, except that GNSS ZTD observation usage was extended to also include receiver stations from the NGAA processing centre, processed with the GIPSY approach.

• B 

  Assessing the impact of different VarBC predictor choices.

• Observation usage as in A2 above, i.e. utilizing one predictor in the form of an offset value for the GNSS ZTD variational bias correction.

• Observation usage as in A2 above, except that the variational bias correction was extended to two predictors: offset value and 1000–300 hPa thickness.

• Observation usage as in A2 above, except that the variational bias correction was extended to two predictors: offset value and TCWV.

• C 

  Assessing the impact of modifying thinning distances for GNSS ZTD.

• Observation usage as in A2, i.e. utilizing one predictor in the form of an offset value for the GNSS ZTD variational bias correction and a GNSS ZTD thinning distance on the order of 100 km.

• Observation usage as in A2, except that a GNSS ZTD thinning distance on the order of 40 km was used.

• D 

  Assessing the potential benefit of general data assimilation improvements on GNSS ZTD utilization for NWP.

• Observation usage as in A2, i.e. utilizing one predictor in the form of an offset value for the GNSS ZTD variational bias correction and a GNSS ZTD thinning distance on the order of 100 km.

• Observation usage as in A2, except that an improved $\mathbf{B}$ matrix was used.

The MetCoOp model domain and the GNSS ZTD observation usage in the operational set-up and in the NGAA ZTD observation usage when applying different thinning distances are illustrated in Fig. , where the left panel shows experiment A1, the middle panel A2 (and A3, B1, B2, B3, C1, D1, D2) and the right panel shows experiment C2. Note that in all experiments of studies A, C, and D only one predictor in the form of an offset value was used. In all experiments in studies B and D a 100 km thinning distance was used. All experiments in studies A, B and C used the operationally used $\mathbf{B}$ matrix and all experiments in studies B, C and D used the NGA1 data set, processed with the Bernese approach.

Results

Verification methods

To evaluate the relative quality of the analyses and subsequent forecasts from the different parallel experiments, we verified them against radiosonde and SYNOP observations within the model domain. The verification was carried out for weather parameters at the surface level and for the upper-air parameters, wind, temperature and humidity. The model data used in the statistics were the analyses and forecasts of up to 24 h. Special emphasis was put on verification of humidity and precipitation. In addition we used the degrees of freedom for signal (DFS) to study the relative impact of observations in the assimilation system . DFS is the derivative of the analysis increments in observation space with respect to the observations used in the analysis system. As proposed by , DFS can be computed through a randomization technique:

Degree of freedom of signal subdivided into various observation types for the four experiments A1, A2, C2 and D2. Results were based on data from eight different data assimilation cycles.

[Figure omitted. See PDF]

Bias and standard deviation of $+$12 h relative humidity (unit: %) forecasts as a function of vertical level for verification against radiosonde observations. Scores are for experiments A1 (red), A2 (blue) and A3 (green).

[Figure omitted. See PDF]

One-month time series of GNSS ZTD (unit: m) from the Onsala receiver station. Analysed (black), background (blue), observation after bias correction (green) and observation before bias correction (red).

[Figure omitted. See PDF]

Bias and standard deviation of $+$12 and $+$24 h relative humidity (unit: %) forecasts as a function of vertical level for verification against radiosonde observations. Scores are for experiments B1 (red), B2 (blue) and B3 (green).

[Figure omitted. See PDF]

$$\begin{array}{cc}{\displaystyle }\mathrm{DFS}& {\displaystyle }={\displaystyle \frac{\partial H{\mathbf{x}}^{\mathbf{b}}}{\partial y}}\approx (\stackrel{\mathrm{̃}}{\mathit{y}}-\mathit{y}){\mathbf{R}}^{-\mathrm{1}}(\mathbf{H}{\stackrel{\mathrm{̃}}{\mathit{x}}}^a-{\mathbf{Hx}}^b)\\ & {\displaystyle }& {\displaystyle }-({\mathbf{Hx}}^a-{\mathbf{Hx}}^b),\end{array}$$ where $\mathit{y}$ is the vector of the observations, $\stackrel{\mathrm{̃}}{\mathit{y}}$ is the vector of perturbed observations, $\mathbf{R}$ is the observation–error covariance matrix, $\mathbf{H}$ is the tangent-linear observation operator for each observation type, ${\mathit{x}}^a$ and ${\mathit{x}}^b$ are respectively the analysis and the background state and ${\stackrel{\mathrm{̃}}{\mathit{x}}}^a$ is the analysis produced with perturbed observations. The previous formulation can be applied to any subset of observations . The absolute DFS represent the information brought into the analyses by the different observation types, in terms of amount, distribution, instrumental accuracy and observation operator definition. They offer an insight into the actual weight given to the observations within the analysis system in terms of self-sensitivity of the observations (i.e. sensitivity at location of observation). However, they do not provide any information on the spatial or cross-correlations between the observations and the analysis. There is also the possibility of estimating the DFS per observation through calculation of relative DFS, by normalizing the absolute DFS by the amount of the observations belonging to a specific subset. Here we have, however, chosen to focus on absolute DFS.

Bias and standard deviation of $+$12 and $+$24 h relative humidity (unit: %) forecasts as a function of vertical level for verification against radiosonde observations. Scores are for experiments C1 (red) and C2 (blue).

[Figure omitted. See PDF]

Kuiper skill score for 12 h accumulated precipitation (left) and $+$6 and $+$18 h cloud forecasts (right) for verification against SYNOP stations in the domain.

[Figure omitted. See PDF]

The different kinds of objective statistical verifications described above were also complemented with a more subjective verification for an individual case study.

Impact on analyses

For the DFS computation, eight cases 4 days apart were chosen: 00:00 UTC (2 June), 03:00 UTC (6 June), 06:00 UTC (10 June), 09:00 UTC (14 June), 12:00 UTC (18 June), 15:00 UTC (22 June), 18:00 UTC (26 June) and 21:00 UTC (30 June). The cases were chosen to reduce the interdependency between the initial conditions, and to obtain data from data assimilation cycles covering different times of the day.

In Fig.  the DFS calculated separately for different observation types and parameters are shown. The values represent the sum over the observations belonging to the same subset of Eq. () calculated for each individual observation. Results are shown for the four experiments A1, A2, C2 and D2. The rest of the experiments all have DFS similar to A2 and are therefore not shown. Comparing the DFS of A1, A2 and C2 showed that the contribution from GNSS ZTD increased with an increasing number of GNSS ZTD observations. A clear interaction with moisture-related observations from IASI and radar can also be seen. The larger DFS of GNSS ZTD after increasing the number of GNSS observations was associated with an increase in DFS from radar-based humidities and a decrease in DFS from the IASI instrument, providing satellite-based humidity information. It is also evident by comparing A2 and D2 that by improving the background error statistics the DFS for GNSS ZTD and also of other observations can be increased. For the DFS scores not shown, the difference in impact on analysis from NGA1 and NGA2 was very small, and the impact on DFS for GNSS ZTD of introducing more predictors in the variational bias correction of GNSS ZTD was also very limited.

Bias and standard deviation of temperature (unit: K) and relative humidity (unit: %) as a function of forecast range. Scores for temperature (a, c) and relative humidity (b, d) at the vertical levels of 300 hPa (a, b) and 850 hPa (c, d). Scores are for experiments D1 (red) and D2 (blue).

[Figure omitted. See PDF]

Statistical verification of forecasts

In Fig.  the scores for verification of $+$12 and $+$24 h relative humidity forecasts of the experiments A1–A3 against radiosonde observations within the domain are shown for different vertical levels. A small positive impact on forecasts can be seen from utilizing NGAA ZTD observations. The positive impact was slightly more pronounced when the NGAA observations were in the form of the NGA1 data set. For variables other than humidity, the impact on the forecast quality was small (not shown).

The impact of utilizing two predictors in the variational bias correction of GNSS ZTD is small, not only in terms of DFS. As another example, Fig.  shows, for one particular receiver station (Onsala), a 1-month time series during the experiment GNSS ZTD of background state equivalent (FG), analysis (AN), observed value before bias correction (OBS RAW) and observed value after bias correction (OBS) for the three different experiments B1–B3. It can be seen that the bias correction worked properly, managing to correct for the systematic difference between the raw observation and the model state equivalents. On the other hand, it was evident that the time evolution of the bias-corrected observations was very similar between the three different runs. The difference between introducing the second predictor in the form of 1000–300 hPa thickness or TCWV was very small.

The small impact of introducing additional predictors in the adaptive bias correction was confirmed also by forecast verification scores. Figure  illustrates this impact on $+$12 and $+$24 h relative humidity forecasts, for verification against radiosonde observations. As for forecasts of other variables (not shown), the impact on relative humidity forecast quality of introducing more predictors was small.

Radar-based accumulated precipitation (a, b, unit mm(3 h)${}^{-\mathrm{1}}$) between 24 June on 23:00 UTC and 25 June 2016 on 02:00 UTC (a), and between 02:00 and 05:00 UTC on 25 June 2016 (b). Forecasted accumulated precipitation (unit: mm(3 h)${}^{-\mathrm{1}}$) based on C1 (c, d) and C2 (e, f) between 00:00 and 03:00 UTC, and between 03:00 and 06:00 UTC on 25 June 2016. The forecasts were initiated from 24 June 2016 at 12:00 UTC.

[Figure omitted. See PDF]

Rain-gauge-observed 3 h accumulated precipitation (unit mm(3 h)${}^{-\mathrm{1}}$) from accumulation period starting on 24 June 2016 at 23:00 UTC (a) and on 25 June 2016 at 02:00 UTC (b).

[Figure omitted. See PDF]

The sensitivity of modifying the thinning distance applied to GNSS ZTD observations is illustrated in Fig. . From the left part of the figure it can be seen that in terms of standard deviation the impact was rather small, except for improved humidity forecasts at the lowest levels when reducing the thinning distance from 80–100 to 40 km. The right part of the figure shows that this improvement was present at forecast ranges up to 36 h. In terms of bias, on the other hand, it can be seen from the left panel that there was an increased positive humidity bias throughout the lower troposphere with reduction of the thinning distance. Again, for forecasts of other variables the impact was small (not shown). An increased humidity bias when reducing the thinning distance was also noticed by . It was speculated whether, in the lack of high-resolution complementary unbiased humidity information, nearby GNSS ZTD receiver stations affect each other during the spin-up phase of predictor coefficients. They only used conventional types of observations in addition to the GNSS ZTD observations. Our study confirmed that the increased bias when reducing the thinning distance was present also when a substantial amount of humidity-related remote sensing observations, such as AMSU-B/MHS, IASI and radar-derived humidities, were assimilated in addition to GNSS ZTD. It should be kept in mind, however, that none of these data sources are assumed to be bias free. For AMSU-B/MHS and IASI a variational bias correction was applied, and for radar-derived moisture information a pre-processing utilization of the model background field was applied . Our results from Sect.  hint that there was a relationship between IASI, radar and the GNSS ZTD impact when modifying the GNSS ZTD thinning distance. However, the interaction of reduction of thinning distances and increased bias needs to be better understood before one can fully benefit from reducing the GNSS ZTD thinning distance. This is one of the aims for future in-depth studies with the MetCoOp data assimilation system.

In addition to improvement of the low-level humidity forecasts when reducing the thinning distance to 40 km, a slight improvement was seen in forecasts of cloud cover and more pronounced improvements in precipitation forecasts, as illustrated in Figure , in terms of the Kuiper skill score. It should, however, be kept in mind that there are some known problems related to precipitation and cloud measurements . Thus, despite the increased bias in humidity related to the reduction in thinning distance, the improvements in terms of standard deviations for humidity forecasts resulted also in improvement in the humidity-related variables of cloud and precipitation. The question whether improvements could also be seen in an individual case is addressed in Sect. .

When investigating the improvements to the system that can be brought by adding new observations and by refinements in the observation handling, it is also useful to get an idea of how much the extraction of information from the new observations, as well as from all the other observations, can be improved by general data assimilation improvements. In our case, the general data assimilation improvements were given by an improved representation of background error statistics. The improved background error statistics had a positive impact on the forecasts, shown in Fig.  for temperature and relative humidity scores. A positive impact was also found on surface pressure forecasts and wind forecasts (not shown).

General data assimilation improvements, like the improved $\mathbf{B}$ matrix presented here, influenced more aspects and observations of the data assimilation system than just GNSS ZTD observations. It is important to keep in mind that such general improvements can also support obtaining more useful information from both newly introduced observation types as well as those that have been in use for some time.

Case study

To investigate whether the modification of thinning distance has any noticeable effect on individual weather conditions, we looked into one particular case in more detail. The individual case selected was a heavy precipitation event that took place over south-western Denmark and southern Sweden during the night/early morning between 24 and 25 June 2016. The upper row of Fig.  shows the radar-derived gauge-adjusted 3 h accumulated precipitation between 23:00 UTC on 24 June and 02:00 UTC on 25 June, and between 02:00 UTC and 05:00 UTC on 25 June. The middle and lower panels show accumulated precipitation forecasts for the runs with 80–100 and 40 km thinning distance, respectively. For this particular case both of the runs had a phase/timing error of roughly 1 h. Therefore, to reduce the effect of the phase/timing error on the verification, the accumulation interval for the precipitation of the forecasts was shifted 1 h relative to the radar-derived precipitation. For the forecasts the accumulation interval was between 00:00 and 03:00 UTC on 25 June, and between 03:00 to 06:00 UTC on 25 June. Between 00:00 and 03:00 UTC the forecasts of the two runs were rather similar, but as the system moved toward the north-east more of the intensity and structure in accordance with observations was retained in the run, with the 40 km thinning distance; however, despite the phase correction, there was a small error in position between 03:00 and 06:00 UTC. Figure  shows 3 h accumulated precipitation from rain gauges for which nonzero precipitation was registered for this particular case. Due to the phase error mentioned in the forecasts we have again chosen to show the accumulation both from 24 June 2016 at 23:00 UTC (left panel) and from 25 June 2016 02:00 UTC (right panel) that was shifted 1 h as compared with the forecasted accumulated precipitation in Fig. . By comparing Figs.  and  one can see that rain gauge observations also support the fact that the forecast of the run with 40 km thinning distance was better.

Conclusions

The processing of GNSS ZTD data from the newly vitalized NGAA processing centre has been described in detail. It is shown that these data have the capability to enhance the NWP forecasts, in particular for humidity, when introduced, in addition to other observations, in the HARMONIE-AROME model. The sensitivity of the forecasts to the two solutions provided for estimating ZTD and to the various settings in the GNSS ZTD data processing has been investigated. The two different methods of estimating ZTD generated very similar results and the impact on the forecasts was therefore also very small. It was also found that the results were rather insensitive to the number of predictors used in the variational bias control. In this study only two predictors were tested at the same time. It might be useful as a next step to test more than two and other parameters, e.g. surface pressure. In contrast to the small impact from the VarBC predictors, the results were rather sensitive to the choice of thinning distance applied. There are potential improvements from reducing the thinning distance of the ZTD observations to make use of more data, but there are also related issues. Reducing the thinning distance resulted in increased humidity forecast biases in the lower troposphere. This may have been due to increased influence from correlation errors; this needs to be investigated further to find the best trade-off between the number of observations and the influences of error correlations. In general the horizontal observation error correlations need to be investigated further, for example by applying techniques proposed by and by taking spatial error correlations into account.

The assimilation of GNSS ZTD in NWP can benefit from more general data assimilation improvements, such as enhanced description of statistical information or improved data assimilation algorithms. In this paper this was highlighted by providing an example in the form of an additional run carried out with what we think is an enhanced description of the background error statistics. Clearly the enhanced description resulted in better use of the GNSS ZTD observations in the NWP system. It is important, however, to keep in mind that such general data assimilation aspects influence not only the GNSS ZTD observation usage but also all other observations. In addition, further developments of the data assimilation algorithms, e.g. the impact on utilization of GNSS ZTD observation in a 4-D variational data assimilation, will be investigated.

Underlying research data are available upon request to [email protected].

The authors declare that they have no conflict of interest.

This article is part of the special issue “Advanced Global Navigation Satellite Systems tropospheric products for monitoring severe weather events and climate (GNSS4SWEC) (AMT/ACP/ANGEO inter-journal SI)”. It is not associated with a conference.

Acknowledgements

The authors would like to thank Jan Johansson at Chalmers Technical University for support with revitalizing the NGAA processing centre. The help from Roger Randriamampianina with computations of DFS is also greatly appreciated. The work was carried out within the framework of the European COST Action ES1206 concerned with “Advanced Global Navigation Satellite Systems tropospheric products for monitoring severe weather events and climate”. The constructive comments and suggestions made by two anonymous reviewers were greatly appreciated. Edited by: Jonathan Jones Reviewed by: two anonymous referees