A variational regularization of Abel transform for GPS radio occultation

In the Global Positioning System (GPS) radio occultation (RO) technique, the inverse Abel transform of measured bending angle (Abel inversion, hereafter AI) is the standard means of deriving the refractivity. While concise and straightforward to apply, the AI accumulates and propagates the measurement error downward. The measurement error propagation is detrimental to the refractivity in lower altitudes. In particular, it builds up negative refractivity bias in the tropical lower troposphere. An alternative to AI is the numerical inversion of the forward Abel transform, which does not incur the integration of error-possessing measurement and thus precludes the error propagation. The variational regularization (VR) proposed in this study approximates the inversion of the forward Abel transform by an optimization problem in which the regularized solution describes the measurement as closely as possible within the measurement’s considered accuracy. The optimization problem is then solved iteratively by means of the adjoint technique. VR is formulated with error covariance matrices, which permit a rigorous incorporation of prior information on measurement error characteristics and the solution’s desired behavior into the regularization. VR holds the control variable in the measurement space to take advantage of the posterior height determination and to negate the measurement error due to the mismodeling of the refractional radius. The advantages of having the solution and the measurement in the same space are elaborated using a purposely corrupted synthetic sounding with a known true solution. The competency of VR relative to AI is validated with a large number of actual RO soundings. The comparison to nearby radiosonde observations shows that VR attains considerably smaller random and systematic errors compared to AI. A noteworthy finding is that in the heights and areas that the measurement bias is supposedly small, VR follows AI very closely in the mean refractivity deserting the first guess. In the lowest few kilometers that AI produces large negative refractivity bias, VR reduces the refractivity bias substantially with the aid of the background, which in this study is the operational forecasts of the European Centre for Medium-Range Weather Forecasts (ECMWF). It is concluded based on the results presented in this study that VR offers a definite advantage over AI in the quality of refractivity.


Introduction
The Abel transform pairs (Abel, 1826) are widely used to reconstruct radially (or spherically) symmetric physical parameters 20 from their line-of-sight projections in a variety of disciplines in engineering and science. In the Global Positioning System (GPS) Radio Occultation (RO) technique, the inverse Abel transform (often called Abel inversion) of measured bending angle in particular has become a cornerstone serving as the standard means of deriving the refractive index. Knowledge of refractive-index structure in the atmosphere is valuable for numerous applications relevant to weather and climate.
The beauty of Abel inversion (AI hereafter) lies in mathematical exactness and conciseness. The shortcoming on the other 25 hand is that it integrates flaws as well as signal components in the measurement. Therefore, the analytical inversion works well for measurements with exceptionally rigorous accuracy and GPS RO is believed to provide such a measurement.
Nonetheless, RO measurements are subject to non-negligible errors arising from diverse sources to which AI is sensitive.
More fundamentally, the premise of AI (i.e., spherical symmetry of the atmosphere) is never strictly fulfilled. Thus, the Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. spherical asymmetry in the atmosphere is one of sources of the measurement error. A major concern in AI is that the measurement error is accumulated and propagated in vertical direction without being controlled. Consequently, a single corrupted piece of the measurement affects not just the location of the particular datum but a wide area that the Abel kernel dictates, and thereby deteriorates the solution even in the region that received RO signals are clean. Previous validation studies suggest that the RO refractive index derived through AI is of high quality, which however does not necessarily mean 5 that AI yields the best possible refractive index in the presence of measurement error.
An alternative to the analytical inversion is finding the inverse solution of the forward Abel transform (FAT) numerically. In our setting, FAT simulates the bending angle given the refractive index. Thus, the numerical inversion (note that it differs from numerical implementations of the analytical AI) attempts to find the refractive index that maps into the measured bending angle through the FAT. Unfortunately, inverse problems of the kind of FAT are generally known to be ill-posed in 10 the presence of measurement error, e.g., the solution may not exist or may not be unique. Therefore, it is challenging or even infeasible at times to obtain useful solutions by applying straightforward algebraic methods to the underlying discrete operators.
Regularization methods such as the Tikhonov regularization (Tikhonov, 1963) are widely used to tackle the ill-posedness in inverse problems (Schuster et al., 2012, and references therein). Based on the notion that it is practically undesirable to fit 15 every detail of noisy measurements in solving the inverse problem numerically, the approach considers an optimization problem in which the solution describes the measurements closely, while the solution's behaviour is regularized as per a priori knowledge. The solution is sought by minimizing a joint function that consists of the data fidelity term (which gauges the discrepancy between measurement and its model counterpart) and the penalty term (which regularizes the solution).
Besides warranting that the inverse problem remains well posed, regularization methods are known useful for reducing 20 adverse effects of measurement error.
Regularization methods have proven useful for a variety of inverse problems in different forms. However, it is seldom straightforward to apply those to real-world problems. First of all, the computational cost, which rises rapidly as the number of model parameters increases, can be expensive when applied to GPS RO. That is because the limb-viewing geometry offers high-resolution measurements; and, one typically attempts to retrieve a profile of refractive index in as much detail as 25 the data permit. More importantly, the regularization must be customized for individual applications in order to be effective.
For instance, different observing systems differ in measurement characteristics and solution's desired behaviour. Thus, the specifics of an optimal regularization differ from one problem to another and those for the Abel transform in GPS RO are underexplored.
In this study, we propose and study a variational regularization (VR hereafter) for the Abel transform in GPS RO. The 30 degree of benefit attainable from using regularization methods depends on the quality of measurement. That is, while comparable to AI in the solution for error-free measurements, regularization methods are more promising for noisier measurements. In order to attain the benefit maximally, however, it is essential to estimate the uncertainties of RO measurement and a priori realistically and take them into account properly. In this study, the data uncertainties are estimated Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. by a reliable method, and they are incorporated into the proposed regularization method through error covariance matrices.
This differs from classical regularization methods, which in general do not make use of the data uncertainty. The quality and vertical coverage of RO measurements differ significantly across occultation events varying with the occultation geometry, atmospheric conditions encountered, and other environmental factors. The variability causes a well-known difficulty in classical regularization methods as to determining the regularization parameter, which is a constant that controls the tradeoff 5 between the data fidelity term and the penalty term. The proposed method on the other hand is highly adaptive to those variations because the error covariances offer a proper scaling for the data fidelity and penalty terms. Moreover, the proposed method reduces the computational cost significantly, solving the underlying inverse problem iteratively by means of the adjoint method, which is a very efficient way of calculating the gradients of the cost function with respect to all control parameters at once. 10 The quality of AI-derived refractive index also depends on how the AI is numerically implemented. However, varying implementations might not cause significant difference in the refractive index because the data resolution of RO bending angle is sufficiently high so that the numerical error due to the discretization of AI is limited in the magnitude. Moreover, RO data processing centres worldwide are now using polished algorithms aiming at attaining the best data quality, which includes well-tuned implementation of the discrete AI. The additional quality improvement in the refractive index obtainable 15 from further polishing the discrete AI is thus likely to be marginal at most. This necessitates a breakthrough that can improve the quality of refractive index to a level beyond AI can tender. This study explores the possibility with the proposed method of variational regularization. Section 2 describes the methods of obtaining refractive index from measured bending angle.
Section 3 demonstrates the performance of the proposed method against AI with a synthetic sounding. Some issues in practical application of the Abel transform are also discussed. In Sect. 4, a real-data validation with respect to two sets of 20 radiosonde observations is given. We conclude with a short summary of our results in Sect. 5.

