A semiempirical error estimation technique for PWV derived from atmospheric radiosonde data

A semiempirical method for estimating the error and optimum number of sampled levels in precipitable water vapour (PWV) determinations from atmospheric radiosoundings is proposed. Two terms have been considered: the uncertainties in the measurements and the sampling error. Also, the uncertainty has been separated in the variance and covariance components. The sampling and covariance components have been modelled from an empirical dataset of 205 high-vertical-resolution radiosounding profiles, equipped with Vaisala RS80 and RS92 sondes at four different locations: Güímar (GUI) in Tenerife, at sea level, and the astronomical observatory at Roque de los Muchachos (ORM, 2300 m a.s.l.) on La Palma (both on the Canary Islands, Spain), Lindenberg (LIN) in continental Germany, and Ny-Ålesund (NYA) in the Svalbard Islands, within the Arctic Circle. The balloons at the ORM were launched during intensive and unique site-testing runs carried out in 1990 and 1995, while the data for the other sites were obtained from radiosounding stations operating for a period of 1 year (2013–2014). The PWV values ranged between ∼ 0.9 and ∼ 41 mm. The method sub-samples the profile for error minimization. The result is the minimum error and the optimum number of levels. The results obtained in the four sites studied showed that the ORM is the driest of the four locations and the one with the fastest vertical decay of PWV. The exponential autocorrelation pressure lags ranged from 175 hPa (ORM) to 500 hPa (LIN). The results show a coherent behaviour with no biases as a function of the profile. The final error is roughly proportional to PWV whereas the optimum number of levels (N0) is the reverse. The value of N0 is less than 400 for 77 % of the profiles and the absolute errors are always< 0.6 mm. The median relative error is 2.0 ± 0.7% and the 90th percentile P90 = 4.6 %. Therefore, whereas a radiosounding samples at least N0 uniform vertical levels, depending on the water vapour content and distribution of the atmosphere, the error in the PWV estimate is likely to stay below ≈ 3 %, even for dry conditions.


Introduction
Although the water vapour (WV) accounts for only 0-4 % of all atmospheric molecules, it is a powerful greenhouse gas, with strong lines of absorption and emission in the infrared (IR).Atmospheric WV also participates in processes affecting the global climate (Elliott and Gaffen, 1995;Ahrens, 2003) and is the principal molecule responsible for atmospheric extinction in IR astronomical observations, especially at wavelengths longer than ∼ 15 microns (far IR), in several bands in the middle and near IR, as well as in the submillimetre and microwave range (Selby and Mampaso, 1991;Hammersley, 1998;García-Lorenzo et al., 2010;Otárola et al., 2010).
The total amount of WV above a particular location is highly variable and can be expressed as the precipitable water vapour (PWV), which is defined as the total water column height when integrated from the surface to the top of the atmosphere with unit cross section1 .PWV is commonly expressed in millimetres, in terms of the height to which that water would stand if completely condensed and collected in a vessel with a cross section of 1 m 2 .
PWV can be measured by radiosounding balloons, radiometers from both ground (Fowle, 1912;Guiraud et al., 1979;Carilli and Holdaway, 1999;Smith et al., 2001) and satellites (Grody et al., 1980;Menzel et al., 1998;Gao and Kaufman, 2003;Deeter, 2007;Wong et al., 2015), sun photometers (Bird and Hulstrom, 1982;Volz, 1983;Plana-Fattori et al., 1998;Firsov et al., 2013), lunar photometers (Barreto et al., 2013), GPS receivers (Bevis et al., 1992(Bevis et al., , 1994)), Fourier transform infrared spectrometers (Kurylo, 1991;Schneider et al., 2006) and others (Schneider et al., 2010).Of these techniques, atmospheric radiosoundings are a direct in situ measurement and one of the most accurate methods of retrieving the PWV.Radiosoundings are also one of the current standards for atmospheric research and are widely used as a valid reference for comparison and calibration.Therefore, a proper error estimation of PWV obtained from radiosoundings is extremely important.Previous works have dealt with the issue of the accuracy of radiosonde measurements of PWV, including experimental errors and differences in sensor output because of variations in wetness in the column or sensor types (see, for example, Miloshevich et al., 2006, or Romero-Campos et al., 2011).Nevertheless, as radiosondes provide a discrete profiling of the atmospheric column, other important factors impacting the final error in PWV measurements, such as the propagation of uncertainties and vertical sampling, have to be taken into account (Liu et al., 2000).These latter error sources have not yet been explored in depth and are the subject of this paper.

