A numerical method to improve the spatial interpolation of water vapor from numerical weather models: a case study in South and Central America

. Commonly Numerical Weather Models (NWM) users can get the vertically Integrated Water Vapor (IWV) value at a given location from the values at nearby grid points. In this study we used a validated and free available Global Navigation Satellite Systems (GNSS) IWV data set to analyze the very well-known effect of height differences. To this aim, we studied the behavior of 67 GNSS stations in Central and South America with the condition of having a minimum of 5 years of data during the period from 2007 till 2013. The values of IWV from GNSS were compared with the respective values from ERA 5 Interim and MERRA-2 in the same period. Firstly, the total set of stations was compared in order to detect in which cases the geopotential difference between GNSS and NWM deserves a correction. Then, an additive integral correction to the IWV values from ERA Interim was proposed. For the calculation of this correction the multilevel values of speciﬁc humidity and temperature given at 37 pressure levels by ERA Interim were used. The performance of the numerical integration method was tested by accurately reproducing the IWV values at each of the grid points surrounding each of the GNSS sites under study. 10 Finally, and considering the IWV GNSS values 6 as a reference, the improvement introduced to the IWV ERAInterim values after adding the corrections is analyzed. In general, the corrections are always recommended but they are not advisable at sea coastal areas or on islands since at least two grid points of the model are usually in the water. In such cases the additive correction could overestimate the IWV.


