Joint 1DVar retrievals of tropospheric temperature and water vapor from Global Navigation Satellite System radio occultation (GNSS-RO) and microwave radiometer observations

Global Navigation Satellite System radio occultation (GNSS-RO) and microwave radiometry (MWR) are two of the most impactful spaceborne remote sensing techniques for numerical weather prediction (NWP). These two techniques provide complementary information about atmospheric temperature and water vapor structure. GNSS-RO provides high vertical resolution measurements with cloud penetration capability, but the temperature and moisture are coupled in the GNSS-RO retrieval process and their separation requires the use of a priori information or auxiliary observations. On the other hand, the MWR measures brightness temperature (T_b) in numerous frequency bands related to the temperature and water vapor structure but is limited by poor vertical resolution (> 2 km) and precipitation.

In this study, we combine these two technologies in an optimal estimation approach, 1D variation method (1DVar), to improve the characterization of the complex thermodynamic structures in the lower troposphere. This study employs both simulated and operational observations. GNSS-RO bending angle and MWR T_b observations are used as inputs to the joint retrieval, where bending can be modeled by an Abel integral and T_b can be modeled by a radiative transfer model (RTM) that takes into account atmospheric absorption, as well as surface reflection and emission. By incorporating the forward operators into the 1DVar method, the strength of both techniques can be combined to bridge individual weaknesses. Applying 1DVar to the data simulated from large eddy simulation (LES) is shown to reduce GNSS-RO temperature and water vapor retrieval biases at the lower troposphere while simultaneously capturing the fine-scale variability that MWR cannot resolve. A sensitivity analysis is also conducted to quantify the impact of the a priori information and error covariance used in different retrieval scenarios. The applicability of 1DVar joint retrieval to the actual GNSS-RO and MWR observations is also demonstrated through combining collocated COSMIC-2 and Suomi-NPP (National Polar-orbiting Partnership) measurements.

1 Introduction

Atmospheric profiles of temperature and water vapor are critical geophysical variables to various weather and climate processes. Humidity in the lower troposphere (LT), especially the planetary boundary layer (PBL), determines the strength and depth of convection and influences circulations and cloudiness, and it indirectly affects atmospheric circulations through non-local influence on infrared irradiances (Stevens et al., 2017). The thermodynamic effect of atmospheric moisture plays an important role in strong positive feedback on global warming and increased subtropical dryness (Sherwood et al., 2010), as well as negative feedback through low cloud formation (Mülmenstädt et al., 2021; Zhou et al., 2016). Also, surface air temperature and atmospheric water vapor content are found in Fujita and Sato (2017) to be connected to extreme precipitation under a warming climate. In particular, Holloway and Neelin (2009) show that the vertical structure of the specific humidity, especially in the free troposphere, is highly correlated to moist convections and rainfalls. To provide accurate measurements of temperature and water vapor, the Global Navigation Satellite System radio occultation (GNSS-RO), passive microwave radiometer (MWR), and hyperspectral infrared (IR) sounders are the key spaceborne sounding techniques for numerical weather prediction and atmospheric science research. In this study, we focus on the combination of GNSS-RO with MWR over IR because of simplicity (lower number of channels) and applicability (better penetration below the thick cloud).

GNSS-RO is a remote sensing technique used to observe the vertical thermodynamic structure from the bending of the occulted GNSS signal ray paths that propagate through the stratified atmosphere (Kursinski et al., 1997). By measuring the time-varying phase of the received signals, the bending angle of each ray path can be computed and inverted to retrieve the vertical profiles of refractivity, which can be used to retrieve temperature and moisture with a priori information. With its limb-sounding geometry, GNSS-RO is capable of providing global coverage with high vertical resolution (∼ 200 m) observations. In addition, the L-band navigation signals used in GNSS-RO can penetrate heavy cloud cover and precipitation and is independent of surface emissivity. Numerous GNSS-RO missions have been launched since 1995, and the number of RO observations has now reached more than 10 000 per day, including the ones from COSMIC-2 (July 2019), METOP-C (November 2018), GRACE Follow-on (May 2018), Sentinel-6 Michael Freilich (November 2020), and commercial CubeSats that have been deployed in recent years.

In addition to GNSS-RO, MWR is also commonly used to profile the atmosphere. MWR utilizes multiple frequency bands, providing information that can be related to the vertical structure of temperature and water vapor in the atmospheric column (Rodgers, 2000). The temperature profile is mostly linked to the oxygen absorption frequency bands located between 50–60 GHz, while the water vapor profile can be retrieved from the brightness temperature measurements of 23.8 and 183 GHz water vapor absorption lines (Liu et al., 2021). Typically, the atmospheric temperature and water vapor profiles are retrieved along with the surface temperature and emissivity simultaneously through an optimal estimation approach (Boukabara et al., 2011). This technology has been extensively utilized in both ground-based and airborne platforms, and numerous MWR instruments have been launched to low earth orbit (LEO) including the Advanced Microwave Sounding Unit (AMSU) and Microwave Humidity Sensor (MHS) on NOAA-18, NOAA-19, and Metop series; Advanced Microwave Scanning Radiometer for Earth Observing System (EOS) (AMSR-E) on Aqua; and Advanced Technology Microwave Sounder (ATMS) on Suomi-NPP (National Polar-orbiting Partnership). Assimilation of observation from AMSU and MHS (Bao et al., 2015), AMSR-E (Kazumori et al., 2008), and ATMS (Bormann et al., 2013) in numerical weather prediction (NWP) models has shown significant positive impacts.

