Consistency and representativeness of integrated water vapour from ground-based GPS observations and ERA-Interim reanalysis

This study examines the consistency and representativeness differences of daily integrated water vapour (IWV) data from ERA-Interim reanalysis and GPS observations at 120 global sites over a 16-year period (1995–2010). Various comparison statistics are analysed as a function of geographic, topographic, and climatic features. A small (±1 kg m⁻²) bias is found in the reanalysis across latitudes (moist in northern and southern midlatitudes and dry in the tropics). The standard deviation of daily IWV differences is generally below 2 kg m⁻² but peaks in the northern and southern storm-track regions. In general, the larger IWV differences are explained by increased representativeness errors, when GPS observations capture some small-scale variability that is not resolved by the reanalysis. A representativeness error statistic is proposed which measures the spatiotemporal variability in the vicinity of the GPS sites, based on reanalysis data at the four surrounding grid points. It allows to predict the standard deviation of daily IWV differences with a correlation of 0.73. In general, representativeness differences can be reduced by temporal averaging and spatial interpolation from the four surrounding grid points. A small number of outlying cases (15 sites) which do not follow the general tendencies are further examined. It is found that their special topographic and climatic features strongly enhance the representativeness errors (e.g. steep topography, coastlines, and strong seasonal cycle in monsoon regions). Discarding these sites significantly improves the global ERA-Interim and GPS comparison results. The selection of sites a priori, based on the representativeness error statistic, is able to detect 11 out of the 15 sites and improve the comparison results by 20 % to 30 %.

1 Introduction

Quantifying the global atmospheric moisture distribution and its variability across timescales remains a challenge to the climate community. Atmospheric reanalyses offer a comprehensive representation of the various components of the hydrological cycle, among which precipitation and evaporation are the dominant terms at the larger space scales and timescales. However, both quantities result from model integrations and are not strongly constrained by observations (Trenberth et al., 2011). The difference of precipitation minus evaporation corresponds to the net vertically integrated atmospheric moisture convergence, a quantity which can also be computed from analysed three-dimensional moisture and wind fields which benefit directly from the assimilation of observations (Trenberth and Fasullo, 2013). However, due to the high spatiotemporal variability of atmospheric moisture, the quality of moisture fields in the reanalyses remains limited, especially in data-sparse areas (Trenberth et al., 2005; Meynadier et al., 2010).

Ground-based Global Positioning System (GPS) integrated water vapour (IWV) observations have been used for some time as an independent validation source for global atmospheric reanalyses over limited regions and periods (Hagemann et al., 2003; Bock et al., 2005, 2016; Heise et al., 2009; Bock and Nuret, 2009) and moist atmospheric process studies (Bastin et al., 2007; Bock et al., 2008; Koulali Idrissi et al., 2012; Means, 2013; Adler et al., 2016; Khodayar et al., 2018). More recently, the value of continuous long time series of GPS IWV data has been investigated for the purpose of studying global and regional climate variability and validating climate models (Nilsson and Elgered, 2008; Vey et al., 2009; Roman et al., 2012; Ning et al., 2013; Chen and Liu, 2016; Wang et al., 2016; Parracho, 2017; Bastin et al., 2019). These studies reported various levels of agreement between GPS and atmospheric models/reanalyses making it difficult to draw general conclusions on the consistency between products. Indeed, the results depend on the model horizontal and vertical resolution, the method employed for the correction of vertical displacement between the model grid points and stations, and the considered geographical area and period of time. The influence of the model horizontal resolution suggests that representativeness differences exist between the model gridded data and station point observations. Such a situation is commonly faced in data assimilation when the station observations capture small-scale variability that is not resolved by the numerical model (Lorenc, 1986; Janjić and Cohn, 2006; Waller et al., 2014). In this context, it is traditional to include the representativeness error into the observation error in addition to the instrument error. For observations of highly variable fields such as humidity, representativeness errors can be considerably larger than instrument error and are state dependent and correlated in time (Janjić and Cohn, 2006). A proper treatment of representativeness errors, especially for humidity observations, is thus expected to improve the assimilation scheme (Waller et al., 2014). To our knowledge, representativeness errors of IWV observations, either from ground-based GPS or satellites, have not been discussed in this context. Representativeness errors arise also when measurements from different instruments are compared. This situation has been discussed for IWV measurements, e.g. by Liou et al., 2001, and Buehler et al., 2012. The representativeness errors represent in this case the effect of measurements not being perfectly co-located in space and time and using different sampling/measurement characteristics (e.g. point measurement vs. average over area/volume, instantaneous vs. time average; Buehler et al., 2012).

In the present study, we seek to analyse differences in daily IWV from ground-based GPS observations and the ECMWF global reanalysis, ERA-Interim (Dee et al., 2011), and identifying the proportion due to representativeness errors. In this context, we consider the GPS IWV observations as the reference and attribute the source of the representativeness errors to the coarse spatial resolution of the reanalysis. This choice is arbitrary and the results could be interpreted the other way round. Since GPS observations and model fields do not represent exactly the same quantity, representativeness errors can also be understood as representativeness differences more generally. Representativeness differences set a limit on the best achievable agreement between global reanalyses and station observations.

Among the motivations of this work, one is to explain the large systematic differences (biases) between GPS and atmospheric models often observed in coastal and mountainous regions (Hagemann et al., 2003; Bock et al., 2005; Parracho et al., 2018). In coastal areas, model grid cells can contain a fraction of IWV over sea not consistent with the GPS observations over land. In mountains, the model IWV can be strongly biased compared to GPS observations made in valleys or uphill. Biases amount typically to −40 % IWV per kilometre of height difference (Bock et al., 2005). Since model values are computed above a smooth orography, which can strongly depart locally from the real topography, a vertical correction must be applied, especially for IWV because the water vapour mixing ratio is the largest in the atmospheric boundary layer. Variation of biases/differences between GPS and models is also observed as a function of latitude and season (Roman et al., 2012; Ning et al., 2013; Parracho et al., 2018). Absolute IWV differences have a tendency to be larger in moister and warmer regions/periods, while relative differences tend to be larger in colder and drier regions/periods, globally. The reasons for these spatial and temporal variations are not clearly understood yet, though it can be guessed that the IWV fields in the model/reanalysis do not have the same quality in all regions of the globe and all periods of the year. Similarly, the GPS measurement and/or processing errors can change over space and time (e.g. during disturbed/severe meteorological events the mapping function errors would be larger; Boehm et al., 2007).