Rationale and objectives
Accurate error estimation of PWV from radiosoundings is essential for regression analyses in comparison or calibration studies.In this sense, special care is needed when radiosoundings from different sources and with differing characteristics are being discussed, or when working in a particularly dry atmosphere.This was the case in Pérez-Jordán et al. (2015), where we used a set of 23 high-vertical-resolution balloons to validate the Weather Research and Forecasting (WRF) mesoscale numerical model for a dry astronomical location (Roque de los Muchachos Observatory, Canary Islands, Spain; ORM hereafter).The possibility of making use of these radiosounding datasets provided a unique opportunity, as no other atmospheric balloons had been launched at the observatory.
A second sample of 42 radiosoundings were also included in Pérez-Jordán et al. (2015) for verification and control.The point selected was the radiosounding operational station located close to sea level in Güímar (GUI, hereafter), in the neighbouring island of Tenerife.The data were downloaded from the repository maintained by the Department of Atmospheric Science of the University of Wyoming2 (WYO hereafter).Both the ORM and GUI-WYO datasets showed strong differences in sampling (∼ 100 levels at GUI-WYO versus more than 2500 at ORM) as a consequence of the reencoding applied at WYO.
In the present work, we are using the same dataset of highvertical-resolution radiosounding flights at ORM and GUI to model the error in PWV calculations from balloons and we propose a method to retrieve it.The dataset at GUI has been extended to 62 profiles.In this study we obtained the data from the Spanish State Meteorological Agency (AEMet) to include all the available operational data points (> 2500).To extend the validity of the work, we have included two equivalent datasets from two locations: Lindenberg (LIN), in Germany, and Ny-Ålesund (NYA) in the Svalbard Islands, Norway.See Sect.1.2 for a detailed description of all the locations and databases.
The method considers two components in the error estimation for PWV obtained from radiosoundings: the propagated uncertainties in the measures and the sampling error.We estimated the uncertainty contribution by means of analytical error propagation through all the levels sampled by the balloon.The propagation includes the covariance with the incoming pressure levels, which was empirically modelled from the autocorrelation function for each location.The sampling component has also been modelled from the highresolution datasets.
Since the method is based on models empirically fitted to local data, we obtained parameters that are valid for the four sites considered.Nevertheless, the proposed method is general and may be extended to any other location by introducing its local data into the models.The results cover a wide variety of conditions originating from the different path lengths (departures from 5 to ∼ 2200 m a.s.l.), PWV content (∼ 1-40 mm) and troposphere conditions (latitudes 28, 52 and 79 • N).The result also includes an estimate of the optimum number of sampled levels needed to fit into the minimum error.
To validate the method for different samplings, the standard (SD) and significant (SIG) levels have been extracted from a subset of 42 radiosoundings from GUI.The SD levels are the basic pressure levels for which the radiosounding data are submitted.They correspond to 1000,925,850,700,500,400,300,250,200,150,100,70,50,30,20 and 10 hPa, although the last ones are not always reported, depending on the top height reached by the balloon (≈ 15 levels).The SIG codification provides ≈ 30 additional measurements to ensure that the linear interpolation between adjacent levels is not missing significant changes in the profiles.The coding rules allow linear interpolations between SIG points to differ from the original measurements by up to ±1 K for temperature and up to ±15 % for relative humidity in the troposphere (up to ∼ 300 hPa) and up to ±2 K for temperature in the stratosphere.Finally, the same profiles, downloaded from WYO, which are encoded to ∼ 100 significant levels, have also been included in the validation.The method was applied to the most pessimistic case, considering only the SD levels, to the SIG levels alone, the SD + SIG and the WYO re-encoded sub-samplings.
The ORM is one of the most competitive astronomical observatories on Earth (see Vernin et al., 2011 for a review); it hosts, among others, the largest optical-IR telescope in the world, the 10.4 m Gran Telescopio Canarias (GTC).The altitude at the ORM ranges from ∼ 2200 to ∼ 2400 m a.s.l., balloons being launched from the lowest level.
The ORM data come from a unique dataset obtained during intensive site-testing assessments performed in April 1990, July 1990 and November 1995 in a joint run carried out by the Instituto de Astrofísica de Canarias and the University of Nice with support of the Centre National de Recherches Météorologiques (Météo-France).The details of these campaigns are described in Vernin andMuñoz-Tuñón (1992, 1994).A total of 23 balloons were launched from the observatory, 22 of them being selected to model the error, after rejecting the one that did not cover the whole troposphere.The profiles provide a dataset of very high vertical resolution with 2500 to 7000 data points per balloon flight.
Because of their latitude and eastern location in the North Atlantic Ocean, the Canary Islands exhibit a vertical troposphere structure with a trade wind thermal inversion layer (IL), driven by subsiding cool air from the descending branch of the Hadley cell.The altitude of the IL ranges on average from 800 m in summer to 1600 m in winter, well below the altitude of the ORM (Dorta, 1996;Carrillo et al., 2016).The IL separates the moist marine boundary layer from the dry free atmosphere, inducing very high atmospheric stability above it.Therefore, the PWV values at both locations in the Canary Islands (GUI and ORM) are strongly uncorrelated with wide differences in their atmospheric conditions.The average mesoscale conditions in the Canary Islands show a relative humidity profile under the IL that drops drastically with increasing altitude.Temperatures below −40 • C are typically reached above 8 km in the troposphere.
GUI is the closest radiosonde operational station to ORM, located ≈ 155 km away and closest to sea level.It is an AEMet station, part of the World Meteorological Organization (WMO) network, code 600183 .The sample from GUI covers 1 year of data with a total of 62 operational radiosoundings, mainly from May 2013 to April 2014.The balloons at GUI are routinely launched twice daily (at 00:00 and 12:00 UTC) by the AEMet.The significant level assignment for GUI, used for validation in Sect.5.4, was carried out by the AEMet following WMO standards (see Sect. 1.1).
The LIN and NYA observatories are certified stations of the GCOS (Global Climate Observing System) Reference Upper-Air Network (GRUAN) 4 and are also part of the WMO network.LIN is a climate reference site of the Deutscher Wetterdienst (DWD) and hosts the GRUAN Lead Centre.The balloons are routinely launched every 6 h.LIN is located in an all-year-wet, maritime temperate climate 5 .NYA is an atmospheric research observatory operated by the Alfred-Wegener-Institut für Polarforschung (AW) and the Helmholtz Centre Potsdam German Research Centre for Geosciences (GFZ).The balloons are launched at an approx-imate rate of one per day.NYA is in a climatic-strategic location at the gateway to the Arctic Ocean.The data compiled from LIN and NYA are statistically equivalent to the dataset from GUI, with 58 and 63 soundings, ranging from May 2013 to April 2014.