Introduction
Water vapor is an abundant natural greenhouse gas of the atmosphere.The knowledge of its variability in time and space is very important to understand the global climate system (Dessler et al., 2008).Most of the regional comparisons of IWV from GNSS are aimed at validating the technique by comparing with radiosonde and microwave radiometers where available.An example of this is the work of Van Malderen et al. (2014) who compared IWV GPS (Global Positioning System) with IWV derived from ground-based sun photometers, radiosondes and satellite-based values from GOME, SCIAMACHY, GOME-2 and AIRS instruments at 28 sites in the northern hemisphere.Because their comparison is oriented to climatology application, they deal with long-term time series (+ 10 years).The authors asseverate that the mean biases of the GPS with the different instruments vary only between -0.3 and 0.5 kg m −2 but there are large standard deviations especially for the satellite instruments.
However, some other comparisons examine the IW V GN SS values with respect to the respective estimates from Numerical Weather Models (NWM).If focusing on the application of the current state-of-the-art reanalysis ERA-Interim from the European Centre for Medium-Range Weather Forecasts (ECMWF), both in local and global scale, some recent papers deserve to be mentioned: Heise et al. (2009) used ground pressure data from ECMWF to calculate IWV from 5-minutes Zenith Total Delay (ZTD) at stations without meteorological data available.The authors validate their results with stations with local measurements of pressure and temperature.They also compare IWV from GPS with respect to IWV from ERA-Interim on a global scale.The authors found that IWV from GPS and ECMWF show well agreement on most stations on the global scale except in mountain regions.Moreover they addressed that temporal station pressure interpolation may result in up to 0.5 kg m −2 IWV uncertainty if a local weather event happened.According to the authors, this phenomenon is observed especially in the tropics and is due to the fact that the ECMWF analysis does not adequately represent the local situation if facing with an increase in the diurnal cycle of surface atmospheric pressure.Buehler et al. (2012) compare IWV values over Kiruna in the north of Sweden from five different techniques (radiosondes, GPS, ground-based Fourier-Transform InfraRed (FTIR) spectrometer, ground-based microwave radiometer, and satellite-based microwave radiometer) with IWV from ERA-Interim reanalysis.The processed GPS dataset covers a ten-year period from November 1996 to November 2006.The authors found a good overall agreement between IWV from GPS and from ERA-Interim being the mean of differences 0.29 kg m −2 and its standard deviation 1.25 kg m −2 .They also point out that ERA-Interim is drier than the GPS at small IWV values and slightly moister at high IWV values (above 15 kg m −2 ).The authors also consider altitudes limits when comparing measurements from different data sets.They proposed an empirical solution by computing linear regression slopes as a function of the height and corrected all measurements to a common reference altitude of 430 meters.Thus, they established a relative bias of -3.5 % per 100 meters that introduced absolute errors below 0.2 kg m −2 .Nevertheless they asseverate that the good performance of the method depends on location and probably on season.Ning et al. (2013) evaluate IWV from GPS in comparison with IWV from ERA-Interim and IWV from the regional Rossby Centre Atmospheric (RCA) climate model at 99 European sites for a 14-year period.Because RCA is not an assimilation model, the standard deviation of RCA-GPS resulted 3 times larger than ERA-GPS.The IWV difference for individual sites varies from -0.21 up to 1.12 kg m −2 for ERA-GPS and the corresponding standard deviation is 0.35 kg m −2 .They investigate the influence of the differences between NWM and GPS in the vertical and horizontal positions.In particular, the authors studied sub-sets of stations with absolute value of height differences smaller than 100 meters.Consequently they do not take into account if there were an over or underestimation by the models.Thus, they found values of the monthly mean IWV differences smaller than 0.5 kg m −2 .Moreover, the authors also highlight that the models overestimate IWV for sites near the sea.Bordi et al. (2014) studied global trend patterns of a yearly mean of IWV from the 20 th century atmosphere model (ERA-20CM) and ERA-Interim, both from ECMWF.The authors highlight a regional dipole pattern of inter-annual climate variability over South America from ERA-Interim data.According to this study, the Andean Amazon basin and Northeast Brazil are characterized by rising and decreasing water content associated with water vapor convergence (divergence) and upward (downward) mass fluxes, respectively.Besides, the authors also compared IWV from ERA-Interim with the values estimated at 2 GPS stations in Bogotá and Brasilia.Such comparison on monthly timescale made known a systematic bias attributed to a lack of coincidence in the elevation of the GPS stations and the model grid points.Tsidu et al. (2015) presented a comparison between IWV from a Fourier Transform InfraRed spectrometer (FTIR, at Addis Ababa), GPS, radiosondes, and ERA-Interim over Ethiopia for the period 2007-2011.The study is focused on the characterization of the different error sources affecting the data time series.In particular, from the study of diurnal and seasonal variabilities, the authors addressed differences in the magnitude and sign of the IWV bias between ERA-Interim and GPS.
They linked this effect with the sensitivity of the convection model with respect to the topography.Wang et al. (2015) performed a 12-year comparison of IWV from 3 third generation atmospheric reanalysis models including ERA-Interim, MERRA and the Climate Forecast System Reanalysis (CFSR) on a global scale.IWV values from the reanalysis models were also compared with radiosonde observations in land and Remote Sensing Systems (RSS) on satellites over oceans.
The authors asseverate that the main discrepancies of the 3 datasets among them are in Central Africa, Northern South America, and highlands.
In this paper, we investigate the differences between IWV from GNSS by using data products from (Bianchi et al., 2016a) and IWV values given by ERA-Interim and MERRA-2.The comparison was performed taking into account the geopotential differences (∆z) between each GNSS station and the corresponding values assigned by the models.We proposed an additive numerical correction to the IWV from NWM and the strategy was tested for ERA-Interim re-analysis model.Section 2 describes the different sets of data used in this study.Following, we give the explanation of the methodology and the presentation of the results obtained after applying the proposed correction to the IWV values from ERA-Interim.

IWV from GNSS
The GNSS data is the main source of information for the spatial and temporal distribution of water vapor.Thus, the main variable considered is the IWV estimated from the delay caused by the troposphere to the GNSS radio signals during its travel from the satellite to the ground receiver.The total delay projected onto the zenith direction (ZTD) is usually split into two contributions: the hydrostatic delay (ZHD, Zenith Hydrostatic Delay) depending merely on the atmospheric pressure and the Zenith Wet Delay (ZWD) depending mainly on the humidity.The IW V GN SS can be obtained from ZWD multiplying it by a function of the mean temperature of the atmosphere.
The reference database of IW V GN SS used in this study comes from a geodetic processing over 136 tracking stations in the American continent located from southern California to Antarctica (see Figure 1), during the 7-year period from January 2007 till December 2013 (Bianchi et al., 2016b).Specifically, the data series of IW V GN SS used here is restricted to those 67 stations with IWV time series spanning more than 5 years (see Figure 1 and Table 1).More details of the steps, models and conventions followed by the geodetic processing to obtain the IW V GN SS values are in Bianchi et al. (2016a).