While both techniques provide high quality observations, each has its own limitations. For example, MWR sounding suffers from poor vertical resolution (> 2 km) due to a limited number of frequency channels and a broad range of weighting functions. MWR measurements from different bands could also be affected by clouds, precipitation, aerosol, absorption by ozone and carbon dioxide, and surface properties, which introduce notable uncertainties in the MWR retrievals. Meanwhile, even though GNSS-RO refractivity and bending angle profiles are closely related to the vertical temperature and moisture structures, they cannot be used to retrieve temperature and water vapor independently without the use of auxiliary information, as explained in Sect. 2 below. In general, atmospheric profiles from a global weather analysis, such as National Centers for Environmental Prediction (NCEP) or European Centre for Medium-Range Weather Forecasts (ECMWF), are utilized as a priori states, which could be biased or erroneous and significantly impact the final temperature and moisture products derived from the RO refractivity. In addition, the RO refractivity retrieval itself could be negatively biased due to ducting (Xie et al., 2006; Ao, 2007), which will be translated into considerable positive bias in the temperature and negative bias in the moisture at the top of the boundary layer. While several bias-correction methods for RO within PBL have been proposed (Xie et al., 2006; Wang et al., 2020), they remain challenging to apply in practice, and the corresponding SI (Système International) traceability for refractivity could potentially be lost.

In this study, a one-dimension variational (1DVar) estimation approach is implemented to combine GNSS-RO and MWR measurements and simultaneously retrieve temperature and water vapor profiles that preserve the high vertical resolution from GNSS-RO. The idea of combining RO with MWR or IR soundings is not new: several methods have been proposed to take advantage of these complementary measurements. For example, von Engeln et al. (2001) applied an optimal estimation method to determine the temperature profile from RO bending angle measurements and oxygen line radiance from the passive microwave limb sounder (MLS). In Borbas et al. (2003), the multivariable regression method is used to estimate the atmospheric states based on the coefficients trained by a set of RO refractivity and MWR or IR brightness temperature observations. The simulation results around the tropopause altitude show that the RO–IR combination can improve the temperature and moisture retrievals by 0.5 K and 2.5 %, respectively, compared to retrievals from IR alone. Ho et al. (2007) took a similar regression approach and expanded the application to the lower troposphere to reduce the impacts of RO refractivity retrieval bias by introducing AIRS observations. For both studies, the regression coefficients are calculated from specific training datasets, which could be biased based on its spatial and temporal distribution. The non-linear behavior between the profile variables (temperature and water vapor) and the observables (brightness temperature and refractivity) could also induce errors in the simple linear regression expressions.

The simulation study conducted in Collard and Healy (2003) shows that by combining RO and nadir sounding (MWR or IR) measurements using a sequential 1DVar algorithm, two complementary observations can be fused together and “contribute the greatest impacts to different parts of atmospheric temperature and humidity fields”. Collard and Healy (2003) performed 1DVar twice to include the background temperature and moisture information with RO refractivity and IR brightness temperature observations separately. While this approach works well for unbiased measurements, it is known that in the lower troposphere GNSS-RO refractivity observation could be negatively biased due to ducting which occurs at the top of the PBL (Ao, 2007), phase unwrapping error caused by low signal-to-noise (SNR) conditions (Wang et al., 2016), and possibly small-scaled refractivity fluctuations (Gorbunov et al., 2015). To reduce the refractivity bias due to ducting, Wang et al. (2017) developed an optimal estimation algorithm to choose the unbiased profile from a family of solutions with collocated AMSR-E total precipitable water (TPW) retrievals. However, the Wang et al. (2017) study focuses on correcting the bias due to ducting with TPW retrieval but does not utilize the full information provided by MWR observations.

In this study, we correct the RO bias and improve the retrieval for both RO and MWR with a more optimized approach by (1) using RO bending angle and MWR brightness temperature measurements instead of RO refractivity and MWR TPW retrieval and (2) overcoming the non-linearity of sequential 1DVar by using the RO and MWR observations simultaneously. The GNSS-RO bending angle observations are simulated by the Abel integral, while the MWR brightness temperature observations are simulated by the radiative transfer model (RTM) that considers atmospheric absorption, as well as oceanic surface reflection and emission. Our results show that this approach can simultaneously remove the GNSS-RO temperature and water vapor retrieval bias, reduce the errors in the a priori profiles used in the retrievals, and capture the small-scale vertical structure that MWR cannot resolve. In Sect. 2, the joint 1DVar retrieval algorithm using both GNSS-RO and MWR observations will be described in detail. Simulations, sensitivity tests, and actual data from collocated COSMIC-2 GNSS-RO and ATMS measurements aboard Suomi-NPP will be discussed in Sect. 3. A conclusion will be given in Sect. 4.

2 Joint retrieval algorithm: observations, forward model, and the 1D variational method

2.1 GNSS-RO

In this study, the retrieved 1D bending angle α(a), where a is the impact parameter, is used as the GNSS-RO observation. It can be related to the atmospheric refractivity through the forward Abel equation (Fjeldbo et al., 1971):

$$\begin{array}{}\text{(1)}& \mathit{\alpha }\left(a\right)=-\mathrm{2}a\underset{a}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}x}}{\displaystyle \frac{\mathrm{d}x}{\sqrt{x^{\mathrm{2}}-a^{\mathrm{2}}}}},\end{array}$$

where n is the refractive index, x=rn(r), r is the distance from the center of curvature to a point along the ray path, and a is defined as the x at the location of tangent points. The 1D bending angle is chosen as the observation in this approach for two reasons. First, it is a more “raw” data product than refractivity, reducing the likelihood of additional processing errors. Bending angles are less vertically correlated and less susceptible to errors arising from the lack of spherical symmetry than refractivity. Second, Abel inversion can only be applied to non-ducting atmospheric profiles, where $\mathrm{d}N/\mathrm{d}r>-\mathrm{157}$ N units km⁻¹. Refractivity retrievals of profiles with refractivity gradient less than this critical value will be biased below the ducting layer (Xie et al., 2006); using the bending angle retrieval instead of refractivity can avoid the negative observational bias. However, it is worth stressing that under ducting conditions the refractivity inversion is ill-posed, and infinite solutions can correspond to the same bending angle measurements. Therefore, with the bending angle alone, the thermodynamic states of the atmosphere cannot be determined, and the MWR observation is utilized in our method to provide essential information needed for constraining the refractivity retrieval under ducting conditions.