Instrumental set-up
The radiosondes used were a Vaisala RS80 at ORM (1990 and1995) and a Vaisala RS92 for the other operational soundings (2013)(2014).In both cases, the sondes were calibrated during the ground check immediately before the balloon release.The standard Vaisala procedure for sonde preparation and systematic error source minimization were followed, together with various corrections for daytime solar heating.For an extensive characterization of both radiosonde set-ups see Miloshevich et al. (2001Miloshevich et al. ( , 2009)).In the particular case of GUI, more details of the corrections applied are given in Romero-Campos et al. (2011).The quality controls for GRUAN (LIN and NYA) are detailed in Immler et al. (2010) and Dirksen et al. (2014).
The uncertainties assumed in this work for both radiosondes are ±0.5 • C, ±5 % and ±1 hPa for temperature, relative humidity and atmospheric pressure.These values are in agreement with the estimates for the RS92 obtained within the extensive network of GRUAN by Dirksen et al. (2014).The equations accept future improvements on sensor uncertainties, such as the inclusion of vertically resolved profiles of uncertainties (Dirksen et al., 2014) that the new binary encodings for radiosonde data, such as BUFR (Binary Universal Form for the Representation of meteorological data), are introducing (Dragosavac, 2007).
Both RS80 and RS92 suffer from time lag in the relative humidity sensor response that increases as temperature drops.The delay is caused by the slow response of the sensor (a film capacitor that acts as a dielectric) at low temperatures.As noted by Miloshevich et al. (2006), this time lag error is not a bias error but rather a smoothing of the profile by an amount that depends on the temperature and local humidity gradient.Such a time lag was studied and characterized by the empirical approximating expression of Miloshevich et al. (2001).The response time increases significantly in the upper troposphere for temp < −68 • C (Vömel et al., 2007), which is reached typically above 10 000 m, where the water vapour mixing ratio drastically drops, reducing the contribution to the total integrated PWV.Besides we can assume a uniform distribution in water vapour in the upper troposphere for the time delay.Under these conditions, the time lag has the effect of smoothing the signal and reducing the dispersion, without significantly moving the average and therefore not impacting on the total PWV in the column.