Abel inversion
The total transpired phase path of a radio wave that propagates through the atmosphere between a transmitter and its receiver can be described by an integral equation (e.g., Wallio and Grossi, 1972) where is the phase path; is the radius from the Earth's curvature center; and, is the refractive index. A change of variable can show that Eq. (1) is equivalent to standard form of the Abel transform (Bracewell, 1978). A geometrical interpretation of Eq. (1) is that the Abel transform of , right-hand side (RHS) term, is the projection of onto the space of the traverse coordinate, , which indicates the nearest approach of the line of sight to the Earth's curvature centre. 30 Being an Abel transform, the analytical inverse of Eq. (1) also exists (Ahmad and Tyler, 1998): Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License.
The inverse Abel transform, Eq. (2), provides a straightforward solution to the reconstruction of the refractive index given measurements of phase path. However, the inverse transform is of limited usefulness for practical applications with realworld data because the differentiation operator in RHS exacerbates small perturbations in (e.g., phase noises or artifacts introduced by arbitrary noise mitigations), causing the derived unstable. 5 Another limitation of Eqs. (1)-(2) is that they are valid for a thin, spherically symmetric atmosphere in which the ray's path can be adequately approximated by the straight line that connects the transmitter and receiver. A variant of Eq. (1) suitable for dense, optically stratified media (Fjeldbo et al., 1971) is: where is the ray's bending angle; and, the impact parameter, , is an invariant measure (known as the Snell's constant) in 10 the Bouger's law (the Snell's law generalized for spherical symmetry). The impact parameter is defined as: = sin = @ @ , where is local zenith angle of the ray and @ is the refractive index at the tangent point @ . With a change of variable, = , Eq. (3) can be rewritten as: The corresponding inverse transform is: (5) Equation (5) has the advantage of not having the differentiation operator in RHS, in addition to accounting for refracting ray paths. Consequently, Eq. (5) is the AI commonly used in GPS RO. A caution is that GPS RO does not measure the bending angle directly. Accordingly, the bending angle must be estimated in some ways in connection with the Doppler shift, which again relates to the phase derivative in either time or frequency domain. Therefore, the AI using Eq. (5) is not entirely free 20 from the differentiation operator. Instead, the procedure of retrieving the refractive index from the phase path, equivalent to Eq. (2), is split into two sequential steps: bending angle estimation and subsequent AI. The differentiation operator is put to use before or within the bending angle estimation, although it does not appear explicitly in Eq. (5).
While real-world measurements cannot be flawless, AI considers them to be complete. Consequently, Eq. (5) integrates not only signal components but also the noises contained in the measured bending angle. Indeed, they are integrated in the same 25 manner without being distinguished from each other. Hence, the quality of derived refractive index is strongly dependent on the signal-to-noise ratio (SNR) of the bending angle. For an effective use of Eq. (5), it is thus crucial to control the measurement noise properly. In practical applications, one may rely on external means of noise mitigation while preparing the bending angle for Eq. (5). The limitation is that the information about the measurement noise apprehended by the pre-Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. mitigation is simply thrown away. Hypothetically, the information could have been utilized in some ways while obtaining the refractive index.
In the presence of measurement error, the AI using Eq. (5) generally does not yield the best possible refractive index; and, corrupted measurement often causes Eq. (5) ineffective or unreliable. An alternative to AI is seeking for the inverse solution of FAT numerically. However, doing so is challenging because the inverse problem of the form of Eq. (4) is generally 5 considered to be ill-posed violating the so-called Hadamard's stability criterion (Hadamard, 1902). The standard approach to addressing the ill-posedness in inverse problems is regularization methods, described in the following.

Tikhonov regularization
Regularization methods have been used for a variety of inverse problems in different forms. Detailed description of the methods is beyond the scope of this study. Here, a sketch of the Tikhonov regularization (TR hereafter) is provided in the 10 context of GPS RO data processing. Amongst the regularization methods, TR is perhaps the most widely used method and is indeed the very method that opened up the concept of regularization. The general purpose of TR is to solve ill-posed inverse problems in which the forward operator = defines a mapping : → , between the solution (model) space and the data (measurement) space . Here, is the state vector consisting of model parameters and is the vector of modeled observation. In our setting, the forward operator is Eq. (4), and, and hold the refractive index and bending angle, 15 respectively. The method solves a minimization problem of a real-valued function such that the minimizer, , is a suitable approximate solution: where is the data fidelity term; is the penalty (also called regularization) function; and, is the regularization parameter.
The data fidelity term measures the misfit between measured and modelled observations, whereas the penalty term 20 weighs the degree of irregularity in the solution. The regularization parameter controls the relative contribution of the two terms to . A widely used form of is the squared S norm, = − @ S , where @ contains measured bending angle. A popular choice for aimed at noise reduction is a seminorm, = Ñ S , where Ñ is the gradient operator.
Determining the tradeoff between the two terms, a decreasing steers Eq. (6) toward (4). Therefore, the solution of TR approaches AI as vanishes. 25 Noises in measured bending angle give rise to high-frequency artifacts in the retrieved refractive index. Although the shown above suppresses such high-frequency modes effectively, it is difficult to choose a proper that can optimally control the degree of solution's smoothness. The difficulty is not different from what one may encounter designing an optimal filter. Moreover, the particular does not account for any height dependency of the refractive-index gradient, which is very strong in the real atmosphere. For instance, Ñ near the Earth's surface can be a few orders greater than that in the stratosphere. 30 Another reason that hampers the settling on is that usually has considerable latitudinal and seasonal variations. Likewise, Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. the shown above tacitly assumes that the uncertainty of measured bending angle is unvarying in space and time. In reality, the quality of RO measurement is influenced by environmental factors. For instance, the tracking of GPS signal performs less optimally under certain unfavourable atmospheric conditions. The influence of the environmental factors varies with height, latitude, and season, etc., as do the atmospheric structure and processes. The regularization parameter is supposed to take into account the spatiotemporal variations because neither nor does. However, is typically a constant (i.e., 5 scalar) and it does not have any structural information inside regarding the spatiotemporal variations.
Given that only noisy measurements are available in the real world, TR is generally advantageous over AI. For instance, TR allows one to merge noise control into the inversion as shown by the penalty function . Nonetheless, the TR in the form described above is most useful for general observations and media whose characteristics are unknown. In order to make the regularization more effective, it is desirable to maximally incorporate all significant prior information that is specific to the 10 problem of interest. In our study, the prior information is about the character of RO measurement, and about the desired structure and behaviour of the refractive index to be retrieved. In addition, regularization methods can be costly when applied to actual RO data. It is thus crucial for practical applications to reduce the cost to a feasible level. We aim at addressing these issues through the variational regularization proposed in this study, as described below.

Proposed method of variational regularization 15
The variational regularization (VR hereafter) proposed in this study searches for the solution, , by minimizing a cost function defined as follows: where T is a priori of ; and are error covariance matrices of T and @ , respectively; and, superscripts "T" and "-1" indicate transpose and inverse of a matrix, respectively. Conceptually, the method seeks the optimal solution that replicates 20 the observation as closely as possible (compelled by the second RHS term) in the vicinity of a priori state (constrained by the first RHS term). The background (first RHS) term corresponds to the penalty term of TR in the role and the observation (second RHS) term is equivalent to the data fidelity term. An obvious distinction from TR in the formulation is that VR makes use of the error covariance matrices that conveys the prior information about the uncertainty of @ and T . In our study, and are built for each month on a 5° ´ 5° latitude-longitude grid. This allows VR to account for spatiotemporal 25 variation of the data uncertainty. Specifically, the normalization of − T and @ − with their respective error, which is the square root of the diagonal elements of and , eliminates the dependency of RHS terms on height, latitude, season, and so forth. In addition, the off-diagonal elements of the matrices permit VR to take into account spatial error correlations, which is essential for the background term.
The forward observation operator used here is Eq. (4), which maps the state variable (refractive index) onto the observation 30 (bending angle). It must be emphasized that the state variable (control variable to be precise) of VR resides in the impact S Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. parameter space where the observation exists, as can be conjectured from Eq. (4). This differs from meteorological data assimilations the state variables of which are usually the same as those of the prediction model: temperature, moisture, pressure, wind, and so on. In the data assimilations, the location of the state-vector elements is known only in relation to the model's native grids. For the reason, in order to place measured bending angles in the model space and vice versa, it is necessary to compute the model-side impact parameter, = . The problem is that the observation-side impact parameter, 5 = sin , is hypothetical in that it changes with the ray's direction even at the same point in the Euclidean space. Also, is not a definite measure of the position, relying on the never-perfect model's refractive index. Hence, the modelled impact parameter always contains some errors. Because data assimilations place measured bending angles in the model space by using the erroneous x as the coordinate, the observation's position in the model space cannot be determined precisely not because of its own uncertainty but because of the model's error in the refractive index. The positional error is an extra error 10 source, leading to a suboptimal use of the observation. This is known as errors-in-variables (EIV) problem in which the error in the independent (or coordinate) variable causes an apparent error in the dependent variable (observation).
In the proposed method, the EIV error is attributed to the uncertainty of the background refractive index because the model space coincides with the observation space. In VR, the solution's geometric height keeps changing implicitly during the iterative minimization of the cost function, whereas the impact parameter assigned to the solution remains fixed. Upon the 15 completion of the minimization, the solution's geometric height is determined as: the local curvature radius of the Earth. Therefore, the solution's height in the proposed method is undecided until the solution is acquired, which is the same as in AI. The posterior height determination (PHD hereafter) reduces a substantial portion of retrieval error when viewed in the height coordinate. More details on the topic are presented in the next section.
In VR, the background term is based on the distance between the state vector and the background, − T , meaning that the 20 solution is expected to be around T . In that sense, the quality of the background matters, although the final solution does not depend on every detail of the background. A good background furnishes faster convergence to the solution and assists in attaining the desired global minimum of the cost function. A factor of particular importance in VR is the adequacy of error covariance matrices. That is, the error covariance matrices must be factual, representing well the actual error characteristics of the background and observation that are used for VR. For instance, while a climatology is usable as the background, it is 25 not easy to estimate or define the relevant data uncertainty precisely. In that regard, short-term forecast of numerical weather prediction (NWP) models is one of the best sources of the background. The forecasts are routinely available and of a reasonably good quality. More importantly, the forecast error is largely uncorrelated with observation error. This allows one to construct realistic and by comparing the forecast and observation for a large number of samples. In this study, the approach proposed by Hollingsworth and Lönnberg (1986) is applied to closely located (< 3 h and < 300 km) pairs of RO 30 soundings  A practical difficulty in solving inverse problems of a large size is the computational cost. It is unfeasible to perturb individual elements of the control vector arbitrarily in all directions and sizes, and then search for the very combination of the perturbations that leads to the minimum of the cost function. In order to reduce the cost, the adjoint technique (Lewis and Derber, 1985;Le Dimet and Talagrand, 1986) is employed in this study. The method efficiently computes the steepest 5 gradient of the cost function with respect to all elements of the control vector at once, which is needed for the optimization algorithm. In this study, a quasi-Newton limited-memory algorithm for large-scale optimization (Zhu et al., 1997) is used.
The proposed method makes use of the control-variable transform (Parrish and Derber, 1992), which reduces the conditioning number of . The square root of needed for the transform is computed with the method proposed by Kaiser (1972). 10

Test with a synthetic sounding
It is hypothesized in this study that analytical AI is not necessarily the best way of retrieving refractive index from bending angles with error. We deliberate that meticulous utilization of the prior information via the variational regularization might be beneficial for the retrieval. Meanwhile, a common issue that arises while verifying hypotheses with real-world data is that the verifying data accompany their own error, which often impedes drawing decisive conclusions on the verification. In 15 order to overcome the difficulty, we begin the verification with a synthetic profile of bending angle. The synthetic case provides the known true solution against which the results of inversion methods can be compared with reduced ambiguity.
We intend to consider large-amplitude errors so as to emphasize their influence on the solution's quality and to assess the relative robustness of the inversion methods to the erratic observations. This section is also purported as an extended description of the methods but with a tangible example. 20

Data generation
The tracking of RO signal is often challenging, in particular in the lower tropical troposphere where sharp refractivity gradients frequently exist. As an example, a high-resolution radiosonde sounding, observed at a tropical site in Nauru (0.52° S, 166.93° E) at 12:00 UT on 3 March 2011, shows a complicated structure in the refractivity mainly due to abrupt smallscale moisture variations across multiple inversion layers (Fig. 1). The station was one of the Global Climate Observing 25 System (GCOS) Reference Upper-Air Network (GRUAN) sites until closed in August 2013. The refractive index in the neutral atmosphere can be approximated as (Smith and Weintraub, 1953): where is the refractivity; is temperature in K; is (total) pressure in hPa; f is water vapor pressure in hPa; and, / =77.6 hPa K -1 and S =3.73·10 5 K 2 hPa -1 are coefficients.
The discretization of Abel transforms causes significant numerical error when the data resolution germane to the discrete 30 transform is poor. In order to reduce the numerical error, we used the highest data resolution available to us for the Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. radiosonde sounding. Despite an important factor for practical applications, the numerical error is not the focus of this study.
Thus, no extra effort (e.g., use of higher-order numerical schemes) is made in this regard. The sampling rate is one second in time, which corresponds to about 4 meters in height interval. The bending angle simulated with the radiosonde data and using Eq. (4) also presents rapid variations (green line in Fig. 2a). This bending angle is assumed to be absolutely accurate (error-free) and is referred to as the perfect (true) observation. The measurement error of the radiosonde data carried forward 5 into the bending angle is considered as legitimate small-scale variations in the true atmosphere.
We then generate a synthetic observation by adding suppositional measurement errors to the perfect observation. The measurement error in bending angle is assumed to follow a first-order autoregressive process and modeled as: where is again the impact parameter and is the length scale of error correlation. Assuming that the vertical error correlation for the bending angle is very weak, we set to 10 m. The black line in Fig. 2a is the error-added, "measured" 15 bending angle.

Abel inversion of low-pass filtered bending angle
Noises in measured bending angle negatively affects the quality of refractive index unless properly mitigated. It is thus customary to apply data smoothing to the bending angle prior to AI. In order to mimic the practical application of AI, we use a low-pass filtered bending angle (red line in Fig. 2b) for which the fourth-order Butterworth filter (Butterworth, 1930) with 20 a cutoff wavelength of 200 m is used. The step following is obtaining refractivity soundings through AI from the true, measured, and filtered bending angles. It must be mentioned that the resulting true refractivity (i.e., the one derived from the true bending angle via AI) differs slightly from the refractivity used to generate the true bending angle (shown in Fig. 1c).
The discrepancy stems from numerical approximations made for the discrete version of Abel inversion. Figure 2c compares the refractivity error, which is the difference from the true refractivity. The red line indicates the 25 refractivity error when the filtered bending angle is used for AI. The result from "raw" bending angle (black line) is shown as the reference against which the effect of the low-pass filtering can be evaluated. A common problem with any low-pass filtering is that the degree of smoothing is hard to control. An excessive smoothing leads to a loss of observational information, whereas a minor filtering causes insufficient noise attenuation. On top of that, the noises are often nonstationary, meaning that the noise spectra vary with height in accordance with vertically varying atmospheric structure. For 30 instance, the low-pass filtering used in this study tends to reduce the refractivity error below 2 km and above 6 km, where Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. the true bending angle varies rather slowly. The error reduction is relative to the case of unfiltered bending angle (i.e., red line versus black in Fig. 2c). On the contrary, the filtering increases the refractivity error around the local peaks of the bending angle in 2-6 km height range. As different occultation events encounter different atmospheric conditions, it is impractical to design a customized low-pass filter for each occultation, which is adaptive to the local noise spectrum that varies with the height of a given sounding. Another limitation of the sequential approach (filtering followed by AI) is one-5 way flow of the information in the process. That is, AI does not pass any information about the effect of attenuated noise back to the filtering. In the regularization methods on the other hand the penalty term, which acts like a filtering, invokes a reverse communication since the term is minimized iteratively and jointly with the data fidelity term. In VR, the feedback to the control vector is given through the adjoint of observation operator.

The errors-in-variables problem 10
The EIV problem occurs when the independent variable (position in space or time) of measurements is not known perfectly.
Suppose that a particular type of radiosonde system has an offset error in the height. In that case, a flawless temperature sensor of the system is bound to produce a temperature bias, which is the height offset multiplied by the local temperature lapse rate, when monitored at a fixed height. The EIV error also emerges while comparing a measurement with its model, if the independent variable of the measurement is none of the model's coordinates. This is exactly the case for RO bending 15 angle, the location of which is defined in terms of the impact parameter. For the comparison with NWP data, for instance, it is necessary to compute the model-side impact parameter. As explained earlier in Sect. 2.3, the model's refractive index is never perfect and so the resulting impact parameter unavoidably contains errors. The impact parameter error then increases the discrepancy between measured and modeled bending angles, which is interpreted either as model error or measurement error depending on the context. The issue is closely relevant to the data assimilation of RO bending angle, at least for the 20 methods that the control variables are defined in the model space. The physical-space statistical analysis (Cohn et al., 1998) can be an exception but is not considered in this study because model-space methods are currently in prevailing use for operational weather forecasting.
The EIV error in bending angle B can be estimated as: B = 'B ': • : ; : = % , where % is model's error in the refractive index. In this study, the error estimate of 12-h forecasts of the ECMWF in the refractivity ( v T ) at the time and location of the 25 synthetic sounding (obtained through the procedure described in Sect. 2.3) is considered as the representative of % (Fig. 3a).
The corresponding : is shown in Fig. 3b; the black line in Fig. 3c indicates B , where the bending angle gradient is based on the true observation. For the particular synthetic sounding used in this study, : is overall of comparable magnitude with the typical observation error in the bending angle ( B @ ) (dashed red line in Fig. 3c; also described in Sect. 2.3). A point to notice though is that B is exceedingly lager than B @ at the heights of sharp bending angle gradient. Around those heights a 30 closer fit to the observation deteriorates the solution because not the observation but the modeled impact parameter is faulty. to be underutilized in the data assimilation. In contrast, the proposed variational regularization reproduces the true refractivity if perfect observations are given and the method is able to fit them exactly. That is because the proposed method holds the solution and observation in the same space. Simply speaking, the complete reconstruction of bending angle in the impact parameter space is the sufficient condition for the perfect recovery of refractivity.  small-scale details of the true refractivity, especially the local peaks (Fig. 1c). The proposed method (red line) yields a refractivity error considerably smaller than v T . The black line in Fig. 4b is the refractivity error resulting from the AI with 15 the "measured" (rather than the smoothed) bending angle, which is the same as shown in Fig. 1c and overlaid for comparison. The proposed method is smaller than the AI in the refractivity error almost everywhere. As described earlier, the low-pass filtering shown in Fig. 2b is unable to reduce the refractivity error (Fig. 2c) and differing choices of the cutoff wavelength do not make any notable difference in the result. The pre-filtering is ineffective in the error reduction at least for the particular synthetic sounding used in this study, where the true bending angle contains high frequencies that cannot be 20 isolated from the noises. It is difficult to further reduce the remaining errors in the solution of VR since the measured bending angle is severely corrupted by the noises.

Posterior height determination
In our experience, RO refractivity tends to agree with correlative data better than the bending angle (used to derive the refractivity) suggests. For instance, the comparison of RO data with independent verifying data (e.g., short-term NWP 25 forecast or high-resolution radiosonde observation) can be made separately in bending angle and refractivity. The comparison provides ∆ = @ − x , where @ and x indicate the observed (i.e., RO) and modelled (i.e., the verifying data) bending angle, respectively. Likewise, the comparison in the refractivity gives ∆ = @ − x . Meanwhile, it is possible to closely approximate ∆ and ∆ with one another. For instance, ∆ can be propagated into the refractivity by: ∆ B = U ∆ , where U denotes the adjoint of forward Abel transform. What is puzzling us is that ∆ is generally smaller 30 than ∆ B in magnitude. The same is observed in the comparison of error estimates (e.g., those described in Sect. An example inside the small box is shown in Fig. 5b for a detailed depiction, where the solution (denoted P) is smaller than the true refractivity (denoted T) by an error < 0. It can be shown that the height displacement of the solution due to PHD 25 is: ≃ − = −10 )_ . The solution is thus placed at a higher location (~) than the true radius ( ). What important here is that the solution's error in the height coordinate is to be perceived as the difference from (~) rather than from ( ), because the solution's true height will never be known in practice. For the particular example shown in Fig. 5b, the solution's apparent error ~~ is smaller than that the true error = , where overbar denotes the line connecting two points and vertical bars indicate the line length. By linearizing = , the slope of ~ can be shown as: The slope is indeed the critical refractivity gradient at which the ray's curvature radius is equal to the Earth's radius and is thus the threshold for the occurrence of atmospheric ducting. By inspection, the solution's true error is = Z and the , where Z is the critical gradient and is the local refractivity gradient. As ‡v 9 ‡v = 1 −ˆ‰, ~ is smaller than in the absolute size if 2 Z < < 0. In general, PHD reduces the apparent error since is known to be about -40 km -1 in the standard atmosphere (United States, 1946). An example for which the apparent error is larger than the true error is shown at the point denoted as B in Fig. 5a, where the true refractivity increases with height. It must be underscored here that not only the derived refractivity but also the height assigned to it possesses error; and, the 5 errors are negatively correlated. Therefore, the derived refractivity is also subject to the EIV problem. However, in practical circumstances (e.g., when a comparison to other data is made), the errors are attributed entirely to the refractivity and the height is assumed to be free of error. Figure 6 illustrates the response of ‡v 9 ‡v to varying . The black solid line in Fig. 6a is the true refractivity sounding in the impact parameter coordinate and black dashed line indicates a hypothetical solution in the same coordinate, which is set to 10 be 5 % larger than the true solution everywhere. In the example, the gradient of the true refractivity, , ranges from -135 km -1 to 22 km -1 (Fig. 6b). The red lines in Fig. 6a are the same refractivity soundings but seen in the height coordinate. Again, the solid and dashed lines indicate the true refractivity and the solution, respectively. As shown in Fig. 6c, ~ is smaller than 5 % ( ) except for around 3.25 km where is positive. In particular, ~ is significantly smaller than near 2.3 km where ~≈ 0 as ≈ Z . The tracking of RO signals that are affected by super-refraction, which occurs when is 15 negatively large, is known to be more challenging. This usually leads to noisy bending angle. In the heights that ≤ Z , therefore, the quality of RO refractivity is not expected to be the best. Surprisingly, Fig. 6c shows an opposing result: AI and VR have the smallest refractivity error when ≈ Z . That is because PHD purges all the apparent error no matter how big the true error is. On the other hand, PHD increases the apparent error in case of strong sub-refraction ( > 0), which is often observed at the immediate underside of local refractivity peaks. Therefore, the apparent refractivity error depends on the 20 optical structure in the atmosphere as well as quality of the bending angle data. Syndergaard (1999) described the reduction of refractivity error resulting from PHD but without relating it to Z .

Test with real data
In this section, we apply AI and VR to actual RO events and thenceforward compare the results with nearby radiosonde soundings. In doing so, we use two sets of radiosonde data in order to complement each other's weakness. In the following, 25 the data sets used here (including RO, NWP, and radiosonde) are briefly described and next the validation with respect to the radiosonde data is presented.

Data
The GPS RO data used here are made available from the Constellation Observing System for Meteorology, Ionosphere, and Corporation for Atmospheric Research (UCAR). The six COSMIC satellites have been producing 1,000-2,500 globally distributed soundings each day since the launch in April 2006 (Anthes et al., 2008). Kuo et al. (2004) and Schreiner et al. (2011) describe the CDAAC's data processing. VR uses the COSMIC neutral atmospheric bending angle as the input, as does CDAAC's AI. This is to ensure the consistency between the methods in the observation. It must be mentioned though that the particular bending angle, pre-filtered by CDAAC for its AI, is not necessarily the best choice for VR. For instance, a 5 lighter low-pass filtering can be sufficient or better. In addition, we take CDAAC's refractivity as the solution of AI, instead of carrying out an AI ourselves. This is to avoid any potential uncertainty that might arise from the difference in the practical implementation of AI. Here, we simply assume that the CDAAC's refractivity is the optimal solution obtainable via AI and from the CDAAC's bending angle. The COSMIC data used here (version 2013.3520) are available online at http://cdaacwww.cosmic.ucar.edu/cdaac/products.html. 10 The first guess provided for VR is the operational forecasts made by the European Centre for Medium-Range Weather Forecasts (ECMWF), which are on a reduced Gaussian grid (~25 km spacing in latitude and longitude) with 91 vertical levels between the earth's surface and 80 km. Although the model's spatial resolution has been increasing with time (e.g., ~9 km and 137 levels as of March 2017), we used the fixed resolution throughout our study period for the sake of convenience.
We used the empirical model of the US Naval Research Laboratory (NRL), MSIS (Hedin, 1991), above the ECMWF 15 model's top and up to 2,000 km. This is to reduce the systematic errors in high altitudes, which are notable when the upper bound of Abel transform is not high enough. The other set of radiosonde observation is the high vertical resolution data from the radiosonde stations operated by the National Oceanic Atmospheric Administration (NOAA), available online at ftp://ftp.ncdc.noaa.gov/pub/data/ua/rrs-data/.

One of radiosonde data sets used in this study is the
Beginning in 2005, NOAA began transitioning from radiotheodolite balloon tracking to GPS tracking. The data from this new system, called Radiosonde Replacement System (RRS), are recorded at 1-second (corresponds to 4-meter in height) resolution, permitting good representation of small-scale atmospheric structures. A particular advantage of this data set 30 regarding the comparison with RO data is that it provides balloon's height at every recorded moment of the flight. The operational radiosonde data (ORD hereafter) on the other hand have height reports only on the standard pressure levels.
Thus, the heights on significant levels must be estimated based on the measured values of pressure (p), temperature (T), and humidity (U). The reconstructed heights are generally suitable, but at times have larger errors due to the poor data resolution Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. as well as measurement errors in pTU. The pTU-height error is interpreted as a refractivity error when the radiosonde data are compared to RO refractivity because the comparison is to be made in the height coordinate. By using the high vertical resolution radiosonde data (HVRRD), the pTU-height error on significant levels can be avoided.

Validation with operational radiosonde data
ORD has a larger number of soundings and a superior geographical coverage than HVRRD the stations of which belong to a 5 small subset of ORD sites. Hence, ORD provides more soundings that are closely located (collocated hereafter) with RO soundings for a given period of time. It also allows the collocated soundings to be sampled at various locations and under diverse atmospheric conditions. Focusing on tropical regions, we used 24,328 collocated ORD-COSMIC pairs in latitudes between 35°S and 35°N for the period from 17 February 2007 to 7 November 2010 (Fig. 7). Prior to the search of collocated (< 2 h and < 200 km) pairs, ORD soundings that are dubious in quality are discarded through a series of data screening as 10 done by Wee and Kuo (2014). Also, COSMIC soundings flagged as bad by CDAAC are dropped.
A challenging factor in GPS RO over the tropics is the sharp refractivity gradient in the lower troposphere, which could cause atmospheric ducting. The variation of phase path in ducting layers can be rapid enough to surpass the receiver's bandwidth. Also, the strong multipath effects pertinent to the anomalous wave propagation are unfavorable to the signal tracking. In addition to increasing random errors, the atmospheric ducting leads to significant negative biases in the retrieved 15 refractivity all the way down to the surface from the top of the ducting layer (Sokolovski, 2003;Ao et al., 2003). On the whole the ducting event occurs rather infrequently and is apt to happen over specific small areas. Once happens, however, it can give rise to exceptionally large errors in the derived refractivity. Consequently, the retrieved refractivity in the heights affected by the ducting tends to appear as outliers in the comparison with ORD. The difficulty is that while the influence of an intense elevated ducting reaches down to the surface, the derived refractivity does not always deviate discernibly large 20 from ORD. Reduced by PHD, the refractivity error might be small in some height below the ducting layer. In some cases, a large error of ORD can counteract the RO error. Moreover, the refractivity error can be the largest at the surface rather than in the vicinity of the ducting layer. This is because AI accumulates the RO error vertically throughout the height range below the ducting layer. All these make the outlier detection futile. The effect of the anomalous wave propagation to AI is very complicated and is beyond the scope of this study. 25 The behavior of retrieved refractivity under the ducting layer is erratic and it makes the comparison with ORD unreliable.
We therefore attempt to detect the occurrence of ducting in each RO event. Once the ducting is identified, we exclude the heights below the ducting layer from the comparison with ORD. We consider this step as a quality control. To do so, we search for ducting layers in the background refractivity ( < −150 km )/ , slightly relaxed from Z ) starting from 7 km and downward. When a ducting layer is detected for the first time, the bending angles below the layer's top are discarded and 30 unused for VR. We do not use the refractivity retrieved via AI for the detection since it cannot have a gradient smaller than Z by definition. Although it is the best not to use the bending angles below the ducting layer for AI, we decided to stick to the CDACC's refractivity and discard the data below the ducting layer instead. In order to reduce the dependency on the Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. background refractivity, the detection is repeated with the refractivity of ORD. Needless to say, both AI and VR are compared to ORD only within the height range that is common to them. In order to make the comparison with ORD robust to outliers, the lower and upper 1 % of AI soundings in the difference from ORD is discarded at each height level. The VR and background refractivity that correspond to the discarded AI are also rejected. Figure 8a shows the difference from ORD in the mean refractivity. The solid black line (denoted AI) indicates the CDAAC's 5 refractivity; the lines denoted as FCST indicate 12h (solid blue) and 24h (dashed green) ECMWF forecasts; and, the refractivity produced by VR is also shown for which the 12h (solid red) and 24h (dashed magenta) forecasts are used as the background. AI results in a distinct negative bias below 3 km, which increases with decreasing height reaching -1.5 % near the surface. This means that the above-mentioned quality control in regard to the ducting layer is not perfect although it has reduced the maximum bias from -3.5 % (not shown). This might be due to the limitation of the forecasts and ORD for the 10 purpose or there could be other RO-side sources of the systematic error. AI shows small positive deviations from ORD above 4.5 km, which are about 0.2 % at 6.5 km. On the other hand, the ECMWF's 12h and 24h forecasts in that height range deviate very little from ORD. In the lowest 2 km, the forecasts show positive systematic deviations that increase with the forecast lead time. Considering that moisture is the dominant contributor to the variability of refractivity in the height range, the ECMWF's forecast model is likely to have a wet bias in the planetary boundary layer, at least with respect to ORD. VR 15 is similar with the background above 2 km, but shifts toward AI in lower heights. In general, VR lies in between AI and FCST. This is not surprising because VR combines RO and FCST by some means, and AI might be a good indicator of RO in the systematic error. Although overly simplistic, the assumption here is that the bias of RO bending angle can be translated into refractivity bias via the observation operator and the subsequent refractivity bias is close to that shown by AI.
As is naturally expected, the relative weighting given to RO and FCST is though crucial for VR to reduce the bias. 20 In the standard deviation from ORD (Fig. 8b), AI is slightly larger than 24h forecast in the lowest 2 km and is comparable elsewhere. As expected, the random forecast error increases with the lead time. Encouragingly, VR is smaller than the background (FCST) throughout the height range regardless of the forecast lead time. This indicates that VR functions as designed attaining an error variance smaller than those of observation (represented by AI) and a priori. Nevertheless, it is noticeable that AI, FCST, and VR are similar in magnitude and structure; and, the difference among them is very small, less 25 than 0.5 % at most. This is because the bulk of the standard deviation is due to the spatiotemporal distance between COSMIC and ORD soundings, in addition to the representativeness error. Although the collocation threshold used in this study is reasonably tight, the significant horizontal inhomogeneity (especially in the moisture) over the tropical region causes the two nearby soundings to differ considerably from each other. The difference is particularly large if the small-scale features in those two soundings are out of phase. The systematic difference on the other hand is insensitive to the distance. 30 The vertical resolution of verifying observation is also relevant to this issue. As mentioned earlier, ORD has a low resolution and the data points are distributed irregularly in height. Moreover, the depths between adjacent data levels differ substantially from one sounding to another depending on the number and location of significant levels. The significant levels are where the observed atmospheric structure turns or changes abruptly. The rapid change around the significant levels, in Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. conjunction with the limited data resolution, results in sizable error if one attempts to interpolate ORD to a regular height grid. Therefore, all other data (AI, FCST, and VR) are interpolated to the heights of ORD in this study. Afterwards, the data points are binned according to the pressure of each ORD height for the statistical comparison. For standard pressure levels, all data samples of the same pressure are grouped together and assigned to the pressure level. In this case, the data samples refer to the refractivity values from all soundings interpolated to the height of the pressure level of collocated ORD. For 5 significant levels on the other hand, bins are allocated in the middle of adjacent standard pressure levels, accommodating all data samples whose pressure are between the neighboring standard levels. For instance, all significant levels between 700 hPa and 850 hPa belong to the same bin. For the reason, the bins have a saw-toothed distribution of data counts as shown in Fig. 8c. The difficulty here is that the samples in a bin are different in their height. Given that the refractivity varies exponentially with height to a good degree, the height discrepancy increases the samples' spread (standard deviation), which 10 could have been reduced greatly if all the samples were taken at the same height. Namely, a portion of the standard deviation shown in Fig. 8b is attributable to the vertical variation of the true refractivity rather than to actual error. A correction accounting for the inter-sample height difference might be possible, but no attempt in the direction is made in this study.

Validation with high vertical resolution radisosonde data
It is possible to overcome the difficulties due to the limited vertical resolution of ORD by using HVRRD. In this 15 comparison, all data (HVRRD, AI, FCST, and VR) are interpolated to a regular height grid. Thanks to the high vertical resolution of HVRRD, the interpolation of HVRRD does not produce large error unless the reported heights are corrupted.
Besides, by placing all data samples to the same height, the uncertainty due to the inter-sample height discrepancy is eliminated. Figure 9 shows the location of 92 HVRRD stations used in this study. The stations are classified into three As mentioned above, all data including HVRRD are interpolated to a regular grid of 50-meter interval and thus the 25 comparison does not need the data binning used for ORD. Without the blurring due to the binning, the comparison yields a more detailed vertical structure of the statistics. Another distinction from the comparison with ORD is that the data rejection below the top of the highest ducting layer detected from radiosonde refractivity is omitted here. The reason is that the differentiation between very shallow layers (due to the high resolution) intensifies measurement noises and so the vertical refractivity gradient computed with HVRRD oscillates rapidly. Unless a heavy pre-filtering is applied, the refractivity 30 gradient is unreliable to be used for the ducting-layer detection. In addition to the strength of refractivity gradient, the effect of ducting also depends on the depth of the ducting layer. That is, some paper-thin ducting layers can be ignored. The Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. radiosonde-based detection of ducting layer is omitted because the choice of pre-filtering, numerical differentiation, and threshold depth is difficult to be made properly and is arbitrary at best. The background-based data detection is still applied.
Figures 10a-c show the systematic differences from HVRRD in the three regions. In the tropics, the overall feature of AI is similar with that shown by the comparison with ORD (Fig. 8a) such as large negative bias near the surface and moderate positive bias in 6-10 km range. Without the spread due to the data binning, the comparison with HVRRD shows larger 5 maxima of the bias in magnitude. The 24h ECMWF forecast (denoted as FCST) is biased negatively below 5 km and positively in 5-10 km range, which is quite different from the comparison with ORD. In addition, VR is analogous to AI deviating significantly from FCST above 2 km, whereas it is close to FCST in the comparison with ORD. After looking into the difference between the comparisons with ORD and HVRRD, it is found that VR gives a heavier weighting to RO observations on average at the HVRRD stations than at the tropical ORD stations. At the tropical HVRRD stations, RO 10 observations are assumed to be more accurate, while the background soundings are less reliable. The difference in the weighting relates to the geological coverage. The HVRRD stations are smaller in the number and located in small specific areas. For instance, about 10 out of 15 HVRRD stations are in latitudes lower than 20°, whereas the majority of ORD sites are located in the subtropics (Fig. 7).
Compared to AI, VR shows a smaller bias in the lowest 2 km and in the stratospheric heights. The error reduction is 15 pronounced in the low tropical troposphere (Fig. 10a) and Arctic stratosphere (Fig. 10c). In general, however, the difference between AI and VR in the bias is small except for in the lowest 2 km. Although all of AI, FCST, and VR deviate significantly from HVRRD, they show similar patterns from mid-troposphere to 20 km in all "latitudinal" regions.
Especially, AI and FCST are very close to each other over the Arctic even though they are largely independent. In principle, AI and FCST can be weakly correlated in the error since the ECMWF assimilates COSMIC data while producing the initial 20 condition of the forecast. If it is the case, however, HVRRD will exert a greater influence on FCST because the radiosonde data, which are also assimilated, are routinely available at the same location (of FCST) at every 12 hours, whereas COSMIC data are obtainable intermittently at random locations. The nearness between AI and FCST thus indicates that their difference from HVRRD is mostly caused by the systematic error of HVRRD. Despite the uncertainty of HVRRD, VR is generally smaller than AI in the systematic difference from HVRRD. 25 The difference among AI, FCST, and VR is more pronounced in the standard deviation (Figs. 10d-f). In the comparison with ORD ( Fig. 8b), AI is no better than (24h) FCST especially in the lowest 2 km. In the standard deviation from HVRRD, On the other hand, AI is smaller than FCST in the tropospheric heights of all regions. In the stratosphere above 13 km over the US and Alaska, however, AI is slightly larger than FCST. It is found through the error estimation described in Sect. 2.3 that the stratospheric degradation of AI is due to a rapid, unsmooth transition of bending angle estimation, which is from the 30 geometrical optics method applied above 20 km to a wave optics method below. Nonetheless, VR agrees better with HVRRD than any of FCST and AI does in all heights regardless of the geographical area. As a result of the high vertical resolution, HVRRD has small-scale vertical variations in the refractivity, which might not be easy to be observed Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. concurrently by other observing systems. Therefore, the evident error reduction attained by VR in both measures (i.e., systematic and random errors) is very impressive.

Summary and concluding remarks
The refractivity soundings provided by GPS RO have been used widely for weather and climate research. Typically, the refractivity is obtained from the inverse Abel transform (Abel inversion) of measured bending angle (measurement). The 5 foremost problem of Abel inversion (AI) is that the error contained in the measurement is carried forward into the refractivity without being controlled. The measurement error includes artifacts introduced by arbitrary noise mitigations that are applied prior to AI. It is challenging to improve the noise mitigation because the separation between signal component and noise is never easy. In order to address the issue, we seek for an approach that accounts for the measurement error while obtaining the refractivity. Our approach is motivated by two facts: Firstly, model has been used to interpret measurements 10 and to evaluate their fidelity in many geophysical applications; and secondly, the error characteristics of RO measurement and the model are fairly well understood at least statistically. After considerable deliberation we come up with an idea that it is possible and synergetic to combine noise attenuation and the inversion together into an estimation problem. Based on the hypothesis, we have proposed a new method for determining atmospheric refractivity from RO bending angle through a variational regularization (VR) of the forward Abel transform. 15 The proposed method considers the numerical inversion of the forward Abel transform. Instead of solving the inverse problem directly, however, VR turns it into an optimization problem in which the solution (refractivity) is required to describe the measurement as closely as possible while the solution's behavior is regularized as per priori knowledge. The optimization problem is solved via the adjoint technique, which is a very efficient way of calculating the gradients of the cost function with respect to all control parameters at once. The fundamental difference between AI and VR is that while the 20 former considers the measurement to be complete, the latter seeks for a filtered solution for which the fitting to the measurement is used as a weak constraint. The proposed method is a generalization of classical regularization methods. The essential feature of VR in the formulation that differentiates the method from classical regularizations is the use of error covariance matrices, which resolves the issues challenging in the classical methods, e.g., as to determining the regularization parameter and as to dealing with the spatiotemporal variation of system parameters. The covariance matrix of measurement 25 error and that of background (a priori estimate of the solution) error play crucial roles in the regularization. The background error covariance matrix, for instance, spreads information from each measurement sample to neighboring locations through the spatial error correlation. As a result, VR cherishes the observational information that is coherent across multiple data samples and prohibits the solution from being too sensitive to the nearest sample. The filtering effect of the error correlation thus encourages VR to produce a smooth solution. It must be emphasized that the smoothness differs from that achievable 30 through low-pass filtering because the solution must attain, in addition to being smooth, the statistical optimality delineated by the error covariance matrices. The optimality accounts for the spatiotemporal variation of the measurement error. Thus, Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. the weightings given to individual measurement samples differ. AI on the other hand treats all measurement samples equally in their trustworthiness.
We have tested the robustness of the proposed method using a synthetic sounding with error. The known true solution in the setting resolves a long-standing problem in real-data tests, which is the ambiguity in the validation stemming from the uncertainty of verifying data. We have demonstrated the limitation of AI with the focus on the difficulty in the pre-filtering 5 of measured bending angle. It is shown that the filtering is particularly challenging for large-amplitude, non-stationary measurement errors. We have described an important source of measurement error in case of model-based inversion or data assimilation, which is the errors-in-variables (EIV) problem. The EIV error is induced by model's error in the coordinate variable (i.e., impact parameter). It is shown that at the heights that the bending angle varies rapidly, the size of EIV error can be exceedingly larger than the typical magnitude of measurement error. We have explained that while the EIV error 10 makes the data assimilation of bending angle less effective, it can be reduced by placing the control variable of the optimization problem in the impact parameter space. Another point highlighted in the test is the posterior height determination (PHD). In both AI and VR, the geometric height of refractivity is determined according to the definition of the impact parameter but only after the refractivity is decided in the impact parameter space. We have described in detail with examples and illustrations that PHD can reduce the refractivity error substantially. We have utilized the synthetic case in 15 order to articulate why we had decided to define the control variable of VR in the observation space, which is to reduce EIV error and to take advantage of PHD.
We have applied the proposed method to actual COSMIC events and validated the result by comparing it with nearby radiosonde soundings. We used two radiosonde data sets, the operationally collected TEMP-format global data and the high vertical resolution data collected at stations operated by NOAA, in order to complement each other's weakness. The former 20 is lower in the resolution but has a superior geographical coverage, and vice versa. The validation shows that VR is suggestively smaller than AI and the background (ECMWF forecasts), in the standard deviation of refractivity from the radiosonde data. Both radiosonde data sets equally show the smaller standard deviation of VR in all heights and latitudes.
VR also agrees better with the radiosonde data than AI does in the mean of refractivity, especially in the lowest 2 km. We have seen in some heights that VR is slightly larger than the background in the mean difference. Although not certain for 25 sure, the likely cause is the systematic error of the radiosonde data.
Based on the results presented herein it is concluded that VR is a considerable improvement over AI in the quality of refractivity. This suggests that VR is able to enhance the data value of RO bending angle with the aid of prior information.
Our study has an important implication for the data assimilation of GPS RO data. Recent years, most of global NWP centers prefer bending angle to refractivity for data assimilation. Although there are good supporting reasons for the preference, this 30 study finds that the assimilation of bending angle has drawbacks that are often disregarded in previous studies. It appears that VR is very promising as it reduces the EIV error and benefits from the PHD. It is straightforward to assimilate the refractivity produced by VR. In order to take the full advantage of VR and to ensure the consistency with the background, however, it will be desirable to incorporate the regularization into data assimilation methods. The recent version of COSMIC Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2017-228 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 26 July 2017 c Author(s) 2017. CC BY 4.0 License. one-dimensional variational retrieval method (1D-Var), for instance, conducts VR on the fly prior to the actual retrieval.
Alternatively, stand-alone VR that shares the background with data assimilation can be carried out as a RO data preprocessing so as to reduce the computational complexity.   is smaller than the true refractivity T; therefore, P is relocated to ~ from the true radius to a higher location ~ since = )/ . This leads the apparent error ~~ to be smaller than the true error . The slope of ~ is constant, -157 km -1 , regardless of the true height or the size of true error as shown by the dashed lines in (a). See text for more details.