The refractive index n as a function of height can be calculated from the atmospheric temperature and water vapor profiles (Smith and Weintraub, 1953):

$$\begin{array}{}\text{(2)}& N=\mathrm{77.6}{\displaystyle \frac{p}{T}}+\mathrm{3.73}\times {\mathrm{10}}^{\mathrm{5}}{\displaystyle \frac{e}{T^{\mathrm{2}}}},\end{array}$$

where $N=(n-\mathrm{1})\times {\mathrm{10}}^{\mathrm{6}}$ is the refractivity in N units, p is the pressure in millibar (mbar), e is the vapor pressure (in mbar), and T is the temperature (in kelvin). The atmospheric refractivity can be more accurately calculated by considering non-ideal gas with slightly modified constants as used in Radio Occultation Processing Package (ROPP) (Culverwell et al. , 2015; Burrows et al., 2014). Here we neglect the non-ideal gas effect due to its limited impact (∼ 0.1 % of bending angle) on our analysis results.

In practice, Eq. (1) can be implemented by numerical integration with changing variables. However this is a time-consuming process, especially in 1DVar, where the Jacobian at multiple samples must be calculated over numerous iterations. To reduce the 1DVar Jacobian calculation complexity, Eq. (1) is implemented using the approach taken by ROPP where an exponential refractivity is assumed at each layer:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}\ln \left(n\right)}{\mathrm{d}x}}=\sim {\mathrm{10}}^{-\mathrm{6}}{\displaystyle \frac{\mathrm{d}N}{\mathrm{d}x}}.\end{array}$$

Under this assumption, the bending angle integration shown in Eq. (1) can be simplified. When the refractivity decreases with height as in most situations, the bending angle for each interval is given by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathrm{\Delta }\mathit{\alpha }=& {\mathrm{10}}^{-\mathrm{6}}\sqrt{\mathrm{2}\mathit{\pi }a{k}_i}{N}_i\exp \left(k_i(x_i-a)\right)\\ & \left[\mathrm{erf}\left(\sqrt{k_i(x-a)}\right)\right]{|}_{x_i}^{x_{i+\mathrm{1}}},\end{array}\end{array}$$

where erf is the error function and

$$\begin{array}{}\text{(5)}& k_i={\displaystyle \frac{\ln \left(N_i/N_{i+\mathrm{1}}\right)}{x_{i+\mathrm{1}}-x_i}}.\end{array}$$

When refractivity increases with height, we assume a constant refractivity gradient within the layer, and Eq. (4) is modified to

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }\mathit{\alpha }=-\mathrm{2}\sqrt{\mathrm{2}a}{\mathrm{10}}^{-\mathrm{6}}{\displaystyle \frac{(N_{i+\mathrm{1}}-N_i)}{(x_{i+\mathrm{1}}-x_i)}}\left[\sqrt{(x-a)}\right]{|}_{x_i}^{x_{i+\mathrm{1}}}.\end{array}$$

Currently, ROPP processing removes the portion of the profiles below the ducting layer when its refractivity gradient reaches the critical value. Here we generalize the application of Eq. (6) to the ducting cases. When ducting occurs, x=rn(r) is decreasing within the ducting layer when its corresponding refractivity decreases sharply. In this case, k_i would remain negative and the slope of N(r) can be approximated as constant, which is assumed in the derivation of Eq. (6). Therefore, Eq. (6) is still applicable to calculate the bending angle profile within the ducting layer. In general, this would result in a sharp peak in the bending angle profile that can be observed in occultation profiles.

2.2 MWR

The MWR forward operator we use in this research is based on a radiative transfer model developed at the Jet Propulsion Laboratory (JPL) for simulation, testing, validation, and calibration of microwave radiometer measurements – from instrument design through to on-orbit operation. The forward model is valid for frequencies between 6 and 183 GHz and can simulate both imager and sounder configurations. Oceanic surface emission is based on Meissner and Wentz (2012). For atmosphere absorption, a number of functions can be employed. Oxygen absorption is modeled based on the work of Liebe et al. (1993), Rosenkranz (1998), and Tretyakov et al. (2005). Water (vapor and liquid) absorption is modeled with Liebe et al. (1993). Nitrogen absorption is modeled with Rosenkranz (1998). Given inputs of instrument and environmental parameters, the forward model produces the simulated T_b.

In this article, the T_b measurements of all 22 channels from the ATMS instrument are simulated using the forward operator described above. This process can be improved in the future by using only the channels that are most sensitive to the tropospheric temperature and water vapor structure and discarding the rest. For ATMS, we can focus on channels 4 to 9 (51.76–55.5 GHz) that are most sensitive to the tropospheric temperature and channels 17 to 22 (165.5–183.31 GHz) that are most sensitive to water vapor (Shao et al., 2021). By reducing the number of channels from 22 to 12, the computational complexity and surface property dependence can be significantly decreased. However, the thermodynamic information close to the surface will also be reduced. The trade-off between the number of channels used and the corresponding retrieval accuracy needs to be further investigated in future studies.

2.3 1DVar

Here we generalize the RO 1DVar inversion method described in von Engeln et al. (2003) to the RO and MWR joint retrieval algorithm. As in most current RO processing systems, the temperature, pressure, and the water vapor pressure at each level form the state vector x that is being estimated:

$$\begin{array}{}\text{(7)}& \mathit{x}=\left[\begin{array}{c}T\\ p\\ e\end{array}\right],\end{array}$$