The primary goal of this study is thus to analyse the consistency global of daily IWV data from the ERA-Interim reanalysis and GPS station observations, and explain the contribution of representativeness errors/differences. To this purpose, we use simple statistics (mean differences and standard deviations, such as those found in most past studies) to quantify the differences between both datasets. We investigate the dependence of these statistics upon latitude, altitude, and time, as well as mean atmospheric moisture content and its spatial and temporal variability. A representativeness error statistic is introduced which quantifies the spatial variability in the ERA-Interim data at the surrounding grid points and explains to a good degree the observed differences between the reanalysis and the observations. All the statistics are computed over a period of 16 years because we want to characterize the systematic ERA-Interim minus GPS differences and not their changes over time (e.g. due to inhomogeneity and/or changes in the quality in either of the datasets). The changes over time are small in magnitude (Parracho et al., 2018) and have negligible impact on the average statistics computed here. After establishing the contribution of representativeness errors, we address the following specific questions:

1. By which means is it possible to mitigate the representativeness errors?

2. Does horizontal interpolation of model values degrade or increase the representativeness error?

3. Can we separate outlying results (e.g. sites with extreme biases and dispersion) due to enhanced representativeness errors those due to enhanced GPS instrument errors? To tackle this question, the seasonal variation of the comparison statistics and of the atmospheric environment (mean IWV and variability) is also analysed.

4. How efficient is the representativeness error statistic in detecting these outlying sites?

The results from this study are important for various applications where IWV data from reanalyses and GPS observations are used jointly. For example, recent attempts have been made to use the ERA-Interim reanalysis as a reference for detecting breaks in the GPS time series (Vey et al., 2009; Ning et al., 2016; Van Malderen, 2017). Outlying sites should be inspected more carefully to determine if the causes for the discrepancy are rather with GPS instrument errors or with reanalysis representativeness errors. This study may also contribute to a better treatment of ground-based Global Navigation Satellite System (GNSS) observation error in data assimilation (in this case, interpreting the representativeness error as an observation error), e.g. by establishing a parametric model of observation error depending on the spatiotemporal variability of IWV around the GNSS site computed from the model fields.

The paper is organized as follows. Section 2 describes how the IWV data from the two datasets are prepared. Special effort is made to use a procedure that maximizes the consistency between the datasets. Section 3 presents the results of IWV difference statistics and analyses their dependence upon a variety of parameters. General tendencies are derived that describe the consistency between the reanalysis and GPS globally. Section 4 introduces a range check which detects 15 outlying sites for which the IWV differences are especially large. The geographic, topographic, and seasonal characteristics of these sites are analysed and site-specific representativeness errors are highlighted. Section 5 discusses the possibility for detecting outlying sites a priori and concludes the paper.

2 Data and methods

2.1 GPS

In this study, we use the tropospheric delay estimates from the first reprocessing of the International GNSS Service (IGS), referred to as IGS repro1 (Byun and Bar-Server, 2009; IGSMAIL-6298, 2010). It includes results for 456 stations over the period from January 1995 to December 2010. Because we are interested in characterizing the systematic differences between GPS and atmospheric reanalyses, a subset of 120 stations which have the longest time series (16 years) is extracted. The zenith tropospheric delay (ZTD) estimates, which are available with a time sampling of 5 min, are first screened for outliers as described in Parracho et al. (2018) and averaged in hourly bins centred on the round hours (00:00, 01:00 UTC, etc.). Next, the hourly ZTDs are converted to IWV using 6-hourly surface pressure, P_s, and weighted mean temperature, T_m, computed from ERA-Interim pressure-level data (see Appendix B in Parracho et al., 2018). No temporal interpolation is applied here so that only the 1-hourly ZTD estimates matching the times of the reanalysis (00:00, 06:00 UTC, etc.) are converted. Finally, the daily IWV values are computed from five 6-hourly values between 00:00 UTC of the current day and 00:00 UTC of the next day, with weights 1∕8, 1∕4, 1∕4, 1∕4, and 1∕8, respectively. Monthly averages are computed directly from the 6-hourly values within the given month to the condition that at least 60 values are available (similar to Parracho et al., 2018). As already mentioned above, inhomogeneities in the GPS IWV time series due to equipment changes are not corrected here. This does not impact the conclusions since we analyse only overall statistics (means and standard deviations) computed over 16 years but not linear trends, and the data have been screened beforehand. Figure 1 shows the stations used in this study. The GPS coordinates, the altitudes of the reanalysis grid points in the vicinity of the GPS stations, and the number of daily and monthly values for each station are given in Table S1 in the Supplement.

Figure 1 Map showing the 120 GPS stations used in this study. A dynamic map including geographical and technical information for all the GPS sites can be found on http://www.igs.org/network (last access: 22 July 2019). Outlying sites (named in red) are detected using a range check based on IWV difference statistics (ERA-Interim minus GPS). The grey shading shows the surface elevation represented in ERA-Interim, from 0 to 5000 m.

2.2 ERA-Interim reanalysis

ERA-Interim is a modern reanalysis produced by ECMWF using the Integrated Forecasting System (IFS) forecast model and the 4D-Var assimilation system in 12-hourly analysis cycles (Dee et al., 2011). The number of observations has increased from 10⁶ per day in 1989 to 10⁷ per day in 2010. The majority of data, and most of the increase over time, are from satellites. Ground-based GPS data were not assimilated, which make the GPS ZTD and IWV an independent validation dataset. We use ERA-Interim analysis pressure-level data (geopotential, air temperature, and specific humidity) extracted from the Meteorological Archival and Retrieval System (MARS) on a regular latitude–longitude grid with a horizontal resolution of 0.75^∘ × 0.75^∘. For each and every GPS site, 6-hourly ERA-Interim fields are extracted for the four grid points surrounding the GPS station. The IWV contents are computed by integrating the reanalysis specific humidity between the GPS station altitude and the top of atmosphere (1 hPa). For GPS station altitudes located between two pressure levels, the ERA-Interim data at the station are interpolated from the adjacent levels. For stations located below the lowest pressure level (1000 hPa), the reanalysis data are extrapolated. Interpolation and extrapolation are done linearly for specific humidity and temperature, and logarithmically for geopotential, as a function of pressure. To ensure the best spatial matching between GPS and ERA-Interim data, the IWV estimates from the four grid points surrounding the GPS station, IWV₁ to IWV₄, are combined by bilinear interpolation, resulting in a value denoted by IWV_interp. Daily and monthly IWV values are computed afterwards in the same manner as for the GPS IWV data (see above).