PWV from radiosondes
The PWV is obtained from the temperature T ( • C), atmospheric pressure p (hPa) and relative humidity RH (%) measured by the radiosondes (Curry and Webster, 1999).Following the definition given in Sect. 1, the PWV can be expressed by where z is the height in m, and ρ and ρ w are the liquid water and WV densities, both in kg m −3 .The definition of the WV mixing ratio r is where m w and ρ w are the WV mass and density, and m d and ρ d are the corresponding values for dry air.We can assume hydrostatic balance (dp = −ρ d • g • dz) and write Eq. ( 1) in the form where g is the Earth's gravity (m s −2 ) and p s and p t are the pressure at the surface and at the top of the sampled atmospheric column in hPa.We can now apply the ideal gas law and Dalton's law of partial pressures to Eq. (2) to get The coefficient 0.622 is the molecular mass ratio of WV in dry air and e is the partial vapour pressure that can be obtained from the definition of relative humidity as e = e sat • RH 100 (hPa). (5) Finally, the saturation vapour pressure e sat can be expressed as an empirical polynomial fit, following Curry and Webster (1999): The coefficients a i have been taken from Flatau et al. (1992) (see Table 1).
Table 1.Saturation vapour pressure coefficients and associated uncertainties (Flatau et al., 1992).95246610 × 10 −4  a 3 2.65027242 × 10 −4  1.20785253 × 10 −5  a 4 3.02246994 × 10 −6  1.01581498 × 10 −7  a 5 2.03886313 × 10 −8 3.84142063 × 10 −10 a 6 6.38780966 × 10 −11 6.69517837 × 10 −14 3 Error budget for PWV from radiosonde data The integral in Eq. ( 3) is computed as a discrete summation (following a trapezoidal method) over all the levels sampled by the radiosonde.We can write the final error associated with the PWV determination, f , as The last term ( s ) is the sampling error and σ is the uncertainty in the measurement, which can be written as where var f is the final variance, expressed as the sum of the variances at each level, σ 2 i , and cov f is the covariance as a function of the autocorrelation between levels ρ ij .

Uncertainty propagation: variance component
Let σ p,i , σ T ,i and σ RH,i be the instrumental uncertainties associated with the direct measure of the atmospheric pressure p, temperature T and relative humidity RH at a certain level i.The variance of PWV can be obtained by error propagation over all the sampled levels N in Eq. ( 3).Nevertheless, each value, excluding the extremes, is shared by adjacent bins; therefore, the derivation of the variance must avoid the double summation (see details in Appendix A).Assuming p and r to be independent, we can write the variance of the trapezoidal expansion of Eq. (3) as in Eq. (A3): The variances of the water vapour mixing ratio, r (σ 2 r ), and pressure, p (σ 2 p ), may be obtained from the Eqs.(A7) and (A9): 11) where p i and r i are the pressure and mixing ratio at level i.
By a recursive use of the error propagation rules in Eqs. ( 4), ( 5) and ( 6), the following expressions are obtained for the uncertainties σ r,i , σ e,i and σ e sat ,i : where a j and σ a,j are the seven (j = 0, 6) saturation vapour pressure coefficients and the associated uncertainties for each sampled level i (see Table 1).The contribution of the covariances between variables in the propagation of errors was specifically calculated, leading to negligible values less than 10 −3 mm, and is therefore not included in Eqs. ( 9)-( 14).

Uncertainty propagation: covariance component
We may consider the vertical profiles as time series generated by the balloon flights (space series).Each observation of the sensors on board the balloon may be viewed as the only possible realization of the physical process governing the atmosphere at this point and time.Therefore, each measurement at a particular height during the flight will be a random value within the accuracy of the sensor that can be assumed as a semi-Markov process, where the influence on the error of the next levels may be modelled through the covariance (see Eq. 8) as a simplified pure exponential autocorrelation function: where the exponential autocorrelation pressure lag (τ ρ ) was empirically determined for each location after obtaining the www.atmos-meas-tech.net/9/4759/2016/Atmos.Meas.Tech., 9, 4759-4781, 2016 autocorrelations ρ ij between a sub-sample of r i and r j , separated by the lag p i − p j , from all the available radiosoundings.The lag was recursively increased while ρ ij ≥ 0.2.For a fixed time lag, the departure levels i were sequentially shifted up to a top limit such as level + lag, approximately reaching the tropospheric height where the content of water vapour drops.Finally the average autocorrelation decay curve was fitted to an exponential function exp[−(p i − p j )/τ ρ ].In Sect.5.1 we show the results for each location.