where T is the temperature, p is the pressure, and e is the vapor pressure in 200 m sampling. The vertical range of the state vector spanned from 0 to 10 km altitude. We focus on estimating the lower atmosphere because (1) the contribution of the upper atmosphere to the lower troposphere is small due to exponentially decreases of atmospheric refractivity (with < 0.4 % standard deviation above 10 km) and (2) most of the vapor is distributed in the bottom 10 km of the atmosphere. Each state variable in Eq. (7) is sampled every 200 m over this range to approximately match the vertical resolution of GNSS-RO. The state variables can be connected to the bending angle observation y_ro and brightness temperature observation y_mwr with the forward operators H_ro and H_mwr as explained in Sect. 2.1 and 2.2:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{y}}_{\mathrm{ro}}=\left[\mathit{\alpha }\right]=H_{\mathrm{ro}}\left(\mathit{x}\right),\text{(9)}& {\displaystyle }{\mathit{y}}_{\mathrm{mwr}}=\left[{\mathit{T}}_{\mathrm{b}}\right]=H_{\mathrm{mwr}}\left(\mathit{x}\right),\end{array}$$

where α and T_b are the RO bending angle and MWR brightness temperature observations over the range of impact parameter and 22-channel frequency band, respectively. As an optimal estimation approach, 1DVar will look for the best solution x by minimizing the following cost function:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}J=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-{\mathit{x}}_{\mathrm{b}}\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{b}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{y}-H\left[\mathit{x}\right]\right)}^T{\mathbf{O}}^{-\mathrm{1}}\left(\mathit{y}-H\left[\mathit{x}\right]\right),\end{array}\end{array}$$

where $\mathit{y}=({\mathit{y}}_{\mathrm{ro}},{\mathit{y}}_{\mathrm{mwr}})$ represents the RO and MWR observations, $H\left[x\right]=\left(H_{\mathrm{ro}}\right[x],H_{\mathrm{mwr}}[x\left]\right)$ is the forward operator, and x_b represents the background state variables, or a priori, which can be obtained from a global weather analysis such as NCEP, MERRA-2, or ECMWF. B is the error covariance matrix of the a priori x_b. O is the error covariance matrix of y which includes observational and model representation errors. In this research, O is given in a diagonal form with estimated uncertainties of bending angle of ${\mathit{\sigma }}_{\mathit{\alpha }}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$ rad and brightness temperature of ${\mathit{\sigma }}_{T_{\mathrm{b}}}=\mathrm{0.25}$ K. The chosen bending angle uncertainty is comparable to the one of COSMIC-2 at ∼5 km impact height (Todling et al., 2022), which represents the average RO observation uncertainty over the lower troposphere. Actual RO bending uncertainty can be much smaller in the free troposphere (> 4 km) and provides better quality observations (Todling et al., 2022). Therefore, the simplified constant uncertainty setting over impact parameters shows the worst case of combination above the PBL. For simplicity, the state covariance matrix B also has a diagonal form that incorporates the uncertainty of the background state variables (T: 2.5 K, p: 1 %, e: 40 %). The background uncertainties of T, p, and e are chosen according to the ECMWF model used in von Engeln et al. (2003), and a similar amount of uncertainties are also assumed for other weather models. The use of a diagonal form for both covariance matrices implies the independence of measurement errors with respect to height. In reality this is not entirely true, and a sensitivity test with respect to the off-diagonal terms is given in the next section. The uncertainty defined in B and O will affect the optimal solution of temperature and vapor retrievals. More discussions on this will be provided in Sect. 3.1.

The estimated state variable x can be determined by iteratively solving the following formula:

$$\begin{array}{}\text{(11)}& {\mathit{x}}_{n+\mathrm{1}}={\mathit{x}}_{\mathrm{b}}+{\mathbf{G}}_n\left[\left(\mathit{y}-H\left[{\mathit{x}}_n\right]\right)-{\mathbf{K}}_n\left({\mathit{x}}_{\mathrm{b}}-{\mathit{x}}_n\right)\right],\end{array}$$

where x_n is the state variable at the nth iteration, and x_b is the state a priori. K_n is the state Jacobian matrix which can be calculated by perturbing each individual variable in x:

$$\begin{array}{}\text{(12)}& {\mathbf{K}}_n={\displaystyle \frac{\partial H\left(\mathit{x}\right)}{\partial \mathit{x}}}{|}_{\mathit{x}={\mathit{x}}_n}.\end{array}$$

In this work, we perturb T and p by 0.1 K and 0.1 hPa, respectively, and we perturb water vapor by 0.02 hPa. Alternatively, one can use the Jacobian term calculated by the ROPP software to lower the time consumption for numerical differentiation computation. The gain matrix G_n can then be calculated from

$$\begin{array}{}\text{(13)}& {\mathbf{G}}_n={\left({\mathbf{B}}^{-\mathrm{1}}+{\mathbf{K}}_n^T{\mathbf{O}}^{-\mathrm{1}}{\mathbf{K}}_n\right)}^{-\mathrm{1}}{\mathbf{K}}_n^T{\mathbf{O}}^{-\mathrm{1}}.\end{array}$$

Here we set the conditions $T_n-T_{n-\mathrm{1}}<\mathrm{0.1}$ K and $e_n-e_{n-\mathrm{1}}<\mathrm{0.2}$ hPa as the convergence criteria.

3 Results and analysis of the joint retrieval algorithm

3.1 Large eddy simulation

An important goal of this study is to investigate the effectiveness of the joint retrieval on the planetary boundary layer (PBL). To this end, high-resolution profiles obtained from the large eddy simulation (LES) (Kurowski et al., 2023) are used to validate the 1DVar algorithm. Temperature, pressure, and water vapor from LES in three different physical regimes are extracted to simulate the RO bending angle and MWR T_b observations and regarded as the truth:

1. shallow to deep convective cloud transition in the Amazonia region, i.e., the Large-scale Biosphere-Atmosphere (LBA) experiment ​​​​​​​(Gustavo Gonçalves de Gonçalves et al., 2013; Grabowski et al., 2006);