2.3 Comparison method

Daily and monthly time-matched IWV values from GPS and ERA-Interim are compared for each and every station, and overall statistics are computed using the full time series (16 years). The overall statistics reveal the systematic or persistent biases and discrepancies between the two datasets. The goal is to identify the main causes of differences among the representativeness differences, errors in the GPS data, and deficiencies in the reanalysis (e.g. in data-sparse regions). The identification of representativeness differences is made by inspection of a number of statistics and their dependence upon characteristics of the GPS station's environment: moist or dry climate (measured by the mean IWV), strength of temporal variability (measured by the standard deviation of IWV and of its first derivative), and spatiotemporal variability of IWV in the vicinity of the station. The latter is computed from the ERA-Interim IWV values at the four grid points surrounding the GPS stations. The maximum absolute deviation of the four IWV values, denoted δ_maxIWV, can reach values as extreme as 18 kg m⁻² in situations of strong large-scale moisture transport (e.g. in the case of tropical plumes reaching the midlatitudes). When averaged over 1 year, the quantity μ_R = mean (δ_maxIWV) is around 2 kg m⁻² for a typical midlatitude station and grows up to 6 kg m⁻² for stations located in regions of complex topography (e.g. station AREQ in the Andes Cordillera). This quantity is referred to as the “representativeness error statistic” in the following.

All the statistics are defined by equations in Appendix A. The values computed for each station are given in Table S2 in the Supplement. They may be useful to readers who want to make their own statistical analysis of our results and/or detect outlying sites based on different thresholds than those we used in Sect. 4.

3 Analysis of the general tendency of IWV differences

The mean and standard deviation of IWV differences (ERA-Interim minus GPS) for all 120 stations over the 16-year period are shown in Figs. 2 to 5. Figure 2 shows the results as a function of station latitude. The general tendency is depicted by the fitted polynomials (the outlying stations, defined beyond the dotted red lines will be discussed in Sect. 4). The different plots show a clear dependence of the results on latitude. The mean difference (Fig. 2a, c) is positive at northern and southern extratropical latitudes (30–80^∘ N and 30–60^∘ S), while it is negative in the intertropical band (30^∘ S–30^∘ N). This result is consistent with the results of Schröder et al. (2016), who compared ERA-Interim to satellite data. The alternation of positive and negative differences is most likely due to biases in the ERA-Interim reanalysis reflecting the difference in moisture information entering the reanalysis over ocean (mainly microwave satellite data) and land (mainly radiosonde and infrared satellite data) (Dee et al., 2011). Indeed, the tropical GPS stations used here are mostly representative of oceanic areas, while the extratropical GPS stations are mainly continental. Similar biases in ERA-Interim were also highlighted by Trenberth et al. (2011), and Parracho et al. (2018), in comparison to other atmospheric reanalyses. The biases remain small, however (below ±1 kg m⁻² or ±10 %). The absolute standard deviation of IWV differences (Fig. 2b) also shows a latitudinal variation with two peaks, around 30^∘ S and 30^∘ N, and dips around the Equator and towards the poles. The equatorial dip is more marked in the relative standard deviation plot (Fig. 2d) because the mean IWV is the largest at these altitudes (∼40 kg m⁻²; see the dashed blue line in Fig. 2c). The enhanced discrepancy between ERA-Interim and GPS daily IWV estimates in the subtropics coincide quite well with the highest day-to-day variability in both hemispheres (see the superposed blue lines in Fig. 2b, d). This strong day-to-day variability is mainly due to the moisture transport associated with the extratropical cyclones in the Northern and Southern Hemisphere storm tracks (Chang et al., 2002; Pfahl et al., 2014). It is not uncommon to observe IWV variations exceeding 20 kg m⁻² d⁻¹ at GPS sites located in the storm track (Bock et al., 2005, 2016). Increased discrepancy between ERA-Interim and GPS at those sites can be due to the imperfect spatiotemporal location of such large moisture variations in the reanalysis or to a representativeness difference between the GPS observations and the reanalysis. No systematic increase in GPS formal error was found in these situations; i.e. the discrepancy is not due to GPS errors.

Figure 2 (a, c) Mean and (b, d) standard deviation of daily IWV difference (ERAI minus GNSS) for 120 global stations as a function of station latitude. Panels (a) and (b) are in kg m⁻²; panels (c) and (d) are in percent of GNSS IWV. The dashed black lines show polynomial fits of orders 5 and 9 for the mean difference and the standard deviation, respectively. The dashed blue lines show polynomial fits of order 7 for (b) the standard deviation of dIWV∕dt (kg m⁻² d⁻¹); (c) the mean IWV (kg m⁻²); (d) the relative standard deviation of dIWV∕dt (% d⁻¹) computed from GPS IWV data. The dotted red lines show the range-check limits used to detect outlying sites (named stations).

Figure 3 Similar to Fig. 2 but plotted as a function of GPS station altitude. The dashed black lines show linear fits.

Figure 3 shows the mean and standard deviation of IWV differences as a function of altitude of the GPS stations. The mean differences (Fig. 3a, c) show no dependence on altitude, meaning that the method of computation of GPS IWV (from ERA-Interim P_s and T_m estimates) and ERA-Interim IWV (from pressure levels) is highly consistent throughout a large altitude range. The standard deviation (Fig. 3b) shows no dependence on altitude either but the relative standard deviation (Fig. 3d) does. The fitted straight line in Fig. 3d shows that this statistic is increasing quite fast as a function of altitude. This tendency can be explained by larger representativeness errors in the reanalysis humidity field as a function of altitude as also found by Waller et al. (2014) in the Met Office high-resolution UK variable resolution model.

Figure 4 shows the standard deviation of IWV differences, σ_Δ, as a function of a few other parameters which give further insight into possible reasons for the discrepancy between GPS and ERA-Interim. Figure 4a and b indicate that, apart from the outliers, there is a moderate tendency for increased discrepancy with increased mean IWV (i.e. warmer and moister climate) and increased IWV variability (including the seasonal variations). The slope of the tendency is actually steeper at the lower IWV bound (mean IWV < 25 kg m⁻²) corresponding to midlatitude and high-latitude sites, while it vanishes at the upper bound, corresponding to intertropical sites (mean IWV ≥ 25 kg m⁻²). The standard deviation of IWV differences reaches a nearly constant level of σ_Δ≈2 kg m⁻² throughout the Equator and the intertropical band. Figure 4c shows that there is a strong tendency for increased discrepancy with increased spatiotemporal variability around the GPS site measured by μ_R (see Sect. 2.3). This interrelation is actually the strongest among all the tested relations between σ_Δ and other statistics. It indicates that representativeness differences are a major source of discrepancy between GPS and ERA-Interim IWV estimates. Finally, Fig. 4d shows that there is only a small tendency for increased discrepancy with increased GPS formal errors.

