Interpolation uncertainty of atmospheric temperature radiosoundings

This paper is motivated by the fact that, although temperature readings made by Vaisala RS41 radiosondes at GRUAN sites (www.gruan.org) are given at 1 s resolution, for various reasons, missing data are spread along the atmospheric profile. Such a problem is quite common in radiosonde data and other profile data. Hence, (linear) interpolation is often used to fill the gaps in published data products. In this perspective, the present paper considers interpolation uncertainty. To do this, a statistical approach is introduced giving some understanding of the consequences of substituting missing data by interpolated 5 ones. In particular, a general frame for the computation of interpolation uncertainty based on a Gaussian process (GP) set-up is developed. Using the GP characteristics, a simple formula for computing the linear interpolation standard error is given. Moreover, the GP interpolation is proposed as it provides an alternative interpolation method with its standard error. For the Vaisala RS41, the two approaches are shown to give similar interpolation performances using an extensive cross10 validation approach based on the block-bootstrap technique. Statistical results about interpolation uncertainties at various GRUAN sites and for various missing gap lengths are provided. Since both provide an underestimation of the cross-validation interpolation uncertainty, a bootstrap-based correction formula is proposed. Using the root mean square error, it is found that, for short gaps, with an average length of 5 s, the average uncertainty is smaller than 0.10 K. For larger gaps, it increases up to 0.35 K for an average gap length of 30 s, and up to 0.58 K for a gap of 15 60 s.


Introduction
The quality of climate variable profiles in the upper troposphere -lower stratosphere (UTLS) is relevant in various scientific fields. In particular, it is important for numerical weather prediction, satellite observation validation and climate change understanding, including extreme events such as droughts and tornadoes. 20 The GCOS (Global Climate Observing System) Reference Upper-Air Network (GRUAN, www.gruan.org) is a network for reference measurements of UTLS (Seidel et al., 2009;Bodeker et al., 2016). Immler et al. (2010) discussed the concepts of reference measurements, traceability, full metadata description, a proper manufacturer-independent instrument characterization, and the assessment of measurement uncertainties for upper-air observations. In this frame, GRUAN data processing for the Vaisala RS92 radiosonde was developed to meet the above criteria for refer-25 ence measurements (Dirksen et al., 2014). The related data product is characterised not only by the above mentioned metrological requirements but also by high-vertical-resolution. After the introduction of the new Vaisala RS41 radiosonde, GRUAN is currently developing the corresponding data processing for the new instrument (Dirksen et al., 2019).
Although temperature readings made by Vaisala RS41 radiosonde at GRUAN stations are given at 1 s resolution, for various reasons, missing data are sometimes present along the atmospheric profile. If one is led to interpolate the missing measure-30 ments, since an interpolation error is implied, the related uncertainty is to be considered in the uncertainty budget.
The interpolation of atmospheric profiles has been considered in the literature from various points of view. In some cases, interpolation is applied to the measurement uncertainty. For example, considering the AERONET Version 3 aerosol retrievals, Sinyuk et al. (2020) obtain the uncertainty by interpolation of a look up table.
A second and more relevant use of interpolation is related to the measurement itself. In this field, Ceccherini et al. (2018) 35 used interpolation for data fusion of Ozone satellite vertical profiles. Interpolation uncertainty and more generally co-location uncertainty has been computed using simulated profiles. Similarly, in co-location uncertainty of total ozone, Verhoelst et al. (2015) contemplate interpolation in the so-called OSSSMOSE simulator.
In the frame of radiosonde co-location uncertainty, considering relative humidity, Fassó et al. (2014) used a statistical approach based on the heteroskedastic functional regression model. Considering pressure, Ignaccolo et al. (2015) extended the 40 latter approach to a 3D functional regression approach. In these two papers, the interpolation uncertainty is implicitly assessed by means of the model error variance.
The comparisons of radiosonde and satellite data are often based on low-vertical-resolution radiosonde profiles measurements such as the data collected within the network of the Universal Rawinsonde Observation Program (RAOB) because of their global coverage. In some cases interpolation is not required because of the higher vertical resolution of satellite profiles 45 (Sun et al., 2010). In other cases, interpolation is required. For example, Finazzi et al. (2019) considered the harmonisation of the low-vertical-resolution RAOB temperature and humidity radiosonde measurements and the corresponding atmospheric profiles derived from the Infrared Atmospheric Sounding Interferometer (IASI) aboard Metop-A and Metop-B satellites. In this frame spline interpolation of RAOB profiles was indirectly assessed through a comparison with GRUAN radiosonde reference measurements.