2. shallow cumulus cloud in the west Atlantic Ocean, i.e., the Rain In shallow Cumulus over the Ocean (RICO) campaign (Rauber et al., 2007);

3. marine stratocumulus region in the northeast Pacific, i.e., the DYnamics and Chemistry Of Marine Stratocumulus (DYCOMS) campaign (Stevens et al., 2003).

For each case, the RO bending angle and MWR T_b observations are simulated using the forward model described in Sect. 2.1 and 2.2. Here we simulate the bending angle over the 0–30 km range of impact height with 50 m sampling to match the RO bending resolution at the lower troposphere. In this paper, no thermal noise is added to either simulated observation, although the 1DVar processing takes this uncertainty into account with the observation error covariance matrix O (Sect. 2.3). The temperature a priori is derived by smoothing the true temperature profile from LES with a 1 km moving window. In addition, we added a constant −2 K bias to the a priori temperature at every altitude to simulate a potential bias in the background profile and test if the 1DVar algorithm can remove it. The total pressure and water vapor pressure a priori are also defined by smoothing the true profile with 3 and 1 km boxcar filters, respectively. This would remove most of the small-scaled structures in the lower troposphere that typically do not show up in the background profiles. For each case, the implemented 1DVar algorithm is used in three different scenarios: RO only, MWR only, and a joint retrieval using RO and MWR (RO+MWR).

Figure 1 The 1DVar results for LES temperature (a, b) and water vapor (c, d) profiles in the LBA scheme. Note that the RO temperature solution (green) is close to the a priori profile (orange); the two curves almost overlap in (a).

The 1DVar results between the surface and 5 km are shown in Figs. 1 to 3, corresponding to the three campaigns listed above, respectively. Each figure includes the temperature profiles (a), temperature difference from truth (b), water vapor pressure (c), and water vapor difference from truth (d). In the case of LBA (Fig. 1), the true T and e are relatively smooth with little fine-scale vertical structures; the background profiles are close to the true profiles except for the −2 K added to the a priori T. Here we intentionally added −2 K to test the robustness of the joint retrieval under a large amount of background uncertainty. While this magnitude of error in the model is not expected to be common, it could occur at the top of the PBL where a large temperature transition occurred over a short altitude interval.

Figure 1b shows that the temperature solution of the RO-only scenario mostly follows the a priori profile with the −2 K error because the 1DVar is heavily weighted with respect to the T a priori compared to the water vapor a priori. The resulting RO-only water vapor retrieval shown in Fig. 1d is also slightly biased (−1 hPa) near the surface to compensate for the negative T bias and yield the same refractivity. These biases are representative of the sensitivity of the derived T and e to the T a priori errors in a RO 1DVar-retrieval in the lower troposphere. On the other hand, the MWR-only solutions (purple dotted-dashed) appear to be less sensitive to the a priori T error with not even 1 K difference at maximum compared to the truth. The results imply that, unlike RO, the MWR is capable of independently solving for T and e, except for the small-scale structure within 500 m altitude from the surface. Combining RO bending and MWR T_b in the 1DVar framework discussed above, we observe that the retrieved T data (red solid line) in Fig. 1a and b are close to truth, despite the −2 K bias that was added to the a priori, and generate a detailed water vapor retrieval (red solid line) in Fig. 1c and d that is more accurate than either MWR (purple dotted-dashed) or RO (green dotted-dashed) alone.

Figure 2 Same as Fig. 1 but with the RICO case.

In the RICO case (Fig. 2), the water vapor retrieval from the MWR-only scenario (purple dotted-dashed) shows a large error of 2 hPa at 2 km. This is due to the low vertical resolution of the MWR measurements which miss the fine-scale structure below 2 km. As a result, the 1DVar solution for MWR-only shares nearly identical structure with the given a priori profile (orange dotted), which was heavily smoothed with small structure removed and correct only the bias from the T measurements. By contrast, the RO-only solution (green dashed) is able to reduce the water vapor retrieval error relative to the a priori (to within 1 hPa) and better resolve the fine vertical structure. However, as shown in the LBA case, the negative e bias (0.5 hPa at the surface) persists in the lower troposphere due to a −2 K bias in the temperature a priori. By combining both RO and MWR with the proposed 1DVar approach, the T and e biases are corrected for and the small scale structures in water vapor are captured.

Figure 3 Same as Fig. 1 but with the DYCOMS case.

Figure 3 shows the T and e profiles in the DYCOMS case, which is located at the stratocumulus region in the northeast Pacific. The lower troposphere in this region is known to have a sharp transition at the top of the boundary layer as can be seen from the temperature and moisture profiles in Fig. 3a and c at the ∼ 1 km height. The sharp transition of temperature and moisture create a ducting layer (Xie et al., 2006; Ao, 2007) where the RO tangent point cannot be located and the thermodynamic information in the ducting layer is lost. This results in an ill-posed inversion problem where multiple refractivity solutions would correspond to the same bending angle profile, and the standard Abel inversion would resort to a solution without ducting and cause a negative bias of up to 15 % in refractivity inside the layer. This is the reason why bending angle is used in the 1DVar instead of refractivity and could, potentially, avoid this refractivity bias. Nonetheless, without additional information, the solution may not converge to the correct refractivity profile in the family of solutions. As shown in Fig. 3c and d, the water vapor retrieval in the RO-only scenario still contains a large negative bias (−2.5 hPa close to the surface) compared to the true profile. This can be due to the −2 K bias introduced in the a priori T.

