ACPAtmospheric Chemistry and PhysicsACPAtmos. Chem. Phys.1680-7324Copernicus PublicationsGöttingen, Germany10.5194/acp-17-13983-2017Data assimilation of GNSS zenith total delays from a Nordic processing centreLindskogMagnusmagnus.lindskog@smhi.seRidalMartinThorsteinssonSigurdurNingTong1Swedish Meteorological and Hydrological Institute, Norrköping, Sweden2Icelandic Meteorological Office, Reykjavík, Iceland3Lantmäteriet, Gävle, SwedenMagnus Lindskog (magnus.lindskog@smhi.se)24November20171722139831399819June20171August201727September201711October2017This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/This article is available from https://www.atmos-chem-phys.net/17/13983/2017/acp-17-13983-2017.htmlThe full text article is available as a PDF file from https://www.atmos-chem-phys.net/17/13983/2017/acp-17-13983-2017.pdf

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.

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∘ 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.

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∘ 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∘ 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.

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 × 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, yi, is rejected
if it does not satisfy the following inequality:

[H(xb)]i-yi2/σb,i2>L×λ,
where λ=1+σo,i2/σb,i2, L is the rejection limit
and [H(xb)]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. .

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 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 B matrix and all
experiments in studies B, C and D used the NGA1 data set, processed with the Bernese approach.

ResultsVerification 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.

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).

DFS=∂Hxb∂y≈(ỹ-y)R-1(Hx̃a-Hxb)-(Hxa-Hxb),
where y is the vector of the observations, ỹ is the vector
of perturbed observations, R is the observation–error covariance matrix,
H is the tangent-linear observation
operator for each observation type, xa and xb are respectively the
analysis and the background state and 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).

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.

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)-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)-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.

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

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 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 magnus.lindskog@smhi.se.

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