Sampling error
The number of levels included for PWV determinations from atmospheric radiosondes may range from tens (when only standard levels are available) to thousands (full profile).All the high-resolution profiles compiled in this study range from ∼ 2500 to ∼ 7000 data points, dense enough to neglect the sampling error.This circumstance allowed us to empirically obtain an expression for s as a function of the number of sampled levels N, following a recursive sub-sampling process.
Each profile was uniformly sub-sampled at equal intervals by taking one point in two, one in three, etc., to obtain 800 different realizations of PWV for the same profile (any other uniform sub-sampling is valid).The dispersion of the residuals increases logarithmically as the number of levels N decreases (see Fig. 2a for an example at GUI), with the residuals defined as where I and I N are the integral in Eq. ( 3) calculated with all the levels in the profile (N max ) and with the different realizations of N levels after the sub-sampling process.The residuals were grouped in slices for intervals of N to fit a model.The size of the slices was selected following a quasilogarithmic scale to overcome the differences in variance (heteroscedasticity) while conserving the statistical significance (see Fig. 2a for details of the slices and the number of residuals included for each).The sampling error was then obtained as the root-mean-square error (RMSE) of the residuals for each slice.The RMSE was calculated as the square root of the sum of the variance and the squared bias for every interval.
The dependence on N of res N comes from the integrals I N in Eq. ( 16) in the form where the first term, PWV, is the composite trapezoidal sum of Eq. ( 3) sub-sampled to N levels and the second one is the associated error.Therefore, taking N as the middle point of each slice interval in Fig. 2a, we modelled the residual (Eq.16) behaviour by fitting a function A/N + B/N 2 + C to the RMSEs, with the coefficient C = 0, as lim (N→∞) res N = 0. Finally, we applied a gradient-expansion algorithm to compute a nonlinear least squares fit to obtain the coefficients s 0 and s 1 (red line in Fig. 2).
In Sect.5.2 we show the results for each location.Equation (18) assumes that the radiosounding is uniformly sampling the whole atmosphere where the PWV concentrates (mainly the lower and mid-troposphere).The limit depends on each particular site but it may be established in ≈ 20 km.All the radiosoundings in the present study fulfil this condition except for ORM-VOL18, which reached 17 797 m and was included after a visual inspection (see Appendix B for a compilation of all the data).The bias induced by the shortage in height profiling is beyond the scope of this paper, as it would need a detailed model (with larger samples) of the probability of missing high-altitude wet layers, depending on the local troposphere and season.These models may also allow for alternative nonuniform sub-samplings, concentrating the measured levels where a priori a less uniform PWV distribution is to be expected (i.e. that would most impact on the sampling error).

Optimized error
The two components in Eq. ( 7) behave contrary to the number of sampled levels.Because of the exponential decay in the autocorrelation, the closer the levels, the larger the covariance contribution to the final error and the less the information obtained to the final integrated PWV.Therefore, σ increases with the number of levels in the profile.In contrast, the sampling error is in the opposite direction and s increases as the number of sampled levels decreases.This behaviour is described in Fig. 2b, which shows the sampling error s fitted in Eq. ( 18) and the median σ calculated as a function of the sampling levels N for each slice of subsampled data.Therefore, it is always possible to sub-sample the profile uniformly, obtaining f (N, PWV) by use of Eq. ( 7) for different values of N (N ≤ N max ).Hence, the optimized error will result in minimizing f whilst reducing N (see Fig. 2b).
where f is the final error defined in Eq. ( 7), σ is the propagated uncertainty (Eq.8) and s is the sampling error (Eq.18).The optimum number of samples N 0 will be defined as the argument of the minimum in Eq. ( 19).
Finally, we can calculate the individual contribution of the uncertainty and the sampling error to the optimized error by means of 5 Results In this section we apply the method described above to the radiosounding dataset from ORM, GUI, LIN and NYA.In Table 2 we present a summary of the main statistical results for each location.The four sites show different behaviours in both the atmospheric PWV distribution and concentration.ORM and NYA show low contents of PWV, ORM being the driest site.Although the results cannot be taken as a climatic result, owing to the limitations in the sampling periods, they highlight the qualitative differences between the sites and help to discuss the results obtained in the error analyses.
We also discuss the empirical estimation of the parameters needed in the models used by the method at each site.All the numerical results are listed in the Appendix B.

Exponential autocorrelation lags results
Figure 3 shows the results of τ ρ estimated from the autocorrelation decay curves obtained from the radiosounding data at each location, as explained in Sect.3.2.The empirical average curves show good agreement with the exponential fits (red and blue lines, respectively, in Fig. 3).The values of τ ρ are also shown in the figure.The results are mainly related to the average PWV vertical distribution over each site.In the Canary Islands, this distribution is governed by the presence of the IL (see Sect. 1.2), which induces an abrupt decay in the WV that explains the lower values of τ ρ and the maximum lag where the autocorrelation falls below 0.2 for GUI and The background grey lines are the autocorrelations ρ ij between a sub-sample of r i from all the available radiosoundings and the following r j , separated by the lag P = p i − p j .For each fixed lag, the departure levels i were sequentially shifted up to ρ ij = 0.2, with a maximum level such as the level + lag approximately reach the troposphere height where the content of water vapour drops.
ORM.In particular, for the ORM this result is in agreement with the excellent conditions for IR astronomy (Hammersley, 1998) and with the low altitude reported for the tropopause by García-Lorenzo et al. (2004).