Figure 4 Standard deviation of daily IWV difference (ERAI minus GNSS) for 120 global stations, as a function of (a) mean GPS IWV; (b) standard deviation of GPS IWV; (c) mean spatial variability of ERAI IWV from the four grid points surrounding the GPS sites used as representativeness statistics (see text); (d) formal error of GPS IWV estimates. Only three outlying stations are named on these plots for clarity.

Figure 5 shows that time averaging is a means of reducing the representativeness differences, as smaller-scale local features captured by the GPS point observations get smoothed out. The mean differences (Fig. 5a, c) are not impacted by the averaging, as expected. The standard deviation of differences (Fig. 5b, d), on the other hand, decreases for the monthly averages, both in absolute and relative units, at all sites. The median standard deviation of the daily IWV differences (ERA-Interim minus GPS) is 1.2 kg m⁻², while the value for the monthly series is 0.51 kg m⁻². The reduction of standard deviation due to averaging is 2.35, which is smaller than the value of $\sqrt{\mathrm{30}}=\mathrm{5.48}$ that one would expect with independent normally distributed data (when averaging over a mean month of 30 d). This inconsistency can be due to the serial correlation in the IWV differences revealing a dependence of the IWV differences upon the meteorological situation. This point might be further investigated by, e.g. separating the IWV differences in different weather regimes. Another means of reducing the discrepancy due to representativeness differences is to use a reanalysis with higher spatial resolution and improved physics representing the smaller-scale atmospheric processes. We compared, for instance, daily GPS IWV data to the Applications of Research to Operations at Mesoscale (AROME) West-Mediterranean operational analysis of Météo-France (this model has a horizontal resolution of 2.5 km × 2.5 km) and found a median standard deviation of difference of 0.81 kg m⁻² over a period of 2 months; we used the GPS and AROME data from the Hydrological cycle in the Mediterranean experiment (HyMeX) Special Observing Period, 5 September–6 November 2012, described in Bock et al. (2016). The median standard deviation of GPS-ERAI differences over 12 stations in the same domain amounts to 0.98 kg m⁻², so there is clear benefit of higher resolution and more modern physics.

Figure 5 Mean vs. standard deviation of IWV difference (ERAI minus GNSS) for (a, b) daily values and (c, d) monthly values. The median values of mean and standard deviation over all 120 stations are 0.47 and 1.2 kg m⁻² (3.1 % and 8.3 %) for the daily results, and 0.47 and 0.51 kg m⁻² (3.1 % and 3.8 %) for the monthly results, respectively. The dotted red lines show the range-check limits used to detect outlying sites (named stations) in the case of the daily data.

Since representativeness differences impose a strong limitation on the agreement between GPS and reanalysis, one may wonder if the horizontal interpolation from the four surrounding ERA-Interim grid points does not further enhance the differences by mixing information from the different grid points. We investigated this question by computing the statistics for each of the four surrounding grid points. Figure 6 shows the results in comparison to the results obtained with the bilinearly interpolated IWV values. The comparison of the mean values (Fig. 6a and b) emphasizes large variations in the biases at some stations which will be further discussed in Sect. 4. The slight shift of the ensemble of results below the 1:1 line is reflecting the fact that a majority of sites exhibit small positive bias (0.47 kg m⁻² on average) as already noticed in Fig. 2a, c, which is not due to representativeness differences but rather to the type of assimilated data (see Sect. 3). The comparison of standard deviations (Fig. 6c and d) shows unambiguously that, at almost all sites, the results for the bilinearly interpolated IWV values are better than for any one of the four surrounding grid points (almost all results sit above the 1:1 line). This conclusion holds for 112 out of 120 stations for the absolute standard deviation (Fig. 6c) and 111 out of 120 stations for the relative standard deviation (Fig. 6d). It indicates that the temporal variability represented by the bilinearly interpolated ERA-Interim IWV data matches best the temporal variability observed by the GPS (i.e. better than from the nearest grid point in the horizontal or in the vertical dimension). When monthly IWV data are compared (not shown), the conclusions are similar, though the number of sites of improved results drops to 71 out of 120 (both for absolute and relative standard deviations). The drop confirms again that the representativeness differences can be reduced by the temporal averaging.

Figure 6 Scatter plots of (a, b) mean and (c, d) standard deviation of daily IWV difference (ERAI minus GPS) when ERAI IWV is bilinearly interpolated from four surrounding grid points (x axis) versus the spread of the mean (a, b) or standard deviations (c, d) for the four surrounding grid points (y axis). The spread is plotted as vertical error bars from the minimum to maximum values. In panels (c) and (d), vertical bars extending below the 1:1 line indicate sites where at least one of the surrounding grid points is in better agreement with GPS than the bilinearly interpolated values; the corresponding stations are named and indicated by a black dot. The dotted red lines show the range check limits as in previous figures. In panels (a) and (b), some of the sites with statistics outside the limits indicated by the dotted red lines are named as well.

4 Analysis of outlying sites

In the previous section, we have seen that the general agreement between GPS and ERA-Interim is limited by representativeness differences which are enhanced in regions of strong temporal variability (Fig. 2b, d), at higher altitude (where mainly the relative standard deviation of differences is impacted; Fig. 3d), and at sites where the mean spatial variability at the four surrounding ERA-Interim grid points is large (Fig. 4c). The standard deviation of differences, σ_Δ, is actually well predicted by our representativeness error statistic, μ_R, with a linear correlation coefficient of r(σ_Δμ_R)=0.73. This strong correlation suggests that the outlying sites, i.e. sites with the largest discrepancy, may have enhanced representativeness errors (Fig. 4c). To investigate this idea, we will analyse in more detail these sites here. First, let us define range limits for each of the four statistics of differences to separate the acceptable sites (i.e. those satisfying the following conditions) from the outliers:

$$\begin{array}{ll}{\displaystyle }-\mathrm{1}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}& {\displaystyle }<{\mathit{\mu }}_{\mathrm{\Delta }}<\mathrm{2}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}\\ {\displaystyle }-\mathrm{6}\mathrm{\%}& {\displaystyle }<{\mathit{\mu }}_{\mathrm{\Delta }}^r<\mathrm{12}\mathrm{\%}\\ {\displaystyle }{\mathit{\sigma }}_{\mathrm{\Delta }}& {\displaystyle }<\mathrm{2.1}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}\\ {\displaystyle }{\mathit{\sigma }}_{\mathrm{\Delta }}^r& {\displaystyle }<\mathrm{18}\mathrm{\%}.\end{array}$$

The values of the limits were determined from visual inspection of Figs. 2, 3, and 5, and shown as the dotted red lines in these figures. The choice for the limits is subjective because we think that the results are very unpredictable when analysing a global network (due to the variety of climates, equipment, and reanalysis performance). So it is necessary to visually inspect the results and determine the limits beyond which the results do not look “normal”. We believe this approach is quite robust thanks to the combination of several representations of the results such as those shown in Figs. 2–6 (i.e. function of latitude, altitude, mean vs. SD scatter plots, etc.). This methodology can be safely applied to other datasets. The reason why we chose non-symmetric limits with respect to zero for the mean differences is because the distribution is not centred on zero. The result is a detection of 15 outlying sites, some of which exceed the limits in more than one test: three sites have excessive absolute bias (CFAG, KIT3, and SANT); nine sites have excessive relative bias (CFAG, COSO, DAV1, KIT3, MAW1, MCM4, POL2, SANT, and SYOG); eight sites have excessive standard deviation of differences (AREQ, BLYT, CFAG, DHLG, IISC, KIT3, LONG, and SANT); and nine sites have excessive relative standard deviation of differences (AREQ, CFAG, KIT3, MAW1, MCM4, MKEA, POL2, SANT, and SYOG). Three sites have statistics exceeding the limits in all four tests (CFAG, KIT3, and SANT). Two of these sites (CFAG and KIT3) are also characterized among the largest representativeness error statistics (Fig. 4c). The time series of IWV and IWV differences for four of the worst cases (CFAG, KIT3, MCM4, and SYOG) can be found in Fig. B2 of Parracho et al. (2018).

Figure 7 shows the values of the four comparison statistics for the 15 outlying cases for the bilinearly interpolated ERA-Interim values and also from the values at the four surrounding grid points (ordered by increasing distance to the GPS station). The results are grouped by region as outlying sites appear to form several clusters located in specific areas of the globe (see Fig. 1). In addition to the four statistics (Fig. 7a to d), we included the altitudes of the GPS stations, h_GPS, and of the four surrounding grid points (Fig. 7e). The above-chosen range limits are superposed in Fig. 7a to d, and a range limit for the altitudes is indicated as h_GPS±500 m (Bock et al., 2014).

Figure 7 (a, b) Mean and (c, d) standard deviation of daily IWV difference (ERAI minus GNSS) for 15 outlying sites grouped by region: Andes (CFAG, SANT, AREQ), central Asia (KIT3, POL2), India (IISC), western USA (DHLG, BLYT, LONG, COSO), Hawaii (MKEA), and Antarctica (MCM4, SYOG, MAW1, DAV1). In panels (a) to (d), the black bars show results for the bilinearly interpolated ERA-Interim data, and the grey bars the results for the four surrounding grid points, ordered by increasing horizontal distance from the GPS station. Panel (e) shows the altitudes of the GPS stations (black bar) and the altitudes of the four surrounding grid points (grey bars). The dotted red lines show the acceptable range limits, similar to Fig. 2 for panels (a) to (d), and ±500 m around the GPS station's altitude in panel (e).

AREQ, SANT, and CFAG are all located in the Andes Cordillera, with AREQ (h_GPS=2470 m) and SANT (h_GPS=696 m) on the western flank of the mountain range facing the sea, and CFAG (h_GPS=680 m) on its eastern flank. The local topography peaks above 3000, 4000, and 3000 m within a radius of 100 km from these three sites, respectively. The altitudes of the four surrounding ERA-Interim grid points are very variable (Fig. 7e), and for AREQ (SANT), all (some) of them are exceeding the altitude range limit. At AREQ, absolute and relative standard deviations of the interpolated data slightly exceed the limits, with σ_Δ=2.4 kg m⁻² and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r=\mathrm{21}$ %, while the bias is almost zero. Moreover, most of the statistics at the four surrounding grid points exceed the range limits. There is thus a significant representativeness difference between the four grid points which is not surprising given the steep orography and the very different altitudes of the grid points. Three of the grid points are actually located more than 500 m higher than the GPS station. For these grid points, the validity of the lower pressure-level data can be questioned as the atmospheric variables are extrapolated far below the model's surface. The results at SANT have similar issues with biases again correlated with variations in the model topography. At both sites, issues with the GPS measurements were eliminated by verifying their consistency with collocated DORIS measurements (Bock et al., 2014). Compared to AREQ and SANT, CFAG has much worse results and gets actually the worst statistics of all 15 sites: μ_Δ=5.8 kg m⁻², ${\mathit{\mu }}_{\mathrm{\Delta }}^r=\mathrm{35}$ %, σ_Δ=3.7 kg m⁻², and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r=\mathrm{22}$ %. Contrary to the previous sites, the results for the four grid points are very similar, though the biases vary slightly (from 5.9 to 4.1 kg m⁻²), which suggests that the discrepancy at this site may not be due that much to spatiotemporal variability in the IWV field. Problems with the GPS measurements cannot be excluded at this site and should be checked by comparison with independent observations.

Further insight into the nature of the discrepancies is given by inspection of the seasonal variation of the comparison statistics (Fig. 8) and of the atmospheric environment (Fig. 9). Figure 8 shows that at all three sites, the biases and standard deviations vary over the year, in relation with the variation of the mean IWV (μ_W, Fig. 9a) and the day-to-day variability (σ_W and σ_dW∕dt, Fig. 9b and c, respectively). Both the μ_Δ and σ_Δ are peaking when μ_W is peaking, during the austral summer months. The relative differences, ${\mathit{\mu }}_{\mathrm{\Delta }}^r$ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$, and IWV variability, ${\mathit{\sigma }}_W^r$ and ${\mathit{\sigma }}_{\mathrm{d}W/\mathrm{d}t}^r$, are peaking in winter when the mean IWV is low. Inspection of μ_R (Fig. 9d) confirms the strong impact of spatiotemporal variability at all three sites but especially at AREQ, where it is the largest among all sites (peaking at μ_R=6.4 kg m⁻²). It is noticeable that at CFAG the yearly mean and the seasonal cycle of IWV in ERA-Interim are larger than observed by GPS (Fig. 9a), which suggests that a representativeness difference is most likely the explanation rather than GPS measurement issues evoked above.