50
As a common trait of the above literature, interpolation of atmospheric profiles is quite common, but a systematic analysis of interpolation uncertainty per se is not yet available. A general approach to interpolation is the Geostatistics approach (Cressie and Wikle, 2011) which is the same as the Gaussian Process (GP) approach (Rasmussen and Williams, 2006) to a large extent.
Its value is due to the fact that it gives optimal interpolation under some conditions. With some variations, the related optimal interpolation algorithm is based on the autocovariance function characterising also the structure function (Sofieva et al., 2008).

55
This approach is often used to interpolate in a higher dimensional space such as the Euclidean plane, the sphere (Alegria et al., 2017), the three-dimensional Euclidean space or the circular shell representing the atmosphere. Interpolating and forecasting is sometimes overlapping, in particular this happens when the GP is defined, for example, on time cross a sphere (Porcu et al., 2016). Interestingly, it can be shown that the spline interpolation is a special case of the GP interpolation (Kimeldorf and It must be noted that the formula is correctly used if the "true data generation mechanism" is a GP. If the GP is simply "an approximation" an additional term must be added.
In this paper, the uncertainty of the one-dimensional linear interpolation is discussed integrating two approaches. In the first stage, the closed form formula of the linear interpolation uncertainty is presented under the assumption that the observed atmospheric profile is generated by a GP. In the second phase, thanks to the availability of appropriate data, the GP assumption is 65 relaxed and a block-bootstrap estimator is constructed. The approach is valid for any atmospheric profile dataset. Considering the motivating application, which focuses on temperature readings of Vaisala RS41 at GRUAN sites, the objective of this paper is to contribute to the understanding of interpolation uncertainty expressed as a function of missing gap length, missing frequency, altitude and site.
To do this, "good" launches without missing data are used. Each profile is divided in a learning set and a testing set, 70 the latter being used as missing data for interpolation uncertainty assessment. This is done for various missing patterns that resemble observed "bad" launches, which are characterised by many missing measurements. In particular, increasing gap average lengths will be analysed. The testing sets will be extracted using a block-bootstrap cross-validation scheme. Hence although the numerical results are specific to Vaisala RS41 temperature data set, the approach is quite general and may be applied to other sensors.

75
The rest of the paper is organised as follows. Section 2 motivates the paper by discussing the soruces of gaps in data reception and their impact in GRUAN data processing. Section 3 introduces the Gaussian Process (GP) set up used to provide the formal assessment of linear interpolation uncertainty and to introduce the GP interpolation with its standard deviation. Section 4 presents the data sets, which are related to Vaisala RS41 observations at seven GRUAN sites and are used in the empirical analysis. Section 5 describes the re-sampling technique able to simulate random patterns of missing values of different duration.

80
Section 6 describes the cross-validation scheme essential for the uncertainty computations and the model selection, which is discussed in Section 7. Section 8 presents the results, compares the behaviour of the two interpolation techniques and proposes an empirically corrected formula for interpolation uncertainty. Eventually, Section 9 draws some conclusions.

Data processing and interpolation
There are several possible reasons for temporary gaps in data reception. These include the presence of obstacles that may inter-85 fere with radio transmission to the ground station (trees, buildings, local geography), extraordinary meteorological conditions, or instrument-related reasons. Considering an ascent as a trajectory rather than a vertical profile, it seems evident that the probability for the occurrence of data gaps tends to increase with the horizontal distance from the launch site (weaker radio signal), which can significantly exceed the vertical distance depending on wind conditions. The GRUAN Lead Centre conducted a statistical analysis for the occurrence of data gaps in RS41 radiosonde soundings performed at 15 GRUAN stations in the period 90 2014-2019. The results show that gaps occur in more than 20% of the soundings, virtually independent of the height ranges, with the majority (>95%) having less than 15 gaps per 1000 s (=1000 data points). Up to 30 km, gaps >10 s only play a role in about 5% of the ascents, however the occurrence of larger gaps generally increases with height (distance). Figure 1 gives an example for the stratospheric height section between 20 km and 25 km, where 13'667 profiles are included.
The GRUAN data processing is based on the raw data from the physical radiosonde's sensors, namely temperature, relative 95 humidity, positioning data (GNSS), and also pressure if an on-board sensor is present. Corrections to the raw data for known or experimentally evidenced systematic effects are applied. For example adjustments from pre-flight ground checks, corrections of sensor time lags, or solar radiative effects. Some intermediate variables are in turn calculated (e.g., effective air speed or ventilation) as components of the correction algorithms. A number of secondary variables are finally derived, for example altitude, geopotential height, water vapour content, or wind components. At different processing stages, smoothing filters are 100 applied for estimation and separation of noise components of the signal. Through all these steps the regular grid of the measured raw data is maintained, that is, all variables and uncertainties in the product variables are given with the original high resolution.
This procedure inevitably leads to certain technical difficulties if data gaps randomly or intermittently occur. For example, smoothing with certain filter kernel lengths easily may introduce effects which are difficult to handle when running over gaps on the regular grid or -even more -when running into larger gaps comparable to or exceeding the actual kernel length. The same 105 applies to uncertainty estimates to be associated with the averaged (smoothed) values. Another example is related to magnitudes which are calculated cumulatively with height, such as pressure derived from positioning, temperature, and humidity data, or the integrated water vapour content. As a consequence there may be processing-related irregularities or deviations in the profile data and uncertainty estimates, the systematics and extent of which is difficult to predict. Depending on the purpose for which the GRUAN data product are further used (e.g., process studies based on high-resolution data, or average-based long-term 110 studies for climate) such systematics may have different impact.

Interpolation uncertainty
In this section, formulas of the uncertainty for both linear and stochastic interpolation are considered under some stochastic assumptions about the data generation mechanism.
In particular, considering a radiosonde flight, we assume that t = 1, ..., T is the flying time in seconds from take off and y(t) 115 is the observed temperature in Kelvin given by the following measurement error equation (1) In model (1), s(t) is the unobserved "true" temperature with a local dynamics described by a Gaussian Process (GP) characterised by a power exponential autocovariance function (Cressie and Wikle, 2011;Rasmussen and Williams, 2006). Hence, conditionally on some unobserved time-dependent atmospheric conditions denoted by a(t), the GP y(t) has the following 120 autocovariance function: where p = 1, 2, the dependence on a(t) is omitted in the right hand side for notational simplicity, and I is the indicator function, that is I = 1 if t = t and zero else.
In (2), the variance of y(t) is given by where σ s > 0 is the standard error of s(t), and σ ≥ 0 is the measurement uncertainty, σ 2 = E( 2 ). For the instruments installed on the Vaisala RS41 it is known that the sensor-intrinsic "noise" of a temperature sensor is very small (< 0.01 K), hence we expect to find a small σ for the data of this paper. In addition, θ > 0 represents the atmospheric persistence range.
The GP is characterised by the parameter set Ψ = (θ, σ s , σ ), which is assumed to be slowly varying in time, hence charac-130 terising locally the atmospheric conditions a(t): Note that, from the practical point of view, the random error is a Gaussian white noise and σ represents the random uncertainty, while (vertically) correlated errors could be confused with s(t). This point will be considered further in Section 8.

Linear interpolation
Considering an observation gap in the interval (t − , t + ), the linear interpolator at time t, for t − ≤ t ≤ t + , is straightforwardly defined by the following formula where, y ± = y(t ± ), and In general, the squared interpolation uncertainty is defined in terms of the true signal s(t) and may be related to the interpolation Mean Square Error Since, using field observations, only M SE y (t) may be directly estimated, if the measurement uncertainty σ is unknown, estimating u(t) may be an issue, and a statistical modelling approach is important.
Assuming that the true signal s is a GP as above discussed, the Appendix shows that the linear interpolation uncertainty given in Equation (6) may be computed by the following Standard Error (SE) formula: where, with abuse of notation, α = α(t). Note that, SE(t) 2 = u(t) 2 if the GP assumption is satisfied, but two different symbols are used because in Section 8 this assumption will be relaxed. Equation (7) defines a function of t which depends on the atmospheric persistence modelled by γ and the gap size t + − t − .

155
Since γ is not continuous in zero, the same happens to SE(t) at the gap interval borders. Figure 2 considers the case where s(t) is a white noise, that is γ(h) = 0 for h = 0 and γ(0) = σ 2 y . At the gap borders, the interpolation is error free, m(t ± ) = y ± , and the uncertainty is u(t ± ) = σ . For t strictly inside the gap interval, we have where the minimum is achieved in the center of the gap interval. In this particular case, the uncertainty range does not depend 160 on the gap size.
The above thresholds may be overcome in the presence of correlation. In general for a GP with θ > 0 the uncertainty depends both on the GP characteristics and the gap size. As an illustration, using σ s = 0.5K, σ = 0.01K, and θ = 3 s, Figure 3 shows how the interpolation uncertainty depends on the gap size and on the distance from the observations in presence of short correlation range. More interestingly for applications, Figure 4 shows that the linear interpolation uncertainty strongly depends 165 on the correlation range.

Gaussian Process interpolation
The assumption that the temperature profile y(t) is a realisation of a GP may be extended to cover for a non constant mean so that, using some basis functions h(), model (1) is rewritten as with parameter set Ψ = (β, θ, σ , σ s ). In this context, Equation (3) defines the variance of y(t) conditional on h(t) β, namely V ar(y(t)|h(t) β). Let us denote the set of all non missing observations during the radiosonde flight by Y , the matrix of the corresponding basis functions by H, and assume that Ψ is known. Then the GP interpolation of a missing observation at time t * is given by the well known conditional expectation formula

175
where Σ Y ,Y is the covariance matrix of the good observations Y |Hβ, and Σ y(t * ),Y is the covariance vector between the missing observation y(t * )|h(t + ) β and Y |Hβ. In addition to point estimation, the GP approach provides also the interpolation standard error: which can be used as an estimate of the interpolation uncertainty, provided the GP is a good description of the problem under 180 study and Ψ is approximately known.

Data
Two datasets provided by the GRUAN Lead Centre (www.gruan.org/network/lead-centre), and related to the seven GRUAN stations of Table 1, are considered here. One is named F ew_nan in this paper and contains 276 temperature profiles characterised by "little" missing data. The second one, named M any_nan, contains 273 profiles with "many" missing.

185
As a preliminary analysis of the "bad" dataset M any_nan, Figure 5 gives the distribution of the fraction of missing data per launch. The average missing fraction is 0.13, and the average gap length is 3.6 s. These values will be used to set the parameters of the simulated gap patterns of Section 5.  For further interpolation analysis, those profiles in F ew_nan with very little missing gaps are selected. In particular, the L=177 launches which have gaps shorter than 5 seconds, and a total of less then 10 missing values per profile have been used 190 in this paper. The profile duration distribution is depicted in Figure 6, with an average profile duration of about 6000 seconds.
This gives a total of more than one million measurements, which will be amplified using the bootstrap technique of Section 5.

Block-bootstrap cross-validation scheme
The block-bootstrap is a well-known technique (Politis and Romano, 1994;Mudelsee, 2014) for generating synthetic time series replicates, and, in this paper, is used to construct the cross-validation scheme. Let us consider a fully observed temperature 195 profile, without missing values and, hence, measurements y taken every second from take off, t = 1, ..., T : Y = (y(1), ..., y(T )).
This section presents a rule for partitioning each original profile Y as follows 9 https://doi.org/10.5194/amt-2020-161 Preprint. Discussion started: 3 July 2020 c Author(s) 2020. CC BY 4.0 License.  10 https://doi.org/10.5194/amt-2020-161 Preprint. Discussion started: 3 July 2020 c Author(s) 2020. CC BY 4.0 License. where Y L is the learning set -used for fitting -and Y * is the validation set -used for testing and bootstrap-correction. In order to construct the testing set, n G gap sequences of average duration µ G [ s] are extracted from the temperature profile Y 200 and moved to the testing set Y * . Observe that, if the testing size (average) fraction is denoted by f , then n G = T × f /µ G .
The gap scheme is obtained by randomly generating and sorting the n G gap starting points 1 ≤ t * 1 ≤ ... ≤ t * n G ≤ T and by building, for each of them, a gap sequence t * j , ..., t * j + g j where the gap duration g j is a Geometric random variable with mean µ G . In particular the length g j is truncated at t * j+1 −t * j −1 205 to avoid overlapping among different gap sequences. Let the resulting testing set index be denoted by t * . Ignoring above truncation, t * has random dimension n * = n G + n G j=1 g j and expected dimension E(n * ) = T × f . Hence, the partitioning rule in (10) is defined by the testing set Y * = (y(t), t ∈ t * ) and the learning set Y L = (y(t), 1 ≤ t ≤ T, t / ∈ t * ).
We are interested in collecting information about the interpolation error in a dense vertical grid, even if the testing fraction f is small. To do this in the application developed below, the above random extraction process is repeated B times. So that for where j = 1, 2 denotes the linear or the GP interpolation respectively.
.., B is used first to estimate the GP model coefficients Ψ by the maximum likelihood method as explained in the next section and denoted byΨ. Then, the interpolated valuesŷ(t * ) = m 2 (t * |Ψ) are computed for the simulated missing times t * in the test data set, Y * b , and the cross-validation interpolation errors are computed 225 as follows: e = e(t * |s, l, b) =ŷ(t * |s, l, b) − y(t * |s, l).
As a result, the interpolation MSE and the corresponding Root MSE (RMSE) are classified by station, altitude and gap length: where ALT is the atmospheric output layering with 1km resolution and avg(·|s, ALT ) is the average of all the cross-validation 230 terms with alt ∈ ALT , launched from station s and generated using gap size µ G .

Modelling details
The GP interpolator depends on the local structure m 2 (t) = m 2 (t|Ψ a(t) ), where Ψ a(t) is as in Equation (4). In order to make teh local GP modelling feasible and computationally efficient a block-partitioning structure has been assumed. This amounts at dividing the atmosphere in layers, which may be different from the output layers of the previous section. Each atmospheric 235 layer identifies a block and the variance-covariance matrix for the entire profile is assumed block diagonal with a constant parameter set Ψ a in each layer block. This is a special case of the spatial partitioning approach (Heaton et al. , 2019), but continuity at the layer borders is ignored here because borders have been deliberately located far from the gaps.
The GP model selection considered the two autocovariance functions γ in Equation (2), various basis functions h(), and various layering's of the atmosphere to define the appropriate concept of local model Ψ a (t) of Equation (4). For each layer a, 240 local estimation has been performed using the maximum likelihood method. The above alternatives have been optimised using the RMSE applied to the block-bootstrap replicates of Section 5.
Considering the layering problem, the results where little sensitive to layer size variations and a 400" layer size has been used, as it provides both a reasonable computing time and a satisfactory atmospheric adaptation. The exponential autocovariance function with p = 1 resulted in a smaller cross-validation RMSE comparing to the square exponential one (p = 2).

245
The best results for the basis functions have been obtained with a piecewise linear function of time. In this regard, also other predictor set-up have been considered: a piecewise quadratic function of time and vector predictor set-ups including altitude, coordinates and wind. Using these more complex models did not result in relevant improvement of RMSE or, worst, it resulted in problems of singularity of the information matrix at various combinations of stations and layers. Hence, invoking Okham's razor and looking for a robust and general model set-up, we concluded for the simplest piecewise linear function of time.

Results
The bootstrap champaign of this section is aimed at assessing the uncertainty of the linear interpolation, Equations (5) and (7), and of the GP interpolation, Equations (8) and (9). The cross-validation design is based on a 3 × 3 combination of gap sizes µ G and missing fractions f , centred on the characteristics of the M any_nan dataset. In detail, we use µ G = 4, 10, 30, 60 sec and f = 0.05, 0.13, 0.20. Moreover, in order to have uncertainty estimates with a high vertical density, the 2-fold block-255 bootstrap validation of Section 5 is replicated B = 50 times giving a data set with more than 51 million measurements for each combination of µ G and f . Figure 7 depicts the overall linear interpolation uncertainty at each GRUAN station for the above 3 × 3 simulation design using the RMSE. The clustered pattern of the nine curves clearly shows that the missing fraction f has a minor influence on the uncertainty in the range 0.05 − 0.20, which is the range of interest for meaningful practical applications. Hence for the rest 260 of the paper, we consider only the M any_nan dataset missing fraction, f = 0.13. Table 2 summarises the RMSE of both the linear and GP interpolation. Overall, the average interpolation uncertainty is smaller than 0.1K for little gaps (µ G = 4"), increases to about 0.16 K for medium gaps (µ G = 10"), and increases further to 0.35 K and 0.58 K for large and very large gaps (µ G =30" and 60") respectively. Considering jointly Figure 7 and the latter that the two interpolation approaches have a very close RMSE. In fact, not only they have close performances, but, for any practical purpose, they are also exchangeable, since the mean absolute difference between the two is smaller 0.01 K. Hence in the rest of the paper, we do not replicate figures and results for both interpolation methods.   14 https://doi.org/10.5194/amt-2020-161 Preprint. Discussion started: 3 July 2020 c Author(s) 2020. CC BY 4.0 License. In order to re-interpret the GP-based linear interpolation uncertainty formula of Figures 3 and 4, we consider the ensemble of all the estimated local GP model parameters set Ψ from the entire cross-validation exercise. Coherently with the known small intrinsic error declared by Vaisala, the top panel of Figure 9 shows very small values of σ . Moreover, from the second panel of 275 the same figure, we see that the values σ s < 1 are common and in particular σ s = 0.5 used in Figures 3 and 4 is quite plausible.
Eventually, the bottom panel of Figure 9 shows that the correlation range θ may be easily between one and 15 minutes.

Interpolation distance
In general, the connection between the uncertainty curves of Figures 3 and 4 and the cross-validation evidence is worth to be studied. Considering both the gap size and the distance from the observations at various altitudes gives rise to hard-to-manage 280 curve plots and a multiplicity of results. For this reason, the subsequent analysis is based on the "interpolation distance" in sec, which is denoted by d and is given by the geometric mean of the temporal distances of each missing data from the closest observations y − and y + in the notation of Section 3. where avg(·|ALT, d) is the average of all the cross-validation terms with alt ∈ ALT and interpolation distance d. We note that, in order to have high sampling information for both low and high interpolation distances, the graph is obtained by merging the block-bootstrap simulations obtained for µ G = 10 and 30 s. Moreover, the graph is limited to 70" because there is a reduction of cross-validated data, especially at high altitudes. Of course, using the same approach, longer interpolation distances may be easily explored.

290
In addition, Figure 11 depicts the corresponding graph for the linear interpolation quadratic average of SE(t * ) = SE(t * |s, l, b), given by Equation (3) at cross-validation time t * , station s, launch l and bootstrap replication b, namely SE(d|ALT ) = avg(SE(t * |s, l, b) 2 |d, ALT ).
The corresponding graph for the GP-SE of Equation (9) is not reported here because, not only the two interpolation methods are exchangeable, as noticed above, but also their SE's give very close results, with a mean absolute difference between the 295 two smaller than 0.005 K. It may be noted that, although the above two graphs have a similar increasing behaviour, the SE systematically underestimates the interpolation uncertainty. This is due, primarily to the GP model approximation for the present case study and, secondarily, to estimation uncertainty. Hence, we propose a corrected uncertainty estimate given by This semi-parametric bootstrap uncertainty estimate extracts information both from the average cross-validation performance 300 at a certain altitude and interpolation distance and from the single profile behaviour approximated by the GP process.  . The x-axis is given by the geometric mean of the distances of each missing data from the closest "good" data, y − and y + . The graph is obtained my merging the data sets with average gap sizes µG = 10 and 30 s.

Practical aspects
As an illustration of the method, the profile of Sodankylä site on 2017-03-03 12:00 is considered in Figure 12, left panel. This profile has T=4722 measurements and no original missing values. Using the block-bootstrap, 563 measurements have been deleted and considered as pseudo-missing generating gaps between 1 and 24 s to be interpolated. From a practical point of view, 305 such a missing rate and gap lengths can be considered a relatively common case, yet serious, for interpolation. Figure 12, right panel, shows both ± the GP uncertainty (7), and ± the Bootstrap uncertainty (14), computed at the interpolated pseudo-missing values. In doing this computation, formulas (12) and (13) are implemented as lookup tables (LUT) with entries geometric distance and altitude. Figure 13 focuses on the above profile around 22 Km height and shows the interpolation uncertainty of a single point gap, two small gaps and three larger gaps.

310
It follows that the implementation of a GRUAN data processing giving interpolated temperature profiles with their uncertainty requires some efforts which are divided into two different phases. First, a massive GP off-line computation is needed fact, the approach is quite general, and model selection and optimisation are data-driven. Hence similar results may be expected for temperature profiles obtained by other instruments, provided that vertical resolution and instrumental error are comparable to the present case. Further, similar results are also expected for other smooth variables, such as pressure.
On the other side, the interpolation uncertainty could be larger for those ECV which are known to have large variations also in the small scale. For example, relative humidity commonly shows highly intermittent profiles in the troposphere with 335 very large and very fast changing gradients. In these cases, we can expect that the cross-validation uncertainty could be large even for small gaps. In addition, the vertical autocorrelation could have a shorter range and the corresponding GP model could provide interpolation uncertainties close to the white noise case considered in Section 3.
Code and data availability. TEXT The underlying MATLAB code is available from the author upon request. The data are available from the GRUAN Lead 340 Center, www.gruan.org.
Competing interests. The authors declare that they have no conflict of interests