Sampling errors results
Figure 4 shows the best fits achieved to the residuals of the sub-sampling process applied at each location (see Sect. 3.3).
The coefficients obtained are also shown in the figure.The first coefficient s 0 governs the rapid expansion of the sampling errors for sample sizes N 0 .This coefficient is proportional to the average PWV (see Table 2); therefore, the more humid the atmosphere, the faster the sample error grows.The results are in concordance with the median PWV shown in Table 2, the driest sites being ORM and NYA, the ones with lower s 0 .The coefficient s 1 governs the lower part of the curves, for sample sizes N 0 .It may be associated with uniformity in the vertical distribution of PWV.For an inhomogeneous vertical PWV distribution the sampling error is more sensitive to the shortage of measures when the number of levels is below the minimum N 0 .The results show the lowest value for LIN and the largest for GUI, where the IL tend to break the PWV profile abruptly.Table 2. PWV statistics (mm) for the four sites: percentiles (P xx ), median (med), minimum (min), maximum (max) and dispersion (s).The dispersion has been estimated robustly by means of 1.4826 × MAD, where MAD is the median absolute deviation and 1.4826 is the scale factor between MAD and the standard deviation for perfect Gaussian distributions.

Site
n

Optimized error results
We have applied the optimization described in Sect. 4 to all the available radiosondes and locations.For each profile, we have computed Eq. ( 7) recursively while sub-sampling the profile and carrying out the minimization in Eq. ( 19) with the local values of τ ρ , s 0 and s 1 (see Figs. 3 and 4).sites.The first two panels (Fig. 5a and b) show two results for the ORM, while the panels c and d are results for GUI.In both cases, we have selected two profiles with strong differences in the PWV content for each location.Additionally, in Fig. 5e and f, we show the results of the method applied to the same profile in Fig. 5d, but now with the re-sampled data of WYO (N max ≈ N 0 ) and the SD + SIG levels (N max < N 0 ; see Sect.1.1).In all cases, the residuals oscillate inside the range defined by the sampling error model (red line in Fig. 5), while converging.For N max N 0 , the error is dominated by σ , whereas for N max ≈ N 0 the sample error s becomes significant.We sub-sampled the data in the profiles WYO and SD + SIG by extracting N , N −2, N −4, . . .points, uniformly distributed for each iteration, to obtain a sufficient number of realizations close to N max .
Figure 6 shows all the results.The values of and N 0 are plotted together for all the available radiosondes as a function of PWV.The optimized error is roughly proportional to PWV, whereas N 0 is the reverse.The optimum sampling number N 0 is less than 400 for 77 % of the radiosoundings (see Fig. 6 and Appendix B for the numerical values).As we are using two different sonde models (RS80 and RS92), we have indicated the range of values covered for each.No differences or any bias as a function of the profile were observed.
We obtained errors less than ∼ 0.6 mm for all cases; this is significantly lower than the 1.3 mm error published by Liu et al. (2000), also taking into account the sampling ef-Table 3. Relative error ( rel ) statistics (%): percentiles 10th and 90th (P xx ), median (med), minimum (min), maximum (max) and dispersion (s).The dispersion has been estimated robustly by means of 1.4826 × MAD, where MAD is the median absolute deviation and 1.4826 is the scale factor between MAD and the standard deviation for perfect Gaussian distributions.

Min
P 10 Med − s Med Med + s P 90 Max 0.9 % 1.3 % 1.3 % 2.0 % 2.7 % 4.6 % 67.2 % fect in the error, but with no deeper analyses into the dependence on N .Specifically for the sampling component s 0 , we also obtained a lower average value of 0.11 mm with a maximum of 0.29 mm (except for some outliers found in the  B2).
NYA sample), in comparison to 0.50 and ≈ 1 mm, respectively, obtained by Liu et al. (2000) after analysing the residuals between the smoothed data (the standard output) and the much denser real-time records from 50 PWV radiosoundings in Hong Kong.Figure 7 and Table 3 summarize the statistical results for the relative errors ( rel hereafter) for the whole dataset of 205 radiosoundings.Relative errors behave in a slightly opposite way to the absolute optimized errors : the drier the atmosphere, the larger the rel (see Fig. 7a).The median relative error is 2.0 ± 0.7 %, and the 90th percentile P 90 = 4.6 %.Around 10 outlier points have been obtained in the relative errors for NYA, mainly concentrated in the winter period (in the extreme case the relative error reaches 67.2 %).These values are not valid error estimates and define the limits of the model based on empirical data.Visual inspection of the mixing ratio profiles corresponding to those dates shows abrupt irruptions of extremely humid layers in the upper levels.This feature is related to the atmospheric circulation in the extreme conditions during the winter months at the Arctic latitude of NYA and goes far beyond the average conditions for which the parameters τ ρ , s 0 and s 1 were estimated for that location.Therefore, for these points, the simplified pure exponential autocorrelation model fails and a new model, based on knowledge of the specific local atmospheric conditions, must to be implemented.
The complete list of rel is also included in Appendix B. These results reduce by more than a half the uncertainty of ≈ 5 % (≈ 15 % for extremely dry conditions) published by Schneider et al. (2010) following the mixing ratio uncertainties obtained by Miloshevich et al. (2009) from the Vaisala RS92 RH sensors and highlight the importance of optimized sampling in the PWV determinations from radiosoundings.