To overcome the negative N bias due to ducting, Wang et al. (2017, 2020) used the information of collocated MWR TPW retrievals and RO grazing reflected bending angles as constraints to choose the unbiased refractivity profile. Here, however, we show that combining MWR T_b with the RO bending angles results in a reduced negative bias due to ducting. As shown in Fig. 3c and d, the MWR-only and RO+MWR scenarios are no longer limited by the ducting condition and are able to correct the moisture bias within the boundary layer. The low resolution of the MWR measurement misses the sharp transition at the PBL top, but the combination of RO and MWR preserve the advantages of both measurements. For temperature retrieval (Fig. 3a and b), while the MWR and RO+MWR solutions are able to correct the −2 K bias, all three solutions fail to resolve the sharp change at the transition due to strong smoothing of the background profile. With a smaller covariance of T a priori, the 1DVar temperature solution will depend more on the low-resolution a priori and MWR observations rather than the RO measurements; therefore, the large error at the transition layer should not come as a surprise.

3.2 Sensitivity study

To investigate the sensitivity of the 1DVar solution to the a priori and the measurement covariances, we perform a number of simulations using a radiosonde profile from the MAGIC campaign (file name: magsondewnpnM1.b1.20121104.120900) (Keeler and Burk, 2012; Lewis , 2016). In this study, we chose a profile that does not contain any ducting conditions to avoid errors due to super-refraction. The temperature, water vapor, and the corresponding refractivity and bending angle profiles are shown in Fig. 4. The temperature and water vapor priors are obtained by a 1 km boxcar filter to the radiosonde observation (RAOB) T and e profiles. In addition to smoothing, constant biases in T of −2, −1, 0, 1, or 2 K were added to the a priori T profile. Alternatively, constant biases in e of −40 %, −20 %, 0 %, 20 %, and 40 % were added to the a priori e profile. For each of these cases, three different scenarios corresponding to RO only, MWR only, and combined RO+MWR were run. The resulting T and e solutions are then compared to the truth by computing the root mean square error (RMSE) over the 0–5 km range. The results are shown in Fig. 5. It is clear that the RMSE for the RO-only solution (green dashed) largely follows the a priori RMSE (orange dotted). This is expected since no independent T information is provided, and the 1DVar has to rely on the a priori to estimate the T profile. By contrast, the MWR-only RMSE (purple dotted dashed) is considerably lower (∼ 0.5 K compared to the ∼ 1.8 K) and is less sensitive to the T bias. The combined RO+MWR scenario (red solid) reduces the T RMSE even further to 0.4 K and is nearly independent of the T bias. On the other hand, the e RMSE of the MWR-only scenario shown in Fig. 5b mainly follows the a priori due to its low resolution. While the RO-only e solution has lower RMSE, it is sensitive to the T bias because temperature and water vapor are coupled. The combination of RO and MWR is able to uncouple the T and e and make the RMSE almost always below 0.3 hPa, independent of the T bias. By contrast, the water vapor bias has a smaller impact on retrieved temperature (Fig. 5c) or water vapor (Fig. 5d). This is mainly due to the fact that the water vapor standard deviation (set to 40 %) is much larger than the temperature standard deviation (2.5 K), causing the 1DVar solutions to match better with the temperature a priori rather than the moisture a priori. Hence the MWR-only moisture solution is not strongly biased, although its structure mainly follows the water vapor a priori. Overall, the RO+MWR combined solutions show the lowest sensitivity to the temperature and water vapor a priori errors.

Figure 4 (a) Temperature, water vapor, and the refractivity profiles calculated using Eq. (2); (b) the bending angle from the forward Abel as described in Sect. 2.1.

Figure 5 Sensitivity test of the (a, c) temperature and (b, d) vapor retrieval with respect to the temperature (a, b) and vapor pressure (c, d) a priori.

Figure 6 Sensitivity test of the (a, c) temperature and (b, d) vapor retrieval with respect to the T_b covariance (a, b) and the bending angle covariance (c, d).

Now we turn our attention to the sensitivity of the retrievals to the MWR and RO measurement error covariances. The results are shown in Fig. 6 for five different levels of uncertainty in T_b (0.25, 0.50, 1.00, 1.50, 2.00 K) and bending angle (0.0008, 0.0012, 0.0016, 0.002, 0.0024 rad). In all cases, the combination of RO and MWR observations have a positive impact on the solution, causing the RMSE to reduce. As expected, a larger covariance reduces the value of measurements, leading to larger retrieval RMSE values. This positive, near-linear trend can be seen in both MWR T_b (Fig. 6a, b) and RO bending angle measurements (Fig. 6c, d). In this test, no noise is added to the measurements, so the relationship is only driven by the covariance values. It can be observed that the temperature RMSE (Fig. 6a, c) of the RO+MWR solution (red solid lines) increased by 0.22 K when T_b covariance increases from 0.25 to 2 K, while it only adds 0.05 K and is saturated when bending angle covariance increased by 3 times. This result illustrates the higher-temperature information content from MWR T_b relative to RO bending angle between 0 to 5 km as expected from previous studies (e.g., Collard and Healy, 2003). Note that T_b covariance should have no effects on estimates for the RO-only scenario. Similarly, error covariance in bending should have no effects on estimates for the MWR-only scenario, as confirmed in Fig. 6c and d. The MWR/RO combined solution has the lowest temperature RMSE value in both scenarios (Fig. 6a, c), while its moisture RMSEs are almost the same as the one of RO-only solutions (Fig. 6b, d). As expected, the results show that when the a priori is not biased, the MWR observation provides little information content to the estimation process.

Figure 7 Sensitivity test of the (a, c) temperature and (b, d) vapor retrieval with respect to the temperature (a, b) and vapor pressure (c, d) covariance.

The third sensitivity test shown in Fig. 7 quantifies the retrieval RMSE with respect to a change of temperature (1.5, 2.5, 3.5, 4.5, 5.5 K) and moisture (20 %, 40 %, 60 %, 80 %, 100 %) in the a priori covariance. A higher a priori covariance implies a greater reliance on observations in the 1DVar estimation process and reflects different behaviors for different measurements. For the RO-only case (green curve in Fig. 7a), the temperature RMSE decreases under larger temperature a priori covariance. This may be due to a better resolved temperature profile, especially at the top of the PBL, where the strong constraint from the a priori is relaxed. On the other hand, the MWR-only case (purple curve in Fig. 7a) shows the opposite trend, indicating that most of the vertical information of the MWR-only retrieval is derived from the temperature a priori. In fact, when the temperature a priori covariance is larger than 4.5 K, its corresponding RMSE exceeds the one from the a priori itself (orange). However, the large variance observed in the temperature RMSE does not propagate to the water vapor RMSE (Fig. 7b) and shows a relatively stable curve when the temperature a priori covariance changes.