Figure 8 Seasonal variation of (a, b) mean and (c, d) standard deviation of daily IWV difference (ERAI minus GNSS) for 15 outlying sites. The grey bars show the statistics computed for each month (January to December, from left to right) over the 16-year period. The dotted red lines show the range check limits.

Figure 9 Seasonal variation of daily IWV data: (a) mean IWV; (b, e) absolute and relative standard deviation of IWV; (c, f) absolute and relative standard deviation of IWV derivative; (d, g) absolute and relative mean spatial variability of ERAI IWV from the four grid points surrounding the GPS sites. The grey (blue) bars show GPS (ERA-Interim) statistics computed for each month (January to December, from left to right) over the 16-year period.

The next two sites, KIT3 (h_GPS=659 m) and POL2 (h_GPS=1755 m), are located in Uzbekistan and Kyrgyzstan, respectively, close to the Alai/Tien Shan mountain range. They both show large difference statistics, with μ_Δ, ${\mathit{\mu }}_{\mathrm{\Delta }}^r$, σ_Δ, and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$ exceeding the limits for KIT3 and ${\mathit{\mu }}_{\mathrm{\Delta }}^r$ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$ for POL2 (Fig. 7a to d). Considering the individual grid points, they almost all also exceed the limits, with large variations both in the bias and standard deviation at KIT3, with somewhat smaller differences at POL2. These variations can again be related to variations in the grid point altitudes, some of which exceed the range limits (Fig. 7e). The difference statistics at these sites exhibit large seasonal variations, with μ_Δ, ${\mathit{\mu }}_{\mathrm{\Delta }}^r$, and σ_Δ peaking in boreal summer (Fig. 8) when μ_W and μ_R are peaking (Fig. 9). The representativeness error statistics peaks are particularly marked at these stations, with KIT3 showing the largest monthly values among all sites (${\mathit{\mu }}_{R,i}=\mathrm{8.8}$ kg m⁻² in August; Fig. 9d). At this site, the GPS measurements were also verified with collocated DORIS measurements (Bock et al., 2014), confirming that representativeness differences between ERA-Interim and GPS IWV data are the main reason for this discrepancy. Interestingly, it can be noticed that the peak in IWV during summer is significantly larger in ERA-Interim compared to GPS (Fig. 9a), suggesting excessive moisture transport into this region in the reanalysis, possibly connected with the too-smooth topography in the model.

The next five sites belong to two geographical regions: IISC in India, and DHLG, BLYT, LONG, and COSO in California, USA, which are all characterized by small discrepancies with only one statistic exceeding the range limits (σ_Δ for the first four, and ${\mathit{\mu }}_{\mathrm{\Delta }}^r$ for COSO). At all five sites, the variation of statistics among the four grid points are small (Fig. 7a to d), as are the variations of the altitudes (Fig. 7e). Station IISC shows a small seasonal variation in the bias and standard deviation (Fig. 8) which might be linked to the variation in IWV temporal variability (Fig. 9b, c, e, f) and spatiotemporal variability (Fig. 9d and g) that show peaks in spring and autumn, i.e. during transitions seasons between the summer monsoon (June to October) and the cooler winter season (December to March). It has been shown previously that monsoon transition periods are accompanied by strong spatial and temporal variability in IWV, which is difficult to represent in atmospheric reanalyses (Bock et al., 2008; Bock and Nuret, 2009; Meynadier et al., 2010; Means, 2013).

The four outlying Californian sites can be separated into two groups: DHLG, BLYT, and LONG, located south of the Sierra Nevada mountain range, in a region of moderate topography, and COSO located in the Basin and Range Province, a narrow valley at the southern exit of the Sierra Nevada. The higher altitude (1485 m) and more complex topographic environment of COSO enhances the representativeness differences. Interestingly, all four sites show a step-like variation of the mean IWV and variability (Fig. 9a, b, c) peaking in July–August–September associated with the North American monsoon (Adams and Comrie, 1997; Means, 2013). This feature contrasts with the Indian monsoon observed at IISC where variability was enhanced during the transition seasons and not during the monsoon. At DHLG and BLYT, the biases actually reverse signs in July–August (Fig. 8a, b) and the standard deviation peaks at σ_Δ>4 kg m⁻² (Fig. 8c). Figure 9b and c show that ERA-Interim underestimates IWV variability at these sites which suggests that GPS observations capture some small-scale moisture variability not represented in the reanalysis.

The next site, MKEA (h_GPS=3730 m), is located on the Mauna Kea volcano on the island of Hawaii. Due to small size of the emerged land area (approximately 10⁴ km²), the imprint of the island is almost inexistent in the reanalysis topography (Fig. 7e). Hence, it is not surprising that the comparison statistics are bad (although only ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$ is exceeding the range limits). The relative IWV differences are huge (Fig. 9e and f) when computed with respect to the low GPS IWV content of this high-altitude site.