IWV from NWM
The values of columnar Integrated content of Water Vapor (IWV) as reanalysis products from ERA-Interim (Dee et al., 2011) and MERRA-2 (GMAO, 2015;Bosilovich et al., 2015;Gelaro et al., 2017) were evaluated in this study.The horizontal resolutions are 0.25 • × 0.25 • for ERA-Interim and 0.625 • × 0.50 • for MERRA-2.Because ERA-Interim data is given 4 times a day, in order to perform the comparison and even if MERRA-2 gives hourly data, we pick up IWV data from MERRA-2 every 6 hours at 0, 6, 12 and 18 hours of Universal Time (UT).Thus, to be able to carry out the comparison, MERRA-2 was only partially evaluated.
ERA-Interim is the global atmospheric reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF).It covers the period from 1979 up to today and supersedes the ERA-40 reanalysis.ERA-Interim address some difficulties of ERA-40 in data assimilation mainly related to the representation of the hydrological cycle, the quality of the stratospheric circulation, and the consistency in terms of reanalyzed geophysical fields (Dee et al., 2011).
MERRA-2 is the successor of The Modern-Era Retrospective Analysis for Research and Applications (MERRA) from NASA's Global Modeling and Assimilation Office (Rienecker et al., 2011).MERRA-2 is improved because it contains less trends and jumps linked to changes in the observing systems than MERRA.Additionally, MERRA-2 assimilates observations not available to MERRA and reduces bias and imbalances in the water cycle (Gelaro et al., 2017).Moreover, the longitudinal resolution of MERRA-2 data is changed from 0.667 • in MERRA to 0.625 • whereas the latitudinal resolution remains unchanged (0.5 • ) (Bosilovich et al., 2015).
To this application we used two different kind of data sets: The 2-D values of the IWV from both re-analysis models along with the correspondent geopotential invariant.We also used three 3-D data sets from ERA-Interim: the air temperature (T ), the specific humidity (q) and the geopotential (z).These variables are given in 37 levels of atmospheric pressure from 1 to 1000 hPa.
We will denominate as z i , T i , p i and q i to the value of the before mentioned variables at a given i level and in a given instant .This set of data will be used for the calculation of the integral correction that will be developed in the following section.