Validation for poorly sampled radiosonde data (N max < N 0 )
To validate the model, we have compared the errors obtained for the high-resolution radiosoundings of GUI with the same profiles with different re-samplings (WYO, SD + SIG, SIG and SD).The ∼ 15 SD levels are the minima and mandatory levels reported by the radiosondes and may be considered the most pessimistic scenario (N max N 0 ).In this case, the sampling component dominates ( ≈ s ) and the error can be directly estimated from Eq. ( 18) with N = N max .For SIG and SD + SIG the situation is gradually more favourable, but still (N max < N 0 ).Finally, in the WYO data, the sampling approximates to N 0 .The results of the comparison6 are shown in the Fig. 8.In the four cases, the error of the less sampled series includes the PWV and error estimates of the next series up to the high-resolution profiles.Therefore, the modelled errors statistically represent the differences between the calculated PWV and the best available estimate.The relative errors for the GUI-SD series range between ∼ 5 and ∼ 30 %.Therefore, Fig. 8 also evidences the ability of the standard levels to reproduce the tendencies in long-term PWV monitoring programmes.

Conclusions
We have considered a semiempirical approximation to the error propagation for PWV determinations from radiosonde profiles.To do this, we have considered two terms in the error estimation: the uncertainties in the measures σ and the sampling error s (the whole atmospheric column is assumed to be uniformly sampled).Moreover, the uncertainty contribution is separated into variance and covariance components.The variance has been estimated by means of analytical error propagation through all the levels sampled by the balloon, avoiding the correlation between adjacent bins.In contrast, the covariance and the sampling error have been empirically modelled at four locations -in the Canary Islands, Germany and the Svalbard Islands -with different datasets of radiosoundings ranging over 1 year.The covariance model is based on the exponential decay of the autocorrelation function, while the sampling error model is based on the number of samples of the composite trapezoidal formulas for numerical integration and their errors.
The uncertainty σ increases and s decreases as the number of sampled levels grows.Therefore, we have optimized the error by gradually reducing the number of samples N. The optimization leads to the calculation of the minimum error (Eq.19) and may be considered an appropriate estimator of the error in the PWV determination.The value of N in the minimum error is the optimum number of samples N 0 (Eq.20).
The parameters estimated to fit the empirical models are governed by the average PWV content and its distribution.From the values obtained in the four sites studied, the ORM is the driest location with the fastest vertical decay of PWV.
The sub-sampling minimization was applied to all the radiosoundings in the datasets.The results show a coherent behaviour with no differences or bias as a function of the profile.The optimized error is roughly proportional to PWV, whereas N 0 is the reverse.The value of N 0 is less than 400 for 77 % of the profiles and the absolute errors are always < 0.6 mm, with the sampling component s 0 < 0.3 mm.
Two different scenarios arise after the determination of N 0 , whether the actual number of levels in the profile N max is greater or equal than N 0 or not.For N max ≥ N 0 it is always possible to reach N 0 by sub-sampling and therefore to obtain the minimum error.For N max < N 0 the sampling component dominates and the final error can be obtained directly from the model (Eq.18).
The method was validated by the comparison of poorly sampled profiles and high-vertical-resolution data.The results showed that the high-resolution PWV values fall inside the interval defined by the low-resolution PWV values ± the estimated errors.The errors grew by up to > 30 % with poorly sampled profiles for dry atmospheres.The median relative error for high-resolution profiles is 2.0 ± 0.9 %, with 90th percentile P 90 = 4.6 %.These results reduce by more than a half the uncertainties previously reported in the literature.
Therefore, not only the uncertainties define the error in PWV estimations from radiosoundings but also the autocorrelation between levels and the sampling.Here we have proposed that it is possible to optimize the number of sampled levels to minimize the error within the instrumental uncertainty.Whereas a radiosounding samples at least N 0 uniform vertical levels, depending on the WV content and distribution of the atmosphere, the error in the PWV estimate is likely to stay below ≈ 3 % (median + dispersion = 2.7 %) even for dry conditions.