The last group of sites is located in eastern Antarctica (Fig. 1). Unfortunately, four of the five Antarctica sites used in this study suffer from large discrepancies. Three of them have two statistics (${\mathit{\mu }}_{\mathrm{\Delta }}^r$ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r)$ exceeding the range limits (Fig. 7b and d). MCM4 is the worst case and has the largest relative standard deviation among all 15 sites: ${\mathit{\sigma }}_{\mathrm{\Delta }}^r=\mathrm{32}$ %. This station is located in McMurdo Sound, an area with complex landscape, including local low mountain peaks, valleys and glacier corridors, and sea within a radius of 100 km. The other three stations are located close to the coastline backed to the main ice shelf with large surface elevation variations (up to 2000 m within a distance of 100 km). The grid points in ERA-Interim are associated at different altitudes with differences in representativeness leading to IWV biases (Fig. 7b). The marked seasonal variation of ${\mathit{\mu }}_{\mathrm{\Delta }}^r$ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$ (Fig. 8b and d) also confirms a dependence of the IWV differences on the atmospheric state and especially on IWV variability, which is enhanced during the austral winter months (Fig. 9e and f). The winter variability is actually much underestimated in ERA-Interim as seen in Fig. 9e and f at MCM4, SYOG, and MAW1, and, quite surprisingly, the spatiotemporal variability, ${\mathit{\mu }}_R^r$, remains nearly constant in ERA-Interim (Fig. 9g). These differences point to an issue in ERA-Interim IWV contents in Antarctica, especially during austral winter, as also suggested by Parracho et al. (2018), who compared ERA-Interim to the NASA Modern Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) reanalysis. These authors also pointed to some issues in the GPS measurements at MCM4 and SYOG between 2002 and 2006, as well as a break in the IWV series at all sites in Antarctica due to a discontinuity in the GPS processing. The IWV issues in ERA-Interim may be linked to the large surface air temperature biases of the reanalysis diagnosed by Bracegirdle and Marshall (2012), from coastal station observations which are related to its too-smooth orography. In addition, Xie et al. (2016), showed that the replicability of daily and annual variance of surface air temperature in this reanalysis decreases from the coast to the interior of the continent. These results also support the findings of Parracho et al. (2018) that the IWV variability and trends in ERA-Interim reanalysis are more realistic near the coast where in situ observations are assimilated than in the interior where the reanalysis mainly relies on satellite observations and short-term model forecasts. Representativeness differences between GPS and ERA-Interim in Antarctica are thus enhanced by deficiencies in the reanalysis.

5 Discussion and conclusions

In this study, we first analysed the general tendency of IWV difference between ERA-Interim reanalysis and global GPS observations. We found that the mean difference, interpreted as the bias of the reanalysis with respect to the observations, exhibits a latitudinal variation of ±1 kg m⁻², consistent with the fact that different moisture information is entering the reanalysis over ocean and land. As a result, the northern and southern midlatitudes exhibit a moist bias, while the tropics are too dry. This bias is not changing with the altitude of the observation site. The standard deviation of daily IWV differences is generally below 2 kg m⁻² but peaks at the northern and southern storm-track latitudes. This result suggests that GPS observations capture some small-scale variability that is not resolved by the reanalysis. Another indication that the discrepancies are process related is that the relative standard deviation is increasing with altitude (from about 8 % at sea level to 16 % at 2.5 km). More generally, it is shown that discrepancies are due to representativeness differences between the gridded reanalysis field and the GPS point observations. A strong correlation (r=0.73) is found between the standard deviation of IWV differences, σ_Δ, and our representativeness error statistic, μ_R, which measures the spatiotemporal variability in the vicinity of the GPS site based on the analysis of the ERA-Interim IWV data at the four surrounding grid points. However, it is shown that in general (for 112 sites out of 120), bilinearly interpolated IWV values from the four surrounding grid points are in better agreement with the GPS observations than any of the grid points individually. Even if the horizontal resolution of the reanalysis grid is quite coarse (0.75^∘ × 0.75^∘), spatial interpolation does not reduce the representativeness. It is also shown that the standard deviation of IWV differences is further reduced when data are time averaged (e.g. in monthly bins). Indeed, spatial and temporal averaging smooths out the variability due to smaller-scale phenomena and makes the reanalysis and observations more consistent at representing the larger-scale meteorological systems.

In a second part, we analysed in more detail the possible reasons for the very bad comparison results obtained at 15 outlying sites. It is shown that at most of the sites, representativeness errors are the most plausible cause for discrepancies which are enhanced because of local topographic and climatic features. The problematic topographic features include steep orography such as that found for sites in the Andes Cordillera (AREQ, CFAG, and SANT), on the island of Hawaii (MKEA), and close to the Himalaya chain (KIT3 and POL2), as well as coastal sites in Antarctica (MCM4, SYOG, MAW1, and DAV1). The climatic features include large seasonal changes in the total IWV, such as those associated with the Indian monsoon (IISC, KIT3, POL2) or the North American monsoon (DHLG, BLYT, LONG, and COSO), and/or in the IWV synoptic variability (observed at most sites during either the transition seasons, winter or summer, depending on the geographic location). When these 15 stations are eliminated from the dataset, the comparison statistics become ${\mathit{\mu }}_{\mathrm{\Delta }}=\mathrm{0.36}\pm \mathrm{0.49}$ kg m⁻², ${\mathit{\mu }}_{\mathrm{\Delta }}^r=\mathrm{2.7}\pm \mathrm{3.5}$ %, ${\mathit{\sigma }}_{\mathrm{\Delta }}=\mathrm{1.22}\pm \mathrm{0.38}$ kg m⁻², and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r=\mathrm{8.2}\pm \mathrm{3.0}$ % (mean ± standard deviation over the 105 sites). They are significantly improved compared to the initial results including the 120 sites: the standard deviations of μ_Δ and ${\mathit{\mu }}_{\mathrm{\Delta }}^r$ are reduced by 30 %, the means of σ_Δ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$ by 20 % and the standard deviations of σ_Δ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$ are reduced by 40 %. Because the comparison of GPS and ERA-Interim is not relevant at these sites, we recommend not to use ERA-Interim in the homogenization process of these GPS time series (Ning et al., 2016; Van Malderen et al., 2017) unless the homogenization method explicitly models the seasonality in the bias and the non-stationarity of the noise variance (Quarello et al., 2018).

These results lead to a more general question of whether it is possible to eliminate problematic stations a priori, i.e. before the comparison statistics are computed. Inspection of the elevation of the four surrounding grid points with respect to the elevation of the GPS station and with respect to each other provides some indication of possible representativeness errors. Some correlation between IWV biases and altitudes at the individual grid points was found in extreme cases (Fig. 7). A simple a priori check based on the comparison of grid point altitudes to station altitudes would eliminate some of the problematic cases. We compared the statistics with and without selection of sites where the elevation of the grid points differs by more than 500 m from the GPS station. When the selection is applied to the nearest grid point only, 15 stations are eliminated, including four of the outlying sites discussed in Sect. 4. This test is not very efficient. When applied to all four surrounding grid points, 34 stations are eliminated, including 11 of the outlying sites (only CFAG, MCM4, BLYT, and IISC remain then in the dataset). On average, the statistics of the mean differences (μ_Δ and ${\mathit{\mu }}_{\mathrm{\Delta }}^r$) do not change very much in that case, mainly because the stations with the largest absolute and relative biases (CFAG and MCM4) are not eliminated. However, the statistics of the standard deviation of differences (σ_Δ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$) are reduced by about 20 %. However, the benefit is at the expense of a strong reduction of the number of sites (34 stations eliminated). Though altitude differences have been shown to explain discrepancies at certain stations a posteriori in Sect. 4, this altitude check appears too excessive to be applied in a systematic way a priori. We also tested the use of the absolute and relative representativeness error statistics, μ_R and ${\mathit{\mu }}_R^r$, and found that a threshold of 20 % on ${\mathit{\mu }}_R^r$ eliminates 13 stations, including 8 out of the 15 outlying sites, and reduces the error statistics σ_Δ and ${\mathit{\sigma }}_{\mathrm{\Delta }}^r$ by 20 % to 30 % on average. This outlier check is efficient and is thus recommended. However, none of the checks that we tested were able to detect all the 15 outlying sites. Hence, it is also advised to carefully analyse the comparison statistics in order to understand the possible causes of discrepancies and eliminate outlying stations a posteriori on a subjective basis as we have done in this study. This was possible here because the number of stations was small. In more extended networks, an automatic selection method based on, e.g. a clustering algorithm, would be necessary.