It is worth noting that the perturbation of the water vapor a priori covariance does not cause a dramatic change in water vapor results (Fig. 7d). At the same time, higher water vapor a priori covariance would push the RO-only temperature solution closer to the a priori profile, leading to increasing temperature RMSE as shown in Fig. 7c. By contrast, the MWR-only temperature solutions show a negative trend which implies the independence of MWR temperature and water vapor retrieval based on the T_b measurements. As the results show in Figs. 5 and 6, the RO+MWR joint retrieval solutions remain the lowest RMSEs among all scenarios in this sensitivity test, with errors at ∼ 0.4 K (temperature) and ∼ 0.3 hPa (vapor pressure) in all levels. These results show that 1DVar retrievals from combining RO and MWR can better reduce the error propagated from the a priori, and its solutions maintain low sensitivity to observation and a priori covariance.

Figure 8 Sensitivity test of the (a, c) temperature and (b, d) vapor retrieval with respect to the off-diagonal term in the background covariance matrix. The diagonal matrix is used for the first two panels (a, b) and the off-diagonal elements are added in the last two panels (c, d).

So far, we used the diagonal covariance matrix by assuming the model data are independent for all levels, but in reality the temperature and water vapor at neighboring levels may be correlated. To investigate the impact of the off-diagonal elements we repeat the RAOB simulation using the covariance matrix computed by the equation documented in Healy (2001):

$$\begin{array}{}\text{(14)}& {\mathbf{B}}_{i,j}={\mathit{\sigma }}_i{\mathit{\sigma }}_j\exp \left[-{\displaystyle \frac{{\left(z_j-z_i\right)}^{\mathrm{2}}}{l^{\mathrm{2}}}}\right],\end{array}$$

where i and j are the column and row indexes of the background matrix, respectively; σ_i is the standard deviation at the ith level; z is the level altitude; and l is the scale length. In this test we set l as 0.75 km, and the corresponding temperature and water vapor results are shown in Fig. 8. As the figure shows, the off-diagonal elements will smooth the estimated solutions with the correlation between different levels, which is consistent with the results shown in Healy (2001). The smoother profile implies reduced information content in each sample when correlations between the neighboring levels are increased. The enforcement of smoothness in the estimated profiles makes the combination with sharp RO bending angle observations in the lower troposphere difficult. The water vapor profile shows less small-scale structures, and the RO and RO+MWR temperature profile has larger error above the PBL. Therefore, while the estimation results are generally insensitive to the covariance, they could be sensitive to the correlation in the background, and the off-diagonal terms have to be carefully chosen in practice. The method of selecting optimal off-diagonal terms in the background matrix needs to be further investigated in the future.

3.3 Real ATMS and COSMIC-2 data

We applied our joint RO+MWR 1DVar algorithm to Suomi-NPP and COSMIC-2 measurements to assess the applicability of this method to real measurements. Fig. 9 shows the collocated cases between Suomi-NPP and COSMIC-2 in a 6 h period starting on 1 April 2019, 00:00–06:00 UTC. The Suomi-NPP data we used are the L1B calibrated/bias-corrected brightness temperature product (Lambrigtsen, 2018) provided by the Goddard Earth Sciences Data and Information Services Center (GES DISC). This radiance dataset has been calibrated by the in-flight ATMS antenna/receiver systems that measure the radiation from two calibration sources during every scan cycle. One is the cosmic background radiation from space (cold space) and the other is the internal blackbody calibration target (hot target). By taking the radiometric counts measurements from both sources and combining them with instrument errors that have been accurately modeled from ground thermal-vacuum tests, the published T_b measurements are calibrated and bias-corrected. The details of the in-flight calibration process are documented in the “Algorithm Theoretical Basis Document” accompanying the published dataset (Lambrigtsen, 2018). Here we assume the T_b observations are unbiased and can be directly used for 1DVar combination without additional calibration steps.

Figure 9 Collocated cases between Suomi-NPP and COSMIC-2 on 1 April 2019, 00:00–06:00 UTC. The red dots are the tangent point location of each COSMIC-2 RO at its lowest penetration height, while the background color is the (a) 31.4 GHz and (b) 165.5 GHz brightness temperature measurement from the onboard ATMS instrument.

Due to the lower inclination of COSMIC-2 orbits, all of the collocated cases are located between −45 to 45^∘ in latitude. To illustrate the 1DVar algorithm, we chose one case (2020-04-01-03:10c2f4_gps58) in the Atlantic Ocean (17.91^∘ N, 36.75^∘ W). This case is chosen because it is located over the ocean, where surface emissivity can be robustly modeled for the MWR observation. In addition, the collocated MERRA-2 reanalysis shows low ice and liquid water contents, which suggests that the MWR observation is less likely to be impacted by clouds and precipitation. Furthermore, the RO retrieval penetrated sufficiently deep in this case so that the 1DVar can estimate the thermodynamic structure within the PBL. The estimation results using different instrument scenarios are shown in Fig. 10.

Figure 10 An example (2020-04-01-03:10c2f4_gps58) of the actual collocated RO-MWR combination (a) temperature and (b) water vapor retrieval.

Since there is no ground truth in this comparison, we cannot conclude definitely which method provides the best solution. However, by comparing the vertical structures between the different estimates, we can find several clues regarding their information content.