Methodology
As mentioned before, even when both reanalysis models give grid values of the vertical integral of the water vapor, the solution provided by each model is linked to its respective geopotential surface invariant.Nevertheless, elevation differences between geopotential from each model grid and computed from GNSS height must be addressed.
In order to compute the geopotential of the GNSS stations (z GN SS ) we followed the van Dam et al. (2010) algorithm.
Because the geodetic coordinates (φ, λ, h) of the GNSS site are known, we obtained the orthometric height (H) at each GNSS station by correcting the ellipsoidal height with the EGM08 model (Pavlis et al., 2012).Thus, the geopotential is (van Dam et al., 2010) where the radius of the ellipsoid at geodetic latitude φ is, (2) with a = 6378137m.and b = 6356752.3142m.being the semimajor and semiminor axis of the WGS84 ellipsoid, respectively (Hofmann-Wellenhof and Moritz, 2006).Moreover, the value of the gravity on the ellipsoid at geodetic latitude φ can be written as (van Dam et al., 2010).
where e 2 = 0.00669437999014 is the first eccentricity squared of the WGS84 ellipsoid and g E = 9.7803253359m s −2 is the normal gravity at the Equator (Hofmann-Wellenhof and Moritz, 2006) and k s = 0.001931853 (van Dam et al., 2010).
For a given GNSS station, the respective geopotential from each of the 2 reanalysis models resulted from a bi-linear interpolation of the respective invariant static geopotential at the 4 grid points around the GNSS site, referred to as Because the points of the NWM grid surrounding the GNSS station have different geopotential values and those values are in turn different from z GN SS , we propose to correct on each of these 4 points.Thus, if ∆z k refers to the difference between z GN SS and z k N W M , where NWM corresponds to ERA-Interim or MERRA-2.
We will then propose a correction procedure that, compensating ∆z k in each of the 4 grid points, corrects the values of IW V N W M .After "moving" the grid points to z GN SS , a bi-linear interpolation is performed to obtained the corrected value of IWV at the location of the GNSS site.In brief, this procedure is equivalent to lifting (or descending as appropriate) each of the grid points in order to build up a plane at z GN SS .
Prior to the correction, we analyze the performance of ERA-Interim and MERRA-2 with respect to GNSS.Thus, although IW V GN SS is produced every half hour and IW V M ERRA−2 is available hourly, we consider only the epochs when ERA Interim data is available to perform the comparison at 0, 6, 12 and 18 UT.Table 2 shows the mean values of IWV from GNSS (IW V GN SS ) during the period 2007-2013 and its standard deviation for all the stations.We assume that the static geopotential from the NWM at each GNSS site (z N W M ) is obtained from a bi-linear interpolation of the static geopotential at the 4 grid points surrounding it (z k N W M ).Thus, Table 2 shows the geopotential difference (∆z = z GN SS − z N W M ) in addition to the respective differences of the mean values (∆IW In general when regarding at Table 2, we can observe that the best agreements between the average IWV values from GNSS and the corresponding average from the models are where the ∆z are small (for example: CONZ, VITH, SMAR, LPGS, MAPA, SCUB among others).In other words, in general the NWMs represent very well the IWV values (∆IW V < 1.5 kg m −2 ) if |∆z| is small.That means, the geopotential difference is in the order of 500 m 2 s −2 at most.
On the other hand, the difference of the model representation of the IWV with respect to GNSS grows as the height differences (∆z) become larger and this is true for all values of IW V GN SS .As an example we can mention the cases of SANT However, other than these cases that can be considered critical, the differences are also important in those sites with moderate |∆z| (larger that 500 m 2 s −2 ) and IW V GN SS > 20kg −2 (CEFE, BRAZ, RIOD, GUAT).
Note that some MERRA-2 differences values could be a little bigger than ERA Interim ones and this would be expected because of the coarser grid.However, this is not a general rule and some stations are in fact better represented by MERRA-2 with |∆z| and ∆IW V smaller than ERA Interim even if they are located in highlands (see SANT, COPO).
Figure 2 (up) shows the the mean IWV values from GNSS (IW V GN SS ) as a function of geopotential differences (∆z).
Results for MERRA-2 are on the left and the same for ERA-Interim on the right.
It is assumed that these different values of ∆IW V are due to the geopotential difference.Therefore, they commonly carry the opposite sign to the ∆z as we expect.Effectively, ∆IW V is nothing but the difference between mean values from GNSS and NWM, a negative value for negative ∆z indicates an overestimation by the reanalysis model and on the contrary, the underestimation by the model is shown with red dots where ∆z is positive.
However the figure 2 (up) shows that there are some cases where ∆IW V have the same sign as ∆z evidencing that the reanalysis models can overvalue or undervalue the value of IW V and this should not be due to ∆z.This effect can be seen within a rectangle whose size are ∆z = ± 2000 m 2 s −2 and ∆IW V = ± 2 kg m 2 .Such a cloud of points represent about 21% of the total number of stations for ERA-Interim and a 27% for MERRA-2.Figure 2 (down) shows the spatial distribution of ∆IW V for MERRA-2 (left) and ERA Interim (right).The red dots indicate positive differences greater than 2.5 kg m 2 and the intermediate differences between 0 and 2.5 kg m 2 are in a red degrade.Similarly, the blue dots show negative differences less than -2.5 kg m 2 and the differences between -2.5 and 0 kg m 2 are in a blue degrade.
If we take as reference the RNNA station in both maps (see station number 51 in Figure 1), and advancing towards the south along the Atlantic coast, the behavior of both models results similar.Both reanalysis are dryer than GNSS and this same effect is seen in the southern mountainous areas.However, going along the Atlantic coast from RNNA to the north and up to the Amazon river we can see different behaviours of the reanalyses, while ERA Interim continues underestimating IW V GN SS , MERRA-2 shows to be wetter than GNSS.The overall agreement between GNSS and MERRA-2 is -0.39 ± 2.77 and between GNSS and ERA Interim is 0.13 ± 2.52.Being these values the result of an average over all the ∆IW V differences.
These findings show that MERRA-2 resulted wetter than GNSS while ERA Interim is slightly dryer than GNSS in Central and South America.This is in agreement with the findings of Buehler et al. (2012), who found out that the mean value of the differences between IWV from GPS and ERA is 0.28 ± 1.25 for a high latitude location in Sweden, revealing also an underestimation of the reanalysis model.
Finally, the correlation coefficients between IW V GN SS values and the respective ones for both NWM, are higher than 0.95 in most of the GNSS stations (not shown).

Computation of the integral correction
Following we proceed to calculate a correction in order to provide a better estimation of the IW V N W M at the GNSS site.
We start by correcting each of the grid points around the station prior to apply a bi-linear interpolation.The correction will be computed only for one of the two re-analysis models tested.We have chosen ERA Interim over MERRA 2 for the calculation and testing of these corrections not only because ERA Interim has a thinner grid, but also considering the results of Zhu (2014).
Effectively, Zhu (2014) compared several reanalysis projects with independent sounding observations recorded in the Eastern Himalayas during June 2010.Among all the reanalysis models ERA-Interim and MERRA, the predecessor of MERRA 2 were included.The authors analyze temperature, specific humidity, u-wind, and v-wind between 100 hPa and 650 hPa.The authors found that ERA-Interim showed the best performance for all variables including specific humidity, the key variable to produce the integrated water vapor.
Thus, we used air temperature (T i ) and specific humidity (q i ) on 37 atmospheric pressure levels from ERA-Interim data to compute the proposed correction.
Recall that the index k refers to the grid point surrounding the GNSS site while the index i refers to the atmospheric pressure level.As we mentioned before, the GNSS geopotential (z GN SS ) is set as a reference, and the value of the geopotential from ERA-Interim (z k ERAInterim ) at each of the 4 grid points surrounding the GNSS site are generally not the same and could differ up to two orders of magnitude.Commonly, neither z GN SS nor the geopotential at any of the 4 grid points matches the geopotential of the nearby pressure level.Therefore, the values of all parameters in the adjacent levels must be used to interpolate (or extrapolate) pressure, temperature and specific humidity in the unknown geopotential (z GN SS and z k ERAInterim ).Thus, the expression of the pressure at an unknown geopotential z j , where j can be any of the unknowns, with respect to a given reference data level (z 0 ) at i = 0 is (van Dam et al., 2010) where T 0 and p 0 refer to the temperature and pressure values at a reference level z 0 , R = 287.04J kg −1 K is the gas constant and λ = 0.006499 K m −1 is the lapse rate of the temperature, and δz is the geopotential difference between z i and the reference level z 0 .Notice that δz is different of ∆z, where the ∆z refers to the difference between z GN SS and z N W M .
The numerator of Eq. ( 5) is the temperature estimated at the desired geopotential z j assuming that the temperature decreases with altitude according to λ.This expression is used to compute p both in z GN SS and in the four grid points of the model (z k ERAInterim ).Finally, the specific humidity (q) is also estimated at the desired z j by a linear interpolation (extrapolation) from data at the adjacent layers.
After knowing p, T and q at each geopotential, z GN SS and the 4 grid points of z k ERAInterim , we can estimate the necessary corrections to the grid points.Such additive corrections to the IWV values at the grid points are equivalent to move the static geopotential of the grid to the z GN SS .Then, the corrected IW V N W M is obtained at the GNSS site by a bi-linear interpolation of the 4 corrected values.
Each value of IWV provided by ERA-Interim is the result of the numerical integration of the expression (Berrisford et al., 2011).
where g 0 is the standard acceleration of the gravity at mean sea level, q(p) is the specific humidity of the air at the pressure level p and the integral is calculated from the first level (p 1 ) up to the model surface level (p s ), i.e. up to the static geopotential that corresponds to the station.
Thus, the proposed correction can generally be written by: where the NWM is ERA Interim, A corresponds to the highest z (z GN SS or z k N W M ) and B the lowest z.Note that the values of q and p in equation ( 7) can be computed as explained before (equation 5).Thus, q j and p j are q and p at z j and q j+1 and p j+1 are q and p at z j+1 , respectively.The values of p grow downwards resulting p 1 = 1 hPa.and p 37 = 1000 hPa.Assuming that the integral of the water vapor is computed from topside and downwards, if the height of a given point from a model is located lower than the position of the receiver, the model integrates a larger column of water vapor and the opposite if the geopotential value from model is larger than the geopotential of the GNSS receiver.Hence, this quantity has to be additive if or subtractive if opposite and the sign is determined by n (see Eq. ( 7)).
In a given instant, we know the geopotential of the GNSS station and the static geopotential assigned by the NWM to the 4 grid points surrounding it (z GN SS and z k N W M , k = 1,2,3,4).We also know the geopotential at 37 pressure levels (z i ) from 1 hPa till 1000 hPa, as well as specific humidity (q) and temperature (T ) at these levels.We should consider that at any time the pressure value of each level is constant but this is not necessarily the case for geopotential height.
Figure 3 illustrates the application of the correction to an example.We take just 1 of the 4 grid points and lets suppose that both unknowns (z k N W M and z GN SS ) are located between the levels 27 (750 hPa) and 28 (775 hPa).Thus, we could use the available data at levels 27 and 28 along with Eq. ( 5) and the before mentioned considerations to estimate p, t and q at z k N W M and z GN SS .Finally, ∆IW V is computed by means of Eq. (7) for this example.

Results
Before analyzing the results of the correction process explained in the previous section, we will present a validation of the numerical integration method used.To this end, we calculate the values IW V ERAInterim at each grid point by using the numerical integral of the equation 6.The integration limits range from 1 hPa to the static geopotential value assigned by the model to the point (z i ERAInterim ).Table 3 shows the results obtained with this procedure.In each grid point the mean value of the differences (IW V ERAInterim data -IW V ERAInterim calculated) is presented.Standard deviations are also shown.It can be seen that in general the resulting values are very close to zero.
In order to evaluate the improvements introduced by the correction, figure 4 shows the ∆IW V as a function of ∆z after applying the proposed integral correction to ERA-Interim data.At first glance we can see that, regardless of whether ∆z is positive or negative, the differences ∆IW V decrease to 2 kg m −2 for all the stations except in 3 cases where they barely exceed that value.Moreover, the correlation between ∆IW V and the geopotential difference decreases to 0.13 as expected.
In addition, if we focus on the plot area of figure 4 limited for ∆z = ± 1500 m 2 s −2 and ∆IW V = 1.5 kg m −2 (zoom not shown), we can see that most of the stations (94 %) are in there.This shows that the proposed correction decreases ∆IW V even for low stations (small z) which generally have the smallest values of ∆z.
The performance of the proposed correction can also be seen in Figure 5.The plots are arranged in two columns where the left column shows stations with positive ∆z, it means that GNSS station is higher to the location assigned by ERA-Interim.
Accordingly, the model integrates a thicker layer of atmosphere and thus IW V ERAInterim values resulted larger than ones from IW V GN SS .The opposite (∆z is negative) is represented by the sites at the right column.Moreover, the differences in ∆z are presented increasing from top to bottom in each column.
We can see that the most important corrections are at BOGT in Bogotá, Colombia, and SANT in Santiago de Chile, Chile.In these examples the differences (IW V GN SS − IW V ERAInterim ), which can reach up to 7 kg m −2 , are significantly reduced.
However the application of this correction, in some cases, should be precautionary.Effectively, sometimes different shortcomings of the model overlap the height problem and therefore the proposed correction could not work.As an example of this we can mention the case of coastal and/or insular stations where 2 or more grid points will be in the ocean.In all these cases the value of IWV calculated from the bi-linear interpolation will be overvalued.Let's analyze in detail the case of stations near the seashore (for example PARC in Punta Arenas, Chile) where 2 of the 4 grid points are in the ocean (see Figure 6).Also ∆z = -1271.86 in PARC indicating that the geopotential from ERA Interim is larger than the GNSS geopotential and therefore the proposed correction will be additive.Besides this result, the IW V ERA−Interim resulted over-estimated by applying a bi-linear interpolation that uses data points in the ocean.In conclusion, the value (IW V ERA−Interim + correction) will overestimate IW V GN SS .Thus, this is an example where applying the suggested correction may worsen the results.The same situation is presented in RIO2 at the Argentinean Atlantic coast.

Discussion and Conclusions
The effect of different heights when comparing results from several data sources not only affects the determination of IWV but also other parameters.For instance, Gao et al. (2012) studied the height corrections for the ERA-Interim 2m-temperature data at the Central Alps and they also found large biases that must be corrected in mountainous areas.Some other authors, also studied the tropospheric refraction effects on space geodetic techniques by considering this effect.For example, Teke et al.
(2013) performed an inter-technique comparison of ZTD in the framework of 4 continuous Very Long Baseline Interferometry (VLBI) campaigns also including NWM and taking into account the effect of the height differences.
The NWM users commonly utilize the IWV values on a grid and calculate with them the IWV value at the desired place by using some interpolation method.In this work, taking the values of IW V GN SS as reference, we show that there are cases where the IWV values obtained from the NWM have differences of several kg m −2 and these discrepancies are mainly due to the difference in geopotentials.
We analyzed the discrepancies between the vertically Integrated Water Vapor values provided by two re-analysis models (ERA-Interim and MERRA-2) with respect to the IW V GN SS values taken as a reference in the South and Central American continent for the period 2007-2013.The results of this comparison allow us to ensure that MERRA-2 resulted wetter than GNSS while ERA Interim is slightly dryer.In addition, when geopotential differences are moderate or large (|∆z| > 500 m 2 s −2 ) and IW V GN SS > 20 kg m −2 , the discrepancies (∆IW V ) are greater than 2 kg m −2 at about 22% of the stations for both models.
Several authors reported problems related to the elevation correction for data from the reanalysis models.The artificial bias in IWV introduced by this altitude difference was previously reported by Bock et al. (2007); Heise et al. (2009); Van Malderen et al. (2014); Bordi et al. (2014) and Bianchi et al. (2016a).Heise et al. (2009) derived IWV from global GNSS ZTD at almost 300 sites by using ground pressure and temperature values from ECMWF and then compared IWV from GNSS and ECMWF.Similar to our work, they found large discrepancies in mountain regions due to the difference in altitudes that caused errors in the meteorological values estimations.Moreover, the analysis performed in our work is also in agreement with Bordi et al. (2014) who compared between GNSS and ERA-Interim IWV values but on monthly time scale for the period 2002-2012.In this case the authors found significant biases of 6.4 kg m −2 in BOGT and 2.5 kg m −2 in BRAZ and they related them to the different elevations between the correspondent GNSS site and the grid points of the model.Thereby, the values of Bordi et al. (2014) are comparable with the correspondent ∆IW V values estimated by our study (6.87 kg m −2 in BOGT and 2.03 kg m −2 in BRAZ, see Table 2).
In this work we proposed an integral correction that compensates on IWV the effect of the geopotential difference between GNSS and the interpolated grid points in the reanalysis model.The results were tested with the respective ones from ERA-Interim.The correction is computed as the numerical integration of the specific humidity where the integral limit is a pressure difference at ∆z (see Eqs. 7).Taking into account that prior to the correction, the 67% of the stations have ∆IW V < 1.5 kg m −2 for ERA-Interim, the application of the numerical correction improves the results and thus the percentage of stations below 1.5 kg m −2 increased to 94%.
Nevertheless, the application of this correction is not advisable at coasts and insular stations of South and Central America because the overvaluation of the model near the coast is overlapped to the height problem.These results are in agreement with Ning et al. (2013), who also compared IWV from GNSS with values from two NWM (ERA-Interim and Rossby Centre Atmospheric climate model, RCA) in Europe for 14 years.The authors also found that models give IWV values larger than GNSS at the seasides or coasts where the tile of the model includes more than 60% of water.
For this reason, the corrections we propose are always recommended but they are not advisable at sea coastal areas or on islands since at least two grid points of the model are usually in the water.
This research was supported by the National Scientific and Technical Council of Argentina (CONICET) PIP 112-201201-00292 and La Plata National University (UNLP) project 11G/142.We would also like to thank the people, organizations and agencies responsible to collect, compute, maintain and openly provide the observations and the products employed in this work: The European Centre for Medium-Range Weather Forecasts (ECMWF) for providing the ERA-Interim reanalysis data (http://apps.ecmwf.int/datasets/).and the Global Modeling and Assimilation Office (GMAO) from National Aeronautics and Space Administration (NASA, USA) for providing MERRA-2 data Table 1: Geopotential values at the selected GNSS stations.Values of z ERA−Interim and z M ERRA−2 come from a bi-linear interpolation of the 4 gridded values of z around the GNSS site.      1 and 2).5), whereas the values of q(z k ERA−Interim ) and q(zGNSS) resulted from a linear interpolation.The computation of ∆IW V k must be done four times (k = 1,...,4), prior to the bi-linear interpolation that produces ∆IW VERA−Interim for the given instant.
where NWM indicates ERA Interim and MERRA-2.All the averages are computed over the period 2007-2013.

Figure 2 .
Figure 2. (up) Values of ∆IW V (∆IW V = IW V GN SS − IW V N W M ) as a function of the geopotential differences ∆z.Results for MERRA-2 are on the left and the same for ERA Interim on the right.(down) Geographical distribution of ∆IW V for MERRA-2 (left) and ERA Interim (right).

Figure 3 .
Figure 3. Scheme of the applied correction to the IWV from ERA-Interim reanalysis for one of the four grid points (k) at a given instant.Both unknowns (zGNSS, dark grey dashed line and z k N W M , thick green line) are located between the pressure levels 27 (750 hPa) and 28 (775 hPa) indicated with thick dashed lines.Atmospheric pressure (pi), temperature (Ti), specific humidity (qi) and geopotential (zi) are known in the 37 levels from 1 hPa till 1000 hPa.Values of p(z k ERA−Interim ) and p(zGNSS) are calculated from equation (5), whereas the values of q(z k ERA−Interim ) and q(zGNSS) resulted from a linear interpolation.The computation of ∆IW V k must be done four times (k = 1,...,4), prior to the bi-linear interpolation that produces ∆IW VERA−Interim for the given instant.

Figure 4 .
Figure 4. Values of ∆IW V (∆IW V = IW V GN SS − IW V ERAInterim) as a function of ∆z after applying the proposed correction.

Figure 5 .Figure 6 .
Figure 5. Residuals of the difference (IW VGNSS − IW VERAInterim) (GNSS -ERA I, black line) along with residuals of the difference [IW VGNSS − (IW VERAInterim + correction)] (GNSS -ERA IC, blue line).Left column shows stations with positive ∆z, it means that GNSS station is higher to the location assigned by ERA-Interim, and the opposite is at the right column.Mean values of the residuals along with the standard deviations are also provided.The sites are shown according to ∆z decreasing from top to bottom at each column.

Table 1 :
Geopotential values at the selected GNSS stations.Values of z ERA−Interim and z M ERRA−2 come from a bi-linear interpolation of the 4 gridded values of z around the GNSS site.

Table 2 :
Differences of the mean values of IWV (∆IW V in [kg m −2 ]) between GNSS and the NWM for the period 2007-2013 at 67 stations located in South America and Central America.The mean value (IW V ) from GNSS at each site is also given and SD refers to the standard deviation.∆z = z GN SS − z N W M refers to the difference in the geopotential [m 2 s −2 ] at each

Table 2 :
Differences of the mean values of IWV (∆IW V in [kg m −2 ]) between GNSS and the NWM for the period 2007-2013 at 67 stations located in South America and Central America.The mean value (IW V ) from GNSS at each site is also given and SD refers to the standard deviation.∆z = z GN SS − z N W M refers to the difference in the geopotential [m 2 s −2 ] at each

Table 2 :
Differences of the mean values of IWV (∆IW V in [kg m −2 ]) between GNSS and the NWM for the period 2007-2013 at 67 stations located in South America and Central America.The mean value (IW V ) from GNSS at each site is also given and SD refers to the standard deviation.∆z = z GN SS − z N W M refers to the difference in the geopotential [m 2 s −2 ] at each

Table 3 :
Mean values of the difference between IW V ERAInterim data and IW V computed from the numerical integral of the equation 7 at each grid point surrounding the GNSS site.The integration limits range from 1 hPa to the static geopotential value assigned by ERA Interim to the point (z k ERA Interim ).

Table 3 :
Mean values of the difference between IW V ERAInterim data and IW V computed from the numerical integral of the equation 7 at each grid point surrounding the GNSS site.The integration limits range from 1 hPa to the static geopotential value assigned by ERA Interim to the point (z k ERA Interim ).

Table 3 :
Mean values of the difference between IW V ERAInterim data and IW V computed from the numerical integral of the equation 7 at each grid point surrounding the GNSS site.The integration limits range from 1 hPa to the static geopotential value assigned by ERA Interim to the point (z k ERA Interim ).