Aside from the large representativeness errors found at a small number of sites, one should recognize that ERA-Interim and GPS IWV data are generally in good agreement globally, except perhaps in Antarctica where the comparison failed at four sites out of five. One of the remaining error sources not addressed in this study is the temporal consistency of both data sources. Therefore, other statistics are more relevant, such as trend estimates (Schröder et al., 2016; Parracho et al., 2018). The methodology described in this paper can also be applied to assess the consistency and representativeness of other data sources (e.g. climate models, satellite IWV data) and other observation types (e.g. surface humidity, temperature).

Appendix A: Definition of variables and comparison statistics

Throughout this study, the GPS IWV data at a given station are denoted by IWV_GPS and the corresponding ERA-Interim IWV data are denoted IWV_ERAI. When the subscript is not specified, the IWV data may refer interchangeably to GPS and ERA-Interim. When the ERA-Interim IWV data from four surrounding grid points need be distinguished, the subscript i is added, with i=1.4, and the bilinearly interpolated value is then denoted by IWV_interp.

GPS and ERA-Interim IWV data are analysed using the following statistics, where the mean and standard deviation are computed over the number of days (months) of the time-matched daily (monthly) data:

• The mean and standard deviation of IWV are

  $$\begin{array}{}\text{(A1)}& {\displaystyle }{\mathit{\mu }}_W=\mathrm{mean}\left(\mathrm{IWV}\right)\text{(A2)}& {\displaystyle }{\mathit{\sigma }}_W=\mathrm{SD}\left(\mathrm{IWV}\right).\end{array}$$

• The relative standard deviation of IWV is

  $$\begin{array}{}\text{(A3)}& {\mathit{\sigma }}_W^r={\displaystyle \frac{\mathrm{SD}\left(\mathrm{IWV}\right)}{\mathrm{mean}\left(\mathrm{IWV}\right)}}.\end{array}$$

• The standard deviation and relative standard deviation of IWV time derivate as

  $$\begin{array}{}\text{(A4)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{d}W/\mathrm{d}t}=\mathrm{SD}(\mathrm{dIWV}/\mathrm{d}t)\text{(A5)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{d}W/\mathrm{d}t}^r={\displaystyle \frac{\mathrm{SD}(\mathrm{dIWV}/\mathrm{d}t)}{\mathrm{mean}\left(\mathrm{IWV}\right)}}.\end{array}$$

The ERA-Interim representativeness error statistic is based on the maximum absolute difference in IWV from the four surrounding grid points, ${\mathit{\delta }}_{max}\mathrm{IWV}=\left({\mathrm{IWV}}_{\mathrm{ERAI},i}\right)-\left({\mathrm{IWV}}_{\mathrm{ERAI},i}\right)$:

• The absolute and relative mean “representativeness error statistic” is

  $$\begin{array}{}\text{(A6)}& {\displaystyle }{\mathit{\mu }}_R=\mathrm{mean}\left({\mathit{\delta }}_{max}\mathrm{IWV}\right)\text{(A7)}& {\displaystyle }{\mathit{\mu }}_R^r={\displaystyle \frac{\mathrm{mean}\left({\mathit{\delta }}_{max}\mathrm{IWV}\right)}{\mathrm{mean}\left({\mathrm{IWV}}_{\mathrm{ERAI},\mathrm{interp}}\right)}}.\end{array}$$

The ERA-Interim minus GPS differences are analysed using the following statistics:

• The mean and standard deviation of IWV differences are

  $$\begin{array}{}\text{(A8)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{\Delta }}=\mathrm{mean}({\mathrm{IWV}}_{\mathrm{ERAI}}-{\mathrm{IWV}}_{\mathrm{GPS}})\text{(A9)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{\Delta }}=\mathrm{SD}({\mathrm{IWV}}_{\mathrm{ERAI}}-{\mathrm{IWV}}_{\mathrm{GPS}}).\end{array}$$

• The relative mean and standard deviation of IWV differences are

  $$\begin{array}{}\text{(A10)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{\Delta }}^r={\displaystyle \frac{\mathrm{mean}({\mathrm{IWV}}_{\mathrm{ERAI}}-{\mathrm{IWV}}_{\mathrm{GPS}})}{\mathrm{mean}\left({\mathrm{IWV}}_{\mathrm{GPS}}\right)}}\text{(A11)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{\Delta }}^r={\displaystyle \frac{\mathrm{SD}({\mathrm{IWV}}_{\mathrm{ERAI}}-{\mathrm{IWV}}_{\mathrm{GPS}})}{\mathrm{mean}\left({\mathrm{IWV}}_{\mathrm{GPS}}\right)}}.\end{array}$$

In Eqs. (A6) to (A9), the ERA-Interim IWV values can be IWV_interp (as in Figs. 2 to 5) or any one of the IWV_ERAI,i when individual grid points are discussed (as in Figs. 6 and 7). In Sect. 4 of the paper, statistics from Fig. 7 are referred to as μ_R,i, ${\mathit{\mu }}_{R,i}^r$, etc. when the representativeness error estimates from individual grid points are discussed.

The units of the values computed using Eqs. (A1), (A2), (A6), (A8), and (A9) are given in kg m⁻².

The units of the values computed using Eqs. (A3), (A10), and (A11) are percentages when multiplied by 100.

The units of the values computed using Eq. (A4) are given in kg m⁻² d⁻¹ and using Eq. (A5) they are % d⁻¹ when multiplied by 100.