The orange dotted lines shown in Fig. 10a and b are the profiles from the NCEP analysis, used as priors in the 1DVar processing. Note that in this case the COSMIC-2 data have not yet been assimilated in the NCEP analysis (NCEP, 2023), and no RO information from COSMIC-2 is included in the a priori. The green, red, and purple lines are the 1DVar solutions in the RO-only, RO+MWR, and MWR-only scenarios, respectively. The yellow line is the Community Long-term Infrared Microwave Combined Atmospheric Product System (CLIMCAPS) retrieval (i.e., not 1DVar) using the ATMS MWR and Cross-track Infrared Sounder (CrIS) data (Smith and Barnet, 2019), provided by GES DISC (Barnet, 2019). Here, MWR provides additional temperature information as shown in Fig. 10a so that the MWR-only solution (purple) deviates from the a priori (orange) and RO-only (green) temperature retrievals by ∼ 3 K below 2 km. The joint RO+MWR retrieval lies in between the RO-only and MWR-only solutions and has approximately the same temperature solution as the CLIMCAPS solution. On the other hand, the water vapor profiles shown in Fig. 10b demonstrates that the joint RO+MWR solution is able to resolve the small-scale moisture structure throughout the profile. The deviation between the RO-only and the RO+MWR moisture solutions is caused by the temperature information that MWR provides, which in fact matches better with the ATMS moisture profile below 1 km. While the improvements cannot be quantified without ground truth, the results show that the 1DVar joint retrieval combines the information and strengths of both techniques and is able to provide the high-resolution and low-bias thermodynamic profiles that a single technique cannot retrieve.

The effectiveness of combining RO and MWR observations within the framework of the proposed 1DVar approach is limited by several factors. First, the existence of cloud and precipitation can significantly increase the forward model error in the T_b calculation (Errico et al., 2007). Second, the MWR input and forward modeling error could still exist even in clear-sky events due to other factors that affect the surface emissivity, such as surface type, surface temperature, and surface wind speed. In this study, we limited the application of joint retrieval to an RO over the ocean. In addition, a quality control (QC) test based on the difference between the observed T_b and MERRA-2 calculated T_b for each channel was applied to ensure that the MWR measurements are not biased due to the reasons stated above. Statistically, 73 out of 132 collocations found within a 6 h period passed the QC test when an rms threshold of 10 K is applied on all 22 ATMS channels. Alternatively, one can also implement an additional calibration process to further remove these factors from the MWR data and forward modeling errors. We expect this approach will improve the quality and quantity of data available for the joint retrievals, but this requires more studies in the future to statistically validate its effectiveness and uncertainty.

4 Conclusions and discussions

In this article, we described a 1DVar approach combining two complementary measurements to obtain high vertical resolution data and solve for temperature and moisture simultaneously. Simulations were performed where high-resolution LES profiles from three different campaigns were used as truth, and three different scenarios corresponding to RO-only, MWR-only, and the RO+MWR combination were tested. The results show that potential biases in the a priori information used in the 1DVar can be significantly reduced after adding T_b observations from MWR. At the same time, the high-resolution RO bending angle observation provides the needed vertical moisture information. The complicated thermodynamic structure in the lower troposphere, including the one with ducting, can therefore be better resolved with much smaller biases compared to the ones using RO or MWR alone. We also analyzed the sensitivity of the temperature and vapor retrieval in each scenario to the a priori, background covariance, and observation covariance and showed that the RO+MWR combination is the most stable among the three scenarios when the background vertical levels are assumed uncorrelated. Finally, the 1DVar approach is applied to real data from COSMIC-2 and Suomi-NPP observations, and the results show the promise of the 1DVar technique.

The joint retrieval approach is similar to the optimal estimation currently used in data assimilation. Both techniques minimize the cost function, require covariance matrices to define the background and observations uncertainties, and use the forward operator to map the state vector to the observables. Typically, data assimilation for numerical weather prediction (NWP) purposes utilizes a 3DVar or 4DVar approach to account for horizontal and temporal coverage. Here we only estimate the state variable on a single spatial dimension (altitude), assuming a spherically symmetric atmosphere, which reduces the number of the state variables and allows us to better quantify the relative strengths and weaknesses of the information used in the retrieval process. This also allows us to introduce as many vertical levels as needed to capture the highest vertical resolution possible with RO measurements. The proposed 1DVar joint retrieval is less dependent on any given operational NWP model and is not limited by observation QC criteria, a priori contribution, or the vertical grid resolution applied in operational NWP models. Therefore, the jointly retrieved temperature and water vapor profiles can be good candidates for validating weather and climate models. In addition, under the circumstances when individual RO and/or MWR measurements are not included in the data assimilation process due to the internal QC of the NWP systems, the joint retrieval profiles can potentially provide additional data that are more amenable for NWP processing.

Two future improvements are envisioned to make the proposed 1DVar approach more practical and accurate. First, all the observation error covariance matrices we used in this study are diagonal. This implies that errors of T_b of different channels and RO bending angle measurements at different heights are assumed to be independent, which is an assumption that is not perfectly valid. The vertical smoothing in the RO bending angle profile could lead to a high correlation of two or more neighboring samples, especially in lower altitudes. One possible way to address this issue is to use the Desroziers diagnosis (Desroziers and Ivanov, 2001) to examine if the error covariance matrices used are reasonable, which is a topic that will be investigated in future studies. Second, the 1DVar approach combines RO measurements with a single set of MWR T_b measurements that are nearest to the RO location. The RO location is defined as the latitude and longitude of the tangent point of the lowest link in the RO. However, in reality, RO has an extended footprint of several hundred kilometers in the occultation plan and a potential drift of the tangent point out of the occultation plane by as much as tens of kilometers. Taking horizontal variability into account could better represent the MWR T_b observation at the tangent point location for each corresponding altitude. This requires a better understanding of the Jacobian function as a function of height for each channel and the location of the tangent point prior to the joint retrieval process.