Data availability
The radiosonde balloons were launched at the ORM by the IAC and the University of Nice under the direction of Jean Vernin (jean.vernin@unice.fr) and Casiana Muñoz-Tuñón (cmt@iac.es),with support from the Centre National de Recherches Météorologiques5 (Météo-France).Soundings at Güímar are launched and verified by the Spanish Agencia Estatal de Meteorología (AEMet; http://www.aemet.es).The sounding data at Lindenberg and Ny-Ålesund have been compiled and verified by the GCOS Reference Upper-Air Network (GRUAN; http://www.dwd.de/EN/research/international_programme/gruan/home.html).The low-resolution radiosounding data are available through the page of the Dept. of Atmospheric Science at the University of Wyoming (http://weather.uwyo.edu/upperair/sounding.html).

Figure 1 .
Figure 1.Location and height above the sea level of the four radiosounding launching points in the Canary Islands, Spain (ORM and GUI), separated by ≈ 155 km, Lindenberg in Germany, and Ny-Ålesund, in the Svalbard Islands, Norway (source: Wikimedia Commons).

Figure 2 .
Figure 2. (a) Example of sampling error modelling for PWV radiosondes, as a function of the number of levels (see Fig.4to see the model fitting for all the sites).The profiles are gradually sub-sampled to obtain a large number of realizations of PWV (grey circles).The error is then fitted (red solid lines; see Eq. 18) to the RMSE of the residuals (see Eq. 16) in sliced (dotted lines) sampling levels.The error bars show the RMSE per slice and the "x" symbols are the bias.(b) Sub-sampling error optimization.The solid black line is the experimental sample error as a function of the number of levels and the red line is the model (same as in panel a but with trimmed y axis).The blue line is the median of the propagated uncertainties (Eq.8) for all the profiles and the dashed line is the final error (Eq.7).The arrows show the optimized number of samples for the minimum error.

Figure 3 .
Figure3.Results of the empirical autocorrelation exponential decay model applied of the four sites studied and exponential autocorrelation pressure lag (see Sect. 3.2).The background grey lines are the autocorrelations ρ ij between a sub-sample of r i from all the available radiosoundings and the following r j , separated by the lag P = p i − p j .For each fixed lag, the departure levels i were sequentially shifted up to ρ ij = 0.2, with a maximum level such as the level + lag approximately reach the troposphere height where the content of water vapour drops.

Figure 4 .
Figure 4. Results of the sample error model fit for the four sites studied (Eq.18).The symbols and lines are the same than in Fig. 2a.
Figure 5 shows six particular examples, focused on the Canarian

Figure 5 .Figure 6 .
Figure 5. Six graphic examples with the application of the method of sub-sampling error optimization for ORM (a) and (b) and GUI (c) and (d), with different PWV concentrations and N max N 0 .Panels (e) and (f) show the same profile as (d) but from data from the WYO repository (see Sect. 1.1), where N max ≈ N 0 and with the SD + SIG levels (N max < N 0 ).The solid grey line shows the residual between the sub-sampled PWV estimation and the best value with all the levels available.The vertical dotted lines are the optimized number of samples for the minimum error.The PWV and errors in the legend are in millimetres.

Figure 7 .
Figure 7. Relative error ( rel ) statistics.Panel (a) shows all the PWV values (black), plotted together with their absolute (blue) and relative errors (red) in a logarithm scale.The thin horizontal line and shadow show the median relative error and dispersion.The dispersion has been estimated robustly by means of 1.4826×MAD, where MAD is the median absolute deviation and 1.4826 is the scale factor between MAD and the standard deviation for perfect Gaussian distributions.Panel (b) shows the histogram of the distribution of relative errors.The vertical lines show the median value and dispersion range (dashed).The main statistics is in the legend.

Figure 8 .
Figure8.Comparison between the PWV series and errors (shadow) obtained from the GUI data with high-resolution radiosoundings (> 2500 levels) and different sub-extractions of the same profiles.The colour code, from the less to the most sampled series, is in the legend.The labels in the x axis are the references for each balloon flight (see TableB2).

Table B2 .
PWV data and optimized error (see Eq. 19) from the radiosoundings at GUI. See the caption of TableB1for description.

Table B3 .
PWV data and optimized error (see Eq. 19) from the radiosoundings at LIN. See the caption of TableB1for description.

Table B4 .
PWV data and optimized error (see Eq. 19) from the radiosoundings at NYA. See the caption of TableB1for description.