A hybrid ice hydrometeor retrieval algorithm for submillimeter-wave radiometers in support of the PolSIR and PMM missions

This paper presents an ice hydrometeor retrieval algorithm for submillimeter-wave radiometers in support of the upcoming PolSIR (Polarized Submillimeter Ice Cloud Radiometer) and PMM (Precipitation Measuring Mission) satellite missions. The algorithm employs a hybrid Bayesian Monte Carlo Integration (BMCI) and Optimal Estimation Method (OEM) approach, which leverages the strengths of the BMCI method but extends the retrieval capability beyond the limitations of the a priori database when BMCI alone fails to identify enough database cases matching the observations. To address the highly non-Gaussian nature of the a priori statistics, a method using cumulative distribution functions (CDFs) and empirical orthogonal functions (EOFs) is applied to enable effective implementation of the OEM algorithm. With the CDFs/EOFs method, the OEM can maximize the posterior probability density function using the a priori constraint that is largely consistent with that used in the BMCI step, ensuring that the entire retrieval operates under a nearly uniform prior constraint.

The algorithm is first applied to evaluate the PMM-C²OMODO (Convective Core Observations through MicrOwave Derivatives in the trOpics) radiometer using simulated observations. Retrieval accuracies for key microphysical parameters are presented. Also, the retrieval diagnostics, including the vertical resolution, Degrees of Freedom (DoF), and Shannon information content, are analyzed. The algorithm is further applied to real CoSSIR (Configurable Scanning Submillimeter-wave Instrument/Radiometer) observations during the IMPACTS (Investigation of Microphysics and Precipitation for Atlantic Coast-Threatening Snowstorms) campaign, and the results are evaluated against the collocated triple-frequency radar retrievals from CRS (Cloud Radar System) and HIWRAP (High-altitude Imaging Wind & Rain Airborne Profiler) observations. Multiple ice particle habits are examined, and both layer-resolved and column-integrated mass and size parameters are evaluated. Both simulated experiments and real-observation retrievals demonstrate that the hybrid BMCI-OEM method is highly effective in reducing simulated and observed brightness temperature (TB) discrepancies. With an appropriate selection of ice cloud particle habits, the TB simulations closely reproduce the observations. Using the logarithmic difference as the quantitative metric, the CoSSIR-retrieved ice water content, layer-resolved particle diameter, ice water path, and column-averaged mean mass diameter with the hybrid BMCI-OEM algorithm differ from triple-frequency radar retrievals by 3.75, 0.82, 2.53, and 0.55 dB, respectively, representing reductions of 0.21, 0.02, 0.85, and 0.13 dB relative to BMCI-only retrievals. The hybrid Bayesian framework also demonstrates high extensibility to other remote-sensing observations. As more information becomes available through multi-sensor integration or the use of hyperspectral measurements, the hybrid Bayesian algorithm shows increasing potential to better constrain cloud microphysical properties.

1 Introduction

Submillimeter wavelengths comprise an important spectral regime for remotely sensing ice cloud properties (Buehler et al., 2007, 2012; Wu et al., 2024). Compared with traditional microwave frequencies, high-frequency radiometer channels have wavelengths comparable to ice particle dimensions, which makes them more responsive to ice particle scattering. Meanwhile, submillimeter waves have greater penetration depth through thick ice cloud compared with visible (VIS) or infrared (IR) bands. The submillimeter spectrum is therefore critical for bridging the gap in cloud and precipitation remote sensing. High-frequency radiometry for ice-cloud remote sensing has been validated through a series of field campaigns using airborne radiometers including CoSSIR (Compact Scanning Submillimeter Imaging Radiometer, later renamed the Configurable Scanning Submillimeter-wave Instrument/Radiometer, Evans et al., 2005; Adams et al., 2020) and ISMAR (International Submillimetre Airborne Radiometer, Fox et al., 2017) as well as several pathfinder spaceborne missions such as IceCube and the Arctic Weather Satellite (AWS). Notable suborbital campaigns include CRYSTAL-FACE (Cirrus Regional Study of Tropical Anvils and Cirrus Layers – Florida Area Cirrus Experiment, Evans et al., 2005), TC4 (Tropical Composition, Cloud and Climate Coupling, Evans et al., 2012), COSMICS (Cold-air Outbreak and sub- Millimetre Ice Cloud Study, Brath et al., 2018), and IMPACTS (Investigation of Microphysics and Precipitation for Atlantic Coast-Threatening Snowstorms, McMurdie et al., 2022; Yorks et al., 2025). Launched in 2017, IceCube was the first experimental spaceborne submillimeter radiometer operating at 883 GHz. IceCube has demonstrated its anticipated scientific value by capturing a wide dynamic range of ice water path (IWP) and showing sensitivity to ice particle size, while also validating the space technology readiness for submillimeter radiometers (Gong et al., 2021). In August 2024, the European Space Agency (ESA) launched the pathfinder operational satellite AWS into polar orbit, and it has been continuously collecting high-quality data ever since (Eriksson et al., 2025). In the coming years, several submillimeter-wave radiometers are expected to be deployed as part of different satellite missions. Key instruments include ESA’s ICI (Ice Cloud Imager, Eriksson et al., 2020), NASA’s PolSIR (Polarized Submillimeter Ice-cloud Radiometer), and the C²OMODO (Convective Core Observations through MicrOwave Derivatives in the trOpics, Brogniez et al., 2022; Auguste and Chaboureau, 2022) planned for the JAXA (Japan Aerospace Exploration Agency) and NASA PMM (Precipitation Measuring Mission, Iki et al., 2025). These high-frequency radiometers represent major steps toward achieving global ice cloud observations from space.

The retrieval algorithms applied in high-frequency microwave ice cloud studies are fundamentally based on Bayes’ theorem, which specifies knowledge in terms of probability density function (PDF) and combines prior information with observations to constrain the state variables. Among various algorithmic approaches, the Bayesian Monte Carlo Integration (BMCI) (Evans et al., 2002, 2005) and Neural Network (NN) (Wang et al., 2017; Brath et al., 2018) are two representative implementations of the Bayesian framework. Among NN-based methods, the Quantile Regression Neural Network (QRNN) (Pfreundschuh et al., 2018; May et al., 2024; May and Eriksson, 2025) is particularly notable for its ability to provide estimates of retrieval uncertainty together with the retrieval results. Both BMCI and NN algorithms depend on a precalculated retrieval database, which consists of random atmospheric and cloud samples distributed according to the a priori PDF and the corresponding simulated brightness temperatures (TBs). The BMCI method assigns weights to the random samples based on the likelihood functions and integrates over the posterior PDF to obtain the mean state and the standard deviation. The NN approach uses the a priori database as a training dataset to learn the statistical distribution to predict state variables for new observations. Both algorithms have been thoroughly evaluated using simulated and real airborne observations, and they are planned to be operationally implemented in the upcoming ICI mission (Eriksson et al., 2020; May et al., 2024).

The primary limitation of these database-driven methods is the sparse coverage of the measurement space by database cases. Theoretically, both BMCI and QRNN require adequate database samples that match the input TBs within the measurement uncertainty. In practice, however, even a large retrieval database inevitably has regimes in the multidimensional measurement space that are sparsely populated, especially when the dimensionality of the measurement vector is high or the measurement uncertainty is small. A common strategy to mitigate this issue is to expand the size and coverage of the retrieval database. Although this approach increases the likelihood of finding cases similar to the observations, the inherent measurement sparsity issue cannot be fully eliminated, as a precalculated database cannot cover all possible conditions. The quality and representativeness of the database constrain the maximum achievable performance for database-driven approaches.

Within the Bayesian framework, the Optimal Estimation Method (OEM) provides a rigorous approach to estimate state variables under prior constraints. Rodgers (2000) established a comprehensive theoretical foundation for the OEM method, and the method has since been widely applied in many atmospheric remote sensing applications (Austin et al., 2009; Pfreundschuh et al., 2020; Liu and Adams, 2025). Nevertheless, its application to submillimeter-wave radiometry for ice cloud retrievals presents several fundamental challenges. First, OEM assumes a Gaussian a priori PDF, but the actual distributions of atmospheric and cloud variables are typically strongly non-Gaussian. Second, OEM assumes a moderately nonlinear forward model, which makes selection of the initial state of iteration important. Third, as the radiative transfer model involves multiple scattering calculations in clouds, no analytical method exists for the rapid computation of the Jacobian matrix.

To leverage the strengths of the database-driven method while retaining the ability to exceed the limitations of the database when BMCI-alone fails, we adopt the hybrid BMCI and OEM strategy proposed by Evans et al. (2012) and extend the method to a state-of-the-art radiative transfer model. We present a detailed description of the retrieval algorithm and assess it using simulated observations as well as real observations collected during the IMPACTS field campaign. The paper is structured as follows. Section 2 provides a detailed description of the hybrid BMCI and OEM algorithm. Section 3 applies the hybrid algorithm to the evaluation of the C²OMODO sensor using simulated observations. Section 4 presents retrieval examples using real CoSSIR observations collected during the IMPACTS campaign and their validation against collocated triple-frequency radar measurements. Finally, Sect. 5 summarizes the results and discusses the strengths and limitations of the retrieval algorithm, as well as possible directions for future improvements.

2 Hybrid BMCI-OEM retrieval algorithm

2.1 Algorithm overview

The hybrid BMCI-OEM algorithm leverages the strengths of the BMCI method and overcomes the limitations imposed by the TB sparsity of the a priori database. With a high-quality retrieval database, BMCI is expected to serve as the primary method for producing retrievals under most conditions. When BMCI fails due to an insufficient number of database cases matching the input TB observations, OEM is applied to maximize the posterior PDF and extend the retrieval capability beyond the constraints of the a priori database. Several techniques are employed to address the OEM implementation challenges. As a result, the OEM step can minimize the TB discrepancy using prior information that is largely identical to that used in the BMCI step.

Several key steps in the OEM implementation are worth noting. First, the CDF (Cumulative Distribution Function) and EOF (Empirical Orthogonal Function) technique proposed by Evans et al. (2012) is employed to address the OEM’s requirement for a Gaussian a priori PDF. The CDFs capture the complete single-point statistical information of the a priori data, and the EOFs capture the second-order statistics characterized by the correlation matrix. Using this approach, the highly non-Gaussian a priori PDF can be transformed into a multivariate Gaussian distribution in the eigen space. Second, OEM is initialized with the BMCI estimates, which are typically close to the final optimal state. This relaxes the moderate linearity requirement of the forward model and reduces the number of iterations needed. Third, since the radiative transfer model does not have an analytical format for quick Jacobian matrix computation, we use the finite difference method by perturbing the state vector in the eigen space and computing the corresponding TB sensitivities.

2.2 Retrieval database

The retrieval database is a critical component of the hybrid BMCI-OEM algorithm. By establishing a comprehensive retrieval database, the computationally efficient BMCI method can handle most retrieval cases and produce reliable retrievals on its own. A high-quality retrieval database requires that the atmospheric and cloud samples be distributed according to prior knowledge and that the corresponding simulated TB observations densely cover the possible regimes of the measurement space.

2.2.1 A priori profiles

To support the PolSIR and PMM missions, we build a comprehensive retrieval database using CloudSat and CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation) products. Compared with cloud resolving models that rely on simplified and uncertain microphysical parameterizations (Kummerow et al., 2011), CloudSat-CALIPSO products are directly constrained by complementary active measurements and therefore provide a more observation-grounded and globally consistent representation of cloud microphysical structure across a wide range of cloud regimes. At the current stage of development, a balance between the database size and computational cost is achieved by selecting one orbit every three days from the full year of 2007 within the tropical region (30° S–30° N). Approximately 1.5 million all-sky atmospheric and cloud profiles are selected. The ice water content (IWC) profiles are obtained from the 2C-ICE product, which is retrieved using synergistic CloudSat Cloud profiling Radar (CPR) and CALIPSO lidar observations and covers a large dynamic range from thin to dense ice cloud layers (Deng et al., 2015). Since the CPR experiences significant attenuation in optically thick clouds, cases with extremely large IWP are likely underrepresented in the a priori database. The liquid water content (LWC) profiles are obtained from the 2B-CWC-RO product, a radar-only product retrieved from the CPR. The thermodynamic profiles and surface parameters are provided by the auxiliary ECMWF dataset. The vertical resolution of the atmospheric and cloud profiles is set to 240 m.

2.2.2 Radiative transfer model

We develop the radiative transfer model (RTM) using ARTS (Buehler et al., 2018), which is widely recognized as one of the most rigorous RTMs in representing multiple scattering within clouds. The developed forward model includes one ice hydrometeor type and one liquid hydrometeor type. The scattering properties of each hydrometeor type are determined by its single scattering data (SSD) and particle size distribution (PSD). For ice clouds, due to the diversity of particle shapes in nature and the significant variations in their scattering properties, we select five representative particle shapes to approximate the natural conditions. The chosen particle shapes (designated by LargePlateAggregate, LargeBlockAggregate, LargeColumnAggregate, 6BulletRosette, and EvansSnowAggregate) align with the selection in May et al. (2024). Liquid cloud particles are assumed to be spherical. All particles are assumed to be randomly oriented, and their SSDs are obtained from the ARTS scattering dataset (Eriksson et al., 2018). A one-moment PSD scheme is assumed for both ice and liquid clouds. The PSD scheme proposed by Field et al. (2007) for a tropical regime is applied to cloud ice, and the modified gamma distribution scheme is used for cloud liquid. The absorbing gases include O₂, N₂, O₃, and water vapor. The surface emissivity over ocean and land is calculated using the TESSEM2 (Tool to Estimate Sea-Surface Emissivity from Microwave to Submillimeter Waves, Prigent et al., 2017) and the TELSEM2 (Tool to Estimate Land-Surface Emissivities at Microwave frequencies, Aires et al., 2011), respectively. The DISORT (DIScrete-Ordinate-method Radiative Transfer, Stamnes et al., 1988) scattering solver is employed, and the simulations are currently limited to the first Stokes parameter. A pencil-beam is assumed in computing the TB simulations, and instrument parameters will be integrated later when the sensors are finalized.

The ice cloud particle size is computed as part of the radiative transfer calculations. Since both ice and liquid hydrometeor PSD schemes are one-moment, the IWC and LWC values from the 2C-ICE and 2B-CWC-RO products can be directly mapped to particle size parameters at each vertical layer. At each layer, the characteristic particle size D_me is defined as the mass-weighted mean of the volume equivalent diameter D_veq:

$$\begin{array}{}\text{(1)}& D_{\mathrm{me}}={\displaystyle \frac{\sum N\left(D_{\mathrm{veq}}\right)\cdot D_{\mathrm{veq}}^{\mathrm{4}}}{\sum N\left(D_{\mathrm{veq}}\right)\cdot D_{\mathrm{veq}}^{\mathrm{3}}}}\end{array}$$

where D_veq is provided by the ARTS SSD dataset, and N(D_veq) is the number concentration in each particle size interval computed by the ARTS PSD scheme. Subsequently, the column-averaged mean mass diameter D_m for each ice and liquid hydrometeor profile is computed as the IWC- or LWC-weighted mean of D_me.

2.3 BMCI algorithm

Using the precalculated retrieval database, the hybrid BMCI-OEM algorithm first performs the BMCI retrieval. Since the database samples are already distributed according to the a priori PDF, the posterior PDF is determined by the likelihood function, which quantifies how well the database TB simulations match the observations relative to the measurement uncertainties:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& P_{\mathrm{post}}\left(x\right|y_{\mathrm{obs}})\propto \exp \left(-\mathrm{1}/\mathrm{2}{\mathit{\chi }}_y^{\mathrm{2}}\right)\\ & {\mathit{\chi }}_y^{\mathrm{2}}=\sum _{j=\mathrm{1}}^N{\displaystyle \frac{(y_{\mathrm{sim},j}-y_{\mathrm{obs},j}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_j^{\mathrm{2}}}}\end{array}\end{array}$$

where y_sim,j and y_obs,j are the simulated and observed TB in the jth channel, σ_j is the measurement noise characterized by the noise equivalent differential temperature (NEDT), and N is the number of radiometer channels. The retrieval result and associated uncertainty are derived by integrating over the posterior PDF to find the mean state and standard deviation:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& x_{\mathrm{ret}}=\sum _i{x}_i\cdot p_{\mathrm{post}}\left(x_i\right|y_{\mathrm{obs}})\\ & {\mathit{\sigma }}_x=\sqrt{\sum _i(x_i-x_{\mathrm{ret}}{)}^{\mathrm{2}}\cdot p_{\mathrm{post}}(x_i\left|y_{\mathrm{obs}}\right)}\end{array}\end{array}$$

The BMCI algorithm does not require the prior PDF to follow a specific functional form. It is also computationally efficient as no additional radiative transfer calculations are needed. However, its primary limitation is the database sparsity in the measurement space. The successful implementation of the BMCI requires that a minimum number of database cases (25 in this study, consistent with Evans et al., 2005) fall within a specified ${\mathit{\chi }}_y^{\mathrm{2}}$ threshold. Due to the finite size of the prior database, this requirement cannot always be satisfied, particularly when the number of sensor channels is large or the measurement noise for each channel is small. In such cases, the BMCI algorithm theoretically fails. One method to obtain reliable retrievals and uncertainty estimates is to inflate the measurement uncertainty NEDT until the specified number of database cases is reached, after which the posterior PDF is integrated, as shown in Eq. (3). This approach weakens the observational constraint imposed on the retrieval and is therefore not optimal. Alternatively, as developed in this study, the OEM optimization procedure is implemented to further maximize the posterior PDF.

2.4 OEM Algorithm

One major challenge in implementing the OEM algorithm is that the OEM requires the a priori PDF to be Gaussian, but the distribution of cloud variables in the retrieval database is highly non-Gaussian. To address the mismatch, we apply the CDF/EOF technique developed by Evans et al. (2012) to capture the statistical characteristics of the retrieval database and transform the a priori information into the eigen space. The core idea of the CDF method is to map the rank of a state variable at a given layer in the retrieval database to the CDF of the standard normal distribution. Through the rank-based mapping, the atmospheric variables and their Gaussian representations can be mutually transformed. For the retrieval database, the CDFs of the state variables including IWC, LWC, and water vapor at each layer are computed by sorting the parameter values and assigning ranks. To avoid discontinuities in ranking, layers identified as clear are replaced with small random values below the instrument sensitivity threshold. The hydrometeor parameters are represented by their ranking. At each layer, the inverse standard normal CDF is applied to the ranked values to transform them into points following a Gaussian distribution. The CDF mapping process is applied layer by layer to the profiles of IWC, LWC, and water vapor, and the database samples are transformed into Gaussian space. The covariance matrix S_ga between different layers and variables of the Gaussian-distributed points is computed.

The EOF method is further applied to capture the correlations between variables at different levels and reduce the state vector dimensionality and computational cost. The eigenvalues and eigenvectors of the covariance matrix mathbfS_ga are computed. The eigenvectors define a new set of orthogonal coordinates, with the variability of the data along each coordinate characterized by the corresponding eigenvalue. By defining the OEM state vector in eigen space, the cost function to be minimized in the OEM is given by:

$$\begin{array}{}\text{(4)}& J=\sum _{i=\mathrm{1}}^{N_{\text{dim}}}{\mathit{\xi }}_i^{\mathrm{2}}+\sum _{j=\mathrm{1}}^{N_{\text{chan}}}{\displaystyle \frac{{\left(y_{\text{obs},j}-F_j\left[\mathit{G}\right(\mathit{\xi }\left)\right]\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_j^{\mathrm{2}}}}\end{array}$$

where ξ is the state vector, G(x) denotes the function mapping the state vector from the eigen space back to the corresponding atmospheric and cloud variables, and F(x) represents the forward model simulations. The eigenspace state vector ξ consists of the corresponding eigenspace representations of the original IWC, LWC, and water vapor profiles. ξ is first transformed back to Gaussian space ξ_g by scaling each EOF component by the square root of the eigenvalue and then multiplying by the corresponding eigenvectors to rotate back to the original coordinate system. The vector ξ_g is then converted back to atmospheric and cloud profiles. For each ξ_g element, its standard normal CDF value is computed and mapped back to the corresponding geophysical parameters by interpolating the relationship between CDF values and cloud parameters at the same layer in the a priori database. This procedure is applied sequentially to all elements of the Gaussian state vector ξ_g.

The OEM iteration starts with the BMCI-estimated results. When an insufficient number of database cases is found, we inflate the measurement uncertainty NEDT in a step of $\sqrt{\mathrm{2}}$ until reaching the required number of database cases (Evans et al., 2005; Liu et al., 2018), after which the mean state estimate is computed. The state vector in the eigen space corresponding to the BMCI estimate is computed and used as the initial guess for the OEM iteration. As no analytical method is currently available for computing Jacobians that include multiple-scattering effects in ARTS, the Jacobian matrix is computed numerically by perturbing each state element in eigen space and computing the corresponding TB sensitivities. With the Jacobian matrix, the robust Levenberg-Marquardt (LM) method is employed to iteratively update the state vector. The iteration stops when reaching a stable state, and the associated uncertainty estimate is computed as:

$$\begin{array}{}\text{(5)}& S_{\mathrm{p}}={\left({\mathbf{K}}^T{\mathbf{S}}_y^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\end{array}$$

where S_a is the a priori covariance matrix of the state vector ξ, S_y is the covariance matrix of measurement error, and K is the Jacobian matrix of F(G(ξ)). Other retrieval diagnostics such as the averaging kernel matrix are produced as well.

Since Eq. (5) provides uncertainty estimates of the state vector in eigen space rather than for the geophysical variables of interest, we apply a Cholesky decomposition of the posterior covariance matrix to generate an ensemble of random cases around the converged state in eigen space (Evans et al., 2012). The ensemble members are then transformed back to atmospheric and cloud variables, and MCI is performed to compute the uncertainty estimates for the variables of interest, including vertical profiles of ice and liquid cloud microphysics and water vapor, as well as the corresponding column-integrated quantities such as IWP, D_m, liquid water path (LWP), and integrated water vapor (IWV).

3 Simulation experiment

3.1 Retrieval accuracy

In this section, we apply the hybrid BMCI-OEM algorithm to assess the performance of the C²OMODO sensor in the PMM missions using simulated observations. C²OMODO is a satellite mission supported by the French space agency Centre National D'Études Spatiales as a contribution to the JAXA-NASA PMM mission. C²OMODO includes one window channel at 89 GHz, four channels at the 183.31 GHz absorption line, and five channels at the high-frequency 325.15 GHz absorption line. The detailed channel frequencies and corresponding measurement uncertainties NEDT are summarized in Table 1. It should be noted that the listed NEDT characterizations are used only to evaluate the algorithm performance. The true values remain undetermined. The targeted footprint size at nadir for the 89, 183, and 325 GHz channels are 10, 5, and 3 km, respectively.

Table 1 Channel characteristics assumption of the C²OMODO radiometer. The NEDT values correspond to those used in the simulation experiment. The actual values remain undetermined.

In the simulation experiment, independent cloud profiles for testing are obtained in the same way as the a priori database but using CloudSat products from a different year. From the full year of 2008 CloudSat data, atmospheric and cloud profiles from one orbit every three days are first extracted, following the same procedure as for the a priori database. From these, 30 000 all-sky atmospheric and cloud profiles are randomly selected, approximately half of which contain at least one ice cloud layer. The same radiative transfer model used to construct the retrieval database is employed to compute the simulated TB observations. After adding Gaussian noise with a standard deviation set according to the assumed NEDT, the simulated observations are input into the hybrid retrieval algorithm, and the retrieved results are compared with the reference values for evaluation. This simulation experiment uses an idealized configuration in which the identical forward model with the same particle habit selection is applied for both simulating the observations and performing the retrievals. The figures presented in this section are based on the LargeBlockAggregate particle habit assumption by default unless noted otherwise. We have also evaluated the impact of inconsistent microphysics assumptions between the a priori database and the retrievals. However, this is not the focus of the current study and therefore the results are not included. All results presented here assume nadir-only observations: we acknowledge that both the PolSIR and PMM-C²OMODO radiometers are cross-track scanners, and results at off-nadir angles can be different.

Figure 1 presents the retrieval results for the C²OMODO radiometer in the idealized simulation experiment based on the hybrid BMCI-OEM algorithm. The top panels show the retrieval results for IWP, LWP, IWV, and D_m compared with the reference, and the bottom panels display the associated retrieval uncertainties. The C²OMODO sensor channel set demonstrates good retrieval capability for IWP. Since no inconsistencies in particle habit or PSD selections in the forward model are included in this idealized experiment setting, the median retrievals follow the diagonal line until IWP falls below approximately 30 g m⁻², where the retrieved values start to exhibit systematic negative biases for thin ice clouds. The black dotted lines indicate the −50 % and +100 % accuracy bounds that are defined as the range between 0.5 ⋅ IWP_true and 2 ⋅ IWP_true for reference. At high IWP values, the scatter points closely follow the diagonal with minimal spread, and the dispersion increases as IWP decreases. For retrieval uncertainties, when IWP exceeds approximately 100 g m⁻², most uncertainty estimates fall within the 100 % range. For smaller IWP values, the uncertainty estimates increase rapidly. Since C²OMODO has an 89 GHz channel that is sensitive to both ice scattering and liquid emission, the sensor demonstrates some capability for retrieving liquid cloud total mass when LWP is beyond 100 g m⁻². The retrieved LWP values generally exhibit a small negative bias, and the associated retrieval uncertainties remain high. This is expected due to the absence of low-frequency channels that are primarily sensitive to liquid hydrometeor emissions. The C²OMODO sensor also demonstrates the capability of constraining water vapor as both 183 and 325 GHz channels are centered on water vapor absorption lines. The retrieved IWV values closely follow the true values with little scatter throughout the entire dynamic range, and the uncertainty estimates are generally below 5 kg m⁻². For the ice cloud particle size, when the true D_m exceeds roughly 400 um, the median retrieval closely follows the diagonal. For D_m below 400 um, the retrieved values tend to underestimate the truth. The corresponding uncertainty estimates remain mostly below 100 um for D_m larger than 400 um, and they increase substantially for smaller ice particles.

Figure 1 Retrieval results for the C²OMODO radiometer in an idealized simulation experiment using the hybrid Bayesian Monte Carlo Integration (BMCI) and Optimal Estimation Method (OEM) retrieval algorithm. The top panels compare the retrieved ice water path (IWP), liquid water path (LWP), integrated water vapor (IWV), and mean mass diameter (D_m) with the true values, and the bottom panels show the corresponding retrieval uncertainties.

Figure 2 compares the BMCI-only and hybrid BMCI-OEM retrievals to investigate the benefits of the additional OEM procedure. When the BMCI fails to find enough database matches, the BMCI-only results are obtained by inflating the NEDT in a step of $\sqrt{\mathrm{2}}$ until the required number of cases is reached. Approximately 7 % of the all-sky testing samples require the NEDT inflation step, and these are the same set of samples for which the OEM procedure is applied. The left panel compares the ${\mathit{\chi }}_y^{\mathrm{2}}$ values for the BMCI-only and the OEM retrievals. As the OEM iteration is initialized with the BMCI estimate, the comparison also illustrates how the ${\mathit{\chi }}_y^{\mathrm{2}}$ values change before and after applying the OEM step. The figure shows a substantial reduction in ${\mathit{\chi }}_y^{\mathrm{2}}$ after applying the OEM. For the BMCI-only results, many ${\mathit{\chi }}_y^{\mathrm{2}}$ values exceed 10², and some exceed 10³. After the OEM step, the ${\mathit{\chi }}_y^{\mathrm{2}}$ decreases significantly, with most values falling below 10¹ and within the range expected from the measurement uncertainty. One more point to note is that the BMCI-only failures occur over a wide range of IWP values. The middle and right panels compare the retrieval results in terms of IWP and D_m in this idealized OSSE. The results indicate that the OEM step provides moderate but consistent improvements in accuracy relative to the BMCI-only retrievals. The improvements are particularly evident for IWP values greater than 1 kg m⁻² and for large particle sizes. It should be noted that the BMCI retrieves the column-integrated IWP and D_m directly, whereas the OEM first retrieves the vertical profiles and then performs the vertical integration. Since these two retrieval results are not necessarily equivalent, the improvements achieved through the OEM step become particularly noteworthy.

Figure 2 Comparison of the BMCI-only and hybrid BMCI-OEM retrievals in an idealized simulation experiment, with the left panel showing the ${\mathit{\chi }}_y^{\mathrm{2}}$ values that quantify the mismatch between the simulated and observed brightness temperature, the middle panel showing the retrieval accuracy for IWP, and the right panel showing the retrieval accuracy for the mean mass diameter D_m.

3.2 Retrieval diagnostics

The OEM framework also provides diagnostic parameters for assessing the quality and reliability of the retrieval results. In this study, we focus on three key diagnostics: vertical resolution, degrees of freedom (DoF), and Shannon information content. To make a more comprehensive analysis, even in cases where BMCI alone can perform the retrieval, a single OEM iteration is still executed to compute the Jacobian matrix and associated diagnostic quantities. For the cases in which OEM is applied, the final Jacobian matrix and associated diagnostic quantities are used in this evaluation.

To assess the vertical resolution, the averaging kernel matrix is first calculated:

$$\begin{array}{}\text{(6)}& A={\left({\mathbf{K}}^T{\mathbf{S}}_y^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}{\mathbf{K}}^T{\mathbf{S}}_y^{-\mathrm{1}}\mathbf{K}\end{array}$$

The vertical resolution at each level is then computed as the full width at half maximum (FWHM) of the corresponding row of the averaging kernel matrix. By separating the elements of the Jacobian matrix into ice cloud, liquid cloud, and water vapor, the averaging kernel matrix and the vertical resolution are computed for different types separately. Figure 3 shows the median and interquartile-range (IQR) profiles of the vertical resolution for ice, liquid, and water vapor. The median vertical resolution of ice clouds is seen to range between 2 and 3 km, with the IQR spanning 1 to 4 km. The resolution is typically smaller at the cloud top than at the bottom, with a median of about 2 km near the top and around 2.7 km below 11 km. The results are highly consistent with those in May et al. (2025). The C²OMODO channel set also shows sensitivity to liquid clouds, with a median vertical resolution also in the 2–3 km range. For water vapor, the median resolution is approximately 4 km under all-sky atmospheric conditions.

Figure 3 The median and interquartile-range (IQR) profiles of the vertical resolution for the ice hydrometeor, liquid hydrometeor, and water vapor, as derived from the averaging kernel matrix in an idealized simulation experiment.

From the averaging kernel matrix, the DoF is computed to quantify the effective number of independent pieces of information provided by the observations. The DoF is computed as the trace of the averaging kernel matrix in Eq. (6). Figure 4 shows the scatter plot of the DoF as a function of IWP and IWV. The DoF generally increases with both IWP and IWV. When the ice cloud is absent or very thin (IWP $<{\mathrm{10}}^{-\mathrm{2}}$ kg m⁻²), DoF < 4 typically occurs in conditions that are dry (IWV < 40 kg m⁻²). DoF values between 4 and 5 are more common when IWV ranges between 40 and 55 kg m⁻², while DoF > 5 are frequent when IWV exceeds 55 kg m⁻². These low IWP values correspond to the “clear-sky” retrieval capability for the water vapor profile. Specifically, using the assumed C²OMODO channel set, 4–6 DoFs for water vapor can be independently retrieved, consistent with the overall SAPHIR (Sondeur Atmosphérique du Profil d’Humidité Intertropicale par Radiométrie) water vapor retrieval capability on which the C²OMODO instrument is based (Mathur et al., 2013). The DoF is mostly 5–6 when a relatively thin ice cloud is present (10⁻² kg m⁻² < IWP < 10⁻¹ kg m⁻²) and increases to 6–7 for thicker ice clouds where the sensor frequencies are particularly sensitive (IWP between $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{1}}$ and 3 kg m⁻²). The larger DoF observed when an ice cloud is substantial corresponds to the sensor’s ability to retrieve information about the vertical distribution of ice cloud mass and particle size, while the capability to resolve the vertical distribution of water vapor is lost, as indicated by the flattening of the gradient on the ordinate. The DoF decreases to about 4–5 at the highest IWP values, indicating information saturation when the ice cloud becomes too thick.

Figure 4 The scatter plot of the Degrees of Freedom (DoF) as a function of ice water path (IWP) and integrated water vapor (IWV) that are derived from an idealized simulation experiment.

The Shannon information content measures the reduction of uncertainty in the retrieved state compared to the prior PDF and quantifies how much information the observations contribute beyond the prior knowledge. The Shannon information content is computed as:

$$\begin{array}{}\text{(7)}& H=\mathrm{1}/\mathrm{2}{\log }_{\mathrm{2}}\left|{\mathbf{S}}_{\mathrm{a}}{\mathbf{S}}_{\mathrm{p}}^{-\mathrm{1}}\right|\end{array}$$

where S_a and S_p denote the covariance matrices of the a priori and posterior PDFs of the state vector, respectively. The value H is expressed in bits, where n bits indicate that the measurements reduce the state uncertainty by a factor of 2ⁿ. As the information from different channels is partially redundant, the sum of per-channel information is not equal to the total information content when all channels are combined. To quantify the contribution of each channel, we compute the incremental information content by starting with a single channel and adding channels one at a time until all channels are included. The sequence begins with the 183.31–10.7 GHz channel, followed by the remaining 183 GHz channels, then the 325 GHz channels, and finally the 89 GHz channel. Figure 5 presents the incremental Shannon information content for IWP, LWP, and IWV, respectively. For IWP, the information content generally increases with increasing IWP but decreases at the high end. Starting from around 20 g m⁻², the information content rises at an approximately constant rate up to about 1 kg m⁻², beyond which the rate of growth slows. The combined 183 GHz channels provide approximately 5 bits of information at an IWP of 100 g m⁻² and reach a maximum of about 10 bits at an IWP of approximately 3 kg m⁻². Adding a single 325 GHz wing band (325 ± 10.7 GHz) provides approximately 1 bit of information even when IWP is extremely small, and the five 325 GHz channels contribute roughly 4 bits of information for higher IWP. The addition of channel 9 (325 ± 0.8 GHz) produces the largest incremental increase in information content among the 325 GHz channels. This result suggests that a channel closer to the water vapor absorption line at 325 GHz may better sense cloud-top information or more effectively distinguish between cloud and water vapor signals near the cloud top. The 89 GHz channel contributes noticeably only at the high end of the IWP range, probably because of its unique role in separating liquid from ice hydrometeors as precipitation often occurs in this IWP range.

For LWP, channels 1 to 3 show some sensitivity to liquid cloud, but channels 4 to 9 contribute very little information due to strong attenuation by water vapor and overlying ice clouds. The 89 GHz channel produces the largest increase in information content for liquid clouds. For IWV, the 183 GHz channels exhibit clear sensitivity, with each additional subchannel producing a noticeable increase in information content across the entire IWV range. Among the 325 GHz channels, only channel 9, which is closest to the water vapor absorption line, exhibits evident sensitivity, while the other channels contribute very little. The 89 GHz channel is situated in an atmospheric window but is influenced by the water-vapor continuum, and it is observed to provide substantial information for IWV. This Shannon information content analysis demonstrates the scientific merit of the new 325 GHz submillimeter bands and also highlights the importance of including traditional high-frequency channels (89 and 183 GHz) to achieve high-quality all-sky retrievals of ice and liquid hydrometeors as well as water vapor.

Figure 5 The Shannon information content of C²OMODO radiometer for ice water path (IWP), liquid water path (LWP), and integrated water vapor (IWV). The incremental information content is computed by starting with a single channel and adding channels one at a time until all channels are included. The sequence begins with the 183.31–10.7 GHz channel, followed by the remaining 183 GHz channels, then the 325 GHz channels, and finally the 89 GHz channel.

4 Field campaign retrieval validation

In this section, we apply the hybrid BMCI-OEM algorithm to real observations collected during the IMPACTS field campaign to further evaluate the retrieval performance. IMPACTS is a NASA-sponsored field campaign to study wintertime snowstorms associated with East Coast cyclones. The campaign flew a complementary suite of remote sensing and in-situ instruments during the boreal winters of 2020–2023 to collect measurements essential for snow process studies. CoSSIR has been upgraded with new front-end receivers, configurable scanning patterns, and dual-polarization capability across all channels, thereby increasing its flexibility for simulating a variety of spaceborne submillimeter radiometers, including PolSIR and PMM-C²OMODO. The newly upgraded CoSSIR was first deployed aboard the ER-2 aircraft during the final phase of the IMPACTS campaign in the 2022–2023 winter and operated in a hybrid conical and along-track scanning mode. This was also the first time CoSSIR flew over mid- and high-latitude winter weather conditions. The CoSSIR instrument includes four subchannels near 183 GHz, three subchannels near 325 GHz, and one channel at 684 GHz, all of which are dual-polarized. The detailed characteristics of channel frequency and NEDT are listed in Table 2. CoSSIR has an approximately 4.3° beam width, resulting in a surface footprint of about 1.5 km at nadir and approximately 2.5×4.2 km at 53° off-nadir when the ER-2 aircraft cruises at an altitude of 20 km. On the same ER2 aircraft platform, the W-band CRS (Cloud Radar System, McLinden et al., 2021) and the dual-frequency Ka/Ku-band HIWRAP (High-altitude Imaging Wind & Rain Airborne Profiler, Li et al., 2015) radars were also deployed. The sensitivities of the W-, Ka-, and Ku-band radar channels are −30, 1.5, and 7.8 dBZ, respectively. The collocated triple-frequency radar observations provide detailed cloud structure and a reliable reference for validating the hybrid BMCI-OEM algorithm’s performance from CoSSIR observations.

Table 2 Channel characteristics of the CoSSIR radiometer during the IMPACTS field campaign.

4.1 CoSSIR retrieval settings

A major difference from the retrieval implementation in Sect. 3 is the update of the a priori database. Since the retrieval database prepared for PolSIR and PMM-C²OMODO in Sect. 2.2 is specific to the tropical region, we develop a separate a priori database for the IMPACTS field campaign. The construction of the new database follows the same procedure described in Sect. 2.2, but it uses different subsets of CloudSat products. All CloudSat datasets from winter 2006–2019 are examined, and profiles with IWP greater than 10 g m⁻² are selected when they fall within a selected latitude-longitude box over the East Coast of the United States that covers the IMPACTS region. The same radiative transfer model and same selection of the ice cloud particle habits described in Sect. 2.2.2 are applied. With this new a priori database, the hybrid BMCI and OEM algorithm is applied to the real CoSSIR observations to produce retrieval results. Only the nadir TB observations are used to enable a direct comparison with retrievals from the collocated triple-frequency radar observations, with the latter serving as the reference. The atmospheric profiles and surface parameters required for the radiative transfer calculations are obtained from the Weather Research and Forecasting (WRF) model outputs supporting the IMPACTS field campaign.

4.2 Triple-frequency radar retrieval settings

One advantage of the hybrid BMCI-OEM algorithm is its great extendibility to other types of remote sensing observations. By updating the forward model to include a radar simulator and replacing the TB observations with radar reflectivity profiles in the measurement vector, the hybrid BMCI-OEM algorithm is directly adapted for triple-frequency radar retrievals. The same ARTS forward model used for TB simulations is upgraded to include the ARTS built-in single-scattering radar solver for simulating radar reflectivities at different bands.

The triple-frequency radar observation profiles are interpolated to the same 240 m vertical resolution as the a priori database. Since only ice clouds are of interest in this study, we only focus on layers where at least one of the triple-frequency radar reflectivities exceeds its minimum detection limitation and the temperature, according to the WRF temperature profile, is below 273 K. The range bins below the melting layer are excluded. The state vector is defined as the IWC in these layers, and the measurement vector consists of the triple-frequency radar observations in these layers that are above their respective detection thresholds. To obtain the initial profile for the OEM iteration, we build a pre-calculated database containing random IWC samples and their corresponding triple-frequency radar reflectivity simulations. In each layer, the IWC is determined by selecting the database sample whose simulated reflectivities best match the triple-frequency observations. The initial profile is then generated using an onion-peeling approach that performs the database-matching process from the top downward, with the retrieved attenuation at each step used to correct the measured reflectivity of the next layer (Liu and Mace, 2022; Liu et al., 2022). During the radar OEM retrieval, the same $\mathrm{CDF}/\mathrm{EOF}$ method and identical a priori constraints as those used in the CoSSIR retrievals in Sect. 4.1 are applied. The IWC and D_me profiles retrieved from the triple-frequency radar are then integrated to obtain the column-integrated IWP and D_m for the validation of the CoSSIR retrieval results.

4.3 Retrieval example

In this subsection, we select one representative cloud scene and present detailed retrieval results from CoSSIR observations using the BMCI-only and hybrid BMCI-OEM approaches, along with retrievals from the triple-frequency radar observations for a comprehensive comparison. The selected cloud scene occurs between 16:52 and 17:02 UTC on 15 January 2023. This flight leg takes place during the tropical transition phase of an extra tropical storm, resulting in a particularly complex hydrometeor vertical distribution. Figures 6–9 show the results obtained using the LargeBlockAggregate ice particle habit for both the radiometer and radar retrievals.

Figure 6 shows the CoSSIR TB observations in blue for different channels, while the orange and green lines show the TB simulations using retrievals from the BMCI and OEM, respectively. Again, the BMCI-only results are obtained by incrementally inflating the measurement uncertainty until the required minimum of 25 database cases are identified. The OEM iteration is initialized using the BMCI estimates, and these two lines evaluate the effectiveness of the OEM in reducing TB discrepancies. It is evident that the TB agreement improves significantly after applying the OEM step. Using the BMCI-only results, the simulations reproduce the observations of some channels reasonably well, but substantial negative biases remain, especially when the TBs have deep depressions around 16:57 UTC when the aircraft flew above thick clouds. Among the 183 GHz channels, the 182.31 GHz channel shows the largest discrepancy, with agreement improving as the channel moves away from the water vapor absorption line. A similar pattern is seen for the 325 GHz channels, where the ± 0.9 GHz channel exhibits the largest mismatch, and the ± 11.5 GHz channel shows the best agreement. This trend is likely due to the greater sensitivity of channels closer to the center of the water vapor absorption line to ice clouds near the cloud top. The BMCI retrievals, however, exhibit numerous artificial ice clouds near the cloud top that should not be present, as will be discussed in Fig. 8. After applying the OEM, all channels except 325 ± 0.9 GHz show substantial improvement, with simulations closely reproducing the TB observations. Considering that the radiative transfer model has many simplifications and assumptions, achieving such good agreement for both clear-sky and thick-cloud scenes is impressive. The results also suggest that the BMCI-only approach may overestimate the overall cloud ice mass due to the consistent cold biases in the simulated TBs relative to the observations.

Figure 6 The blue lines show the CoSSIR brightness temperature (TB) observations for a representative cloud scene occurred on 15 January 2023 during the IMPACTS field campaign. The orange and green lines show the TB simulations computed from the BMCI-only and OEM retrieval results.

The left-column panels in Fig. 7 show the collocated W-, Ka-, and Ku-band radar reflectivity profiles, respectively, while the right-column panels show the radar reflectivity residuals defined as reflectivity simulations minus observations for the W-, Ka, and Ku-band retrievals. Overall, the simulated radar reflectivity profiles reproduce the observations reasonably well, with most differences falling within ± 4 dB. The W-band results show some negative biases within the frontal band. These biases are expected partly due to the saturation of the W-band in very thick ice clouds, and partly as a compensation effect to minimize the cost function. For Ku band, biases are positive near echo top and become negative and largest above the melting band. As explained above, we do not perform LWC retrievals below the 273 K isotherm, nor do we aim to validate the algorithm for LWP retrievals, primarily due to the limitations of ARTS and other existing RTMs in accurately modeling radiative transfer within the melting layer.

Figure 7 The left-column panels show the collocated W-, Ka-, and Ku-band radar reflectivity profiles corresponding to the same scene as the CoSSIR observations in Fig. 6. The right-column panels present the radar reflectivity residuals that are defined as the simulated reflectivity minus the observations after applying the triple-frequency radar retrieval algorithm.

Figure 8 shows the IWC and IWP retrievals from CoSSIR observations using the BMCI-only algorithm (left column) and hybrid BMCI-OEM algorithm (middle column), compared with the reference retrievals from triple-frequency radar observations (right column). Each column shows the retrieved IWC profiles at 240 m vertical resolution (top), the corresponding IWC uncertainties (middle), and the IWP and associated uncertainty estimate (bottom). Compared to the radar retrievals that provide more detailed and smoothly varying cloud profiles, the radiometer retrievals generally capture the overall cloud structure, but finer internal details cannot be fully reproduced due to the limited information content in the TB observations. As the DoF shown in Fig. 4 mostly ranges from 4 to 7, the radiometer observations have limited capability to constrain the vertical profiles, and a large portion of the constraints comes from the a priori information. For the BMCI-only retrievals, the most notable feature is the vertical extension beyond the expected cloud layers at both the cloud top and bottom. The internal cloud structure is also poorly characterized. Both sides of this flight leg (16:52–16:54 and 17:00–17:01 UTC) contain almost no ice mass, but the BMCI-only algorithm indicates the presence of ice clouds. The OEM optimization shows the ability to moderately remove spurious cloud layers and to better highlight the internal structure. Although the OEM retrievals still show clouds on both sides of the curtain, these clouds are assigned to the correct altitudes with much thinner geometric thickness, making them more physically consistent with the radar reference for the anvil clouds. The uncertainty estimates from both algorithms show a similar pattern, with the cloud cores becoming more tightly constrained and the cloud top and bottom layers containing larger errors. Again, the hybrid BMCI-OEM method correctly retrieves the anvil cloud layer height on both sides, while the BMCI-only method does not. For the IWP retrievals shown in the bottom panels, the major IWP changes after applying OEM occur around 16:57 UTC, coinciding with the largest TB depressions. While the BMCI-only retrievals show better agreement with the radar results during this period than the BMCI-OEM method, a more thorough comparison using a larger set of retrieval cases will be presented in Sect. 4.4. The radiometer-retrieved IWP and the radar retrievals show a similar trend, but notable differences are observed for the thick ice clouds in this selected retrieval example. For the radar retrievals, the uncertainty estimates for both IWC and IWP are small, reflecting the rich information content in the triple-frequency radar observations in constraining IWC.

Figure 8 Ice water content (IWC) and ice water path (IWP) retrievals from CoSSIR observations using the BMCI-only algorithm (left column) and BMCI-OEM algorithm (middle column), with retrievals from the triple-frequency radar observations shown in the right column as the reference. Within each column, the retrieved IWC profiles, the corresponding IWC uncertainties in percentage, and the IWP along with its uncertainty estimate are presented, respectively.

Figure 9 shows the retrieval results for ice cloud particle size, and they are presented in the same format as Fig. 8. As the calculation of D_me depends highly on the IWC due to the one-moment PSD assumption, the resulting particle size retrievals are similar to the IWC results in Fig. 8. The radiometer retrieval results generally follow the radar-derived results, but they do not reproduce the detailed features in the D_me profiles. For the mean mass diameter D_m shown in the bottom panels, the retrieved values reach up to approximately 750 um and are capped at this value. The limitation comes from both the maximum particle size in the SSD dataset and the specific PSD selected for the retrievals. Compared with the BMCI-only results, the OEM retrievals partially improve the positive biases on both sides of this cloud curtain. These large ice particles with very small IWC values are highly unlikely to exist in nature. In the future, we plan to develop customized criteria to filter out such spurious clouds or implement a cloud-flagging algorithm upfront to prevent such occurrences. One thing to note is that the BMCI-only retrievals are computed by directly applying BMCI to the D_m in the retrieval database, while the OEM retrievals are obtained by first computing the particle size D_me profiles and then performing the IWC-weighted calculations. For the radar retrievals, the detailed inner structure is well resolved, and the uncertainty estimates are small for both D_me profiles and the column-averaged D_m parameter.

Figure 9 Same as Fig. 8, but showing the vertical profiles of particle size D_me and the corresponding column-averaged mean mass diameter D_m.

4.4 Retrieval Statistics

In this subsection, we include additional CoSSIR observations and a broader set of ice cloud particle habits to provide a more comprehensive assessment of the radiometer retrievals relative to the collocated triple-frequency radar retrievals. Four deep cloud scenes on 15 January 2023 are selected, and five particle habits described in Sect. 2.2.2 are applied.

The top row of panels of Fig. 10 shows the W-band CRS observations for the four selected cloud scenes observed on 15 January 2023. These observations were collected between roughly 16:00 and 18:30 UTC. The second row of panels shows the ${\mathit{\chi }}_y^{\mathrm{2}}$ values for the BMCI-only retrievals, and the third row shows the ${\mathit{\chi }}_y^{\mathrm{2}}$ values after applying the OEM. Results for all five different ice cloud habits are presented. The results further highlight the merits of OEM in reducing the TB discrepancies. For the BMCI-only retrievals, most ${\mathit{\chi }}_y^{\mathrm{2}}$ values fall within the range of 10²–10³, except in the higher reflectivity cores of the clouds where the ${\mathit{\chi }}_y^{\mathrm{2}}$ values often exceed 10³. The results differ by particle habits, with the LargeColumnAggregate and EvansSnowAggregate habits in general showing higher ${\mathit{\chi }}_y^{\mathrm{2}}$ values than the other three habit assumptions. The improvements after applying the OEM are substantial. In general, the ${\mathit{\chi }}_y^{\mathrm{2}}$ values decrease by one to two orders of magnitude and fall below or close to the expected values due to the measurement uncertainties, given that ${\mathit{\chi }}_y^{\mathrm{2}}$ follows a chi-squared distribution with a mean of M and a standard deviation of $\sqrt{\mathrm{2}M}$, where M is the number of channels. Even in the dense ice cloud regions where the BMCI-only ${\mathit{\chi }}_y^{\mathrm{2}}$ values are extremely large, the OEM step reduces them to very small values. For the LargeColumnAggregate and EvansSnowAggregate, the ${\mathit{\chi }}_y^{\mathrm{2}}$ values do not reduce to comparable low levels, indicating the critical role of particle habit selection in determining retrieval performance.

Figure 10 The top row shows the W-band CRS observations for four selected cloud scenes observed on 15 January 2023. The second row of panels shows the ${\mathit{\chi }}_y^{\mathrm{2}}$ values for the BMCI-only retrievals, and the third row presents the ${\mathit{\chi }}_y^{\mathrm{2}}$ values after applying the OEM. Results for five different ice cloud habits are presented.

Figure 11 compares the overall performance of CoSSIR and triple-frequency radar retrievals, with the first two columns showing layer-level parameters (IWC and D_me) and the third and fourth columns showing column-integrated parameters (IWP and mean mass diameter D_m), respectively. The BMCI-only and BMCI-OEM retrievals derived from CoSSIR are both evaluated against the corresponding triple-frequency radar retrievals. The fifth column shows the median absolute logarithmic error (MedALE) for each variable as a quantitative metric of retrieval accuracy. MedALE is defined as:

$$\begin{array}{}\text{(8)}& \mathrm{MedALE}=\mathrm{median}\left(\left|\mathrm{10}{\log }_{\mathrm{10}}\left({\displaystyle \frac{x_{\mathrm{cossir}}}{x_{\mathrm{radar}}}}\right)\right|\right)\end{array}$$

where x_cossir and x_radar denote the retrievals from the CoSSIR and triple-frequency radar observations, respectively. MedALE is expressed in units of dB. Each row presents the results obtained using a different particle habit selection. It is seen that the retrieval performance varies with particle habit selection. Among the five selected particle habits, the LargePlateAggregate, LargeBlockAggregate and 6BulletRosette exhibit better constrained retrievals, while the LargeColumnAggregate and EvansSnowAggregate show poorer performance. Compared with the BMCI-only approach, the BMCI-OEM retrievals provide moderate improvements for all four ice cloud variables. For the IWC retrievals shown in the first column, the OEM step moderately reduces the positive biases present in the BMCI-only results and reduces the spread of the scatter. The IWC values retrieved from radiometer observations generally follow the radar results, particularly for clouds with large IWC. This indicates that the TB observations are sensitive to layer-resolved ice cloud microphysics. However, as the information content in the TB observations is limited, the radiometer retrievals still show a spread of nearly one order of magnitude relative to the radar results. The MedALE for IWC ranges from 3.36 dB for the 6BulletRosette to 6.12 dB for the LargeColumnAggregate habit when applying the hybrid BMCI-OEM algorithm. For the D_me profile retrievals, the BMCI results are widely scattered. The OEM method reduces the divergence, but the retrieved particle sizes show persistent overestimation for D_me<700 um relative to the triple-frequency radar retrievals. For the IWP retrievals shown in the third column, the OEM tends to consistently improve the biases, and the scattered samples are generally better aligned with the diagonal. The MedALE for IWP shows evident decrease after implementing OEM. For the LargePlateAggregate, LargeBlockAggregate, and 6BulletRosette habits, MedALE decreases from the 3–4 dB range to the 2–3 dB range. Improvements in D_m retrievals are also apparent, as the OEM partially mitigates the BMCI-only bias of overestimating low values and underestimating high values.

Finally, since the LargePlateAggregate, LargeBlockAggregate, and 6BulletRosette habits reproduce the TB observations well (Fig. 10), we average their MedALE values to quantify the retrieval accuracy. After applying the hybrid BMCI-OEM algorithm, the differences between the CoSSIR-retrieved IWC, D_me, IWP, and D_m and the triple-frequency radar retrievals are 3.75, 0.82, 2.53, and 0.55 dB, respectively. The improvements of the hybrid BMCI-OEM algorithm over the BMCI-only approach are 0.21, 0.02, 0.85, and 0.13 dB for these four variables.

Figure 11 The scatterplots of the layer-level parameters including ice water content (IWC) and particle size in the first two columns, and the vertically integrated parameters including ice water path (IWP) and mean mass diameter (D_m) in the third and fourth columns. The BMCI-only and BMCI-OEM retrieval results from radiometer are compared with the triple-frequency radar retrievals. The fifth column shows the median absolute logarithmic error for each variable as a quantitative metric of retrieval accuracy. Each row corresponds to a different assumed ice cloud particle habit.

5 Summary and conclusions

This paper presents the hybrid BMCI-OEM retrieval algorithm that we developed for the submillimeter-wave radiometer in support of the PolSIR and PMM missions. The primary retrieval products are IWP and mean mass diameter D_m, but the algorithm can also retrieve the vertical profile of ice and liquid cloud microphysical quantities and water vapor. The Bayesian algorithm employs a hybrid strategy that combines the advantage of BMCI and OEM. With a comprehensive a priori database, the computationally efficient BMCI method can handle most retrieval cases. For cases where the BMCI fails due to insufficient database cases matching the TB observations, the additional OEM approach is applied to further maximize the posterior PDF and extend the retrieval capability beyond the limitations of the a priori database. With the CDF/EOF technique, the largely identical a priori constraints used in the BMCI step can be applied in the OEM step, even though the original a priori PDF in the retrieval database is highly non-Gaussian. The entire retrieval process therefore operates under reasonably consistent a priori constraints.

The hybrid BMCI-OEM algorithm is first applied to the evaluation of the C²OMODO radiometer using simulated observations. Retrieval results for IWP, LWP, IWV, and D_m are presented, and diagnostic analyses are performed based on the method provided by the OEM framework. Besides the column-integrated variables that are expected to become operational products for the PolSIR and PMM missions, we also explore potential vertical profiling capability using the averaging kernel matrix. The results indicate that the median vertical resolution achievable for ice cloud profiling is approximately 2.5 km. The DoF exhibits a positive dependence on IWP and IWV, and it typically ranges from 5 to 7 when the IWP exceeds 10 g m⁻². The information content analysis indicates that higher IWP values correspond to higher information content. The 325 GHz channels contribute approximately 4 bits of additional information for IWP in the presence of the 183 GHz channels. Although the 89 GHz channel shows little benefit for IWP retrievals, it provides substantial information for the LWP and IWV retrievals.

The BMCI-only and hybrid BMCI-OEM algorithms are then applied to real CoSSIR observations collected during the IMPACTS campaign to further evaluate its performance. The ice cloud parameters retrieved from CoSSIR observations are compared against triple-frequency radar retrievals obtained from collocated observations from the CRS and HIWRAP instruments, both of which are mounted on the same aircraft platform as CoSSIR. By upgrading the measurement vector and forward model, the hybrid BMCI-OEM algorithm is directly applied to the triple-frequency radar retrievals, highlighting the great extendibility of the hybrid Bayesian framework to other remote sensing measurements. The triple-frequency radar retrievals represent our best understanding of cloud profiles. Although we built a comprehensive database using all available CloudSat profiles within a latitude-longitude domain, the BMCI-only results still exhibit significant TB discrepancies, especially when TB is significantly depressed by dense ice clouds. The OEM method is demonstrated to significantly improve the TB agreements, with simulations for almost all channels closely following the observations when certain particle habits are used. The simulation performance is sensitive to the selection of particle habit, as some particle habits fail to reproduce the observed TB. Compared to the triple-frequency radar retrievals, the vertical extent of the BMCI-retrieved profiles tends to exceed the expected cloud layers at both cloud top and bottom, and the internal cloud structure is poorly resolved. The OEM method shows the ability to moderately remove spurious cloud layers and better highlight the cloud core. Statistical analysis indicates that the OEM provides moderate improvements over BMCI-only retrievals for both layer-resolved microphysical parameters and column-integrated variables. The OEM shows the advantage of reducing positive biases in the BMCI-only results and reducing errors. Also, the IWCs retrieved from radiometer observations generally follow the radar retrieved values, indicating the sensitivity of TB observations to layer-level ice cloud microphysics. However, the limited information in TB leads to retrievals varying by nearly an order of magnitude around the radar results. Using MedALE as the evaluation metric, the CoSSIR-retrieved IWC, D_me, IWP, and D_m differ from the radar retrievals by 3.75, 0.82, 2.53, and 0.55 dB, respectively, when applying the hybrid BMCI-OEM method. Compared with the BMCI-only retrievals, the errors are reduced by 0.21, 0.02, 0.85, and 0.13 dB for these four variables.

Several limitations of the current study should be noted. First, the radiative transfer model adopts a pencil-beam approximation and does not consider beamwidth or footprint size. As the scanning angle increases, the footprint size expands, and beam-filling effects become increasingly important. As a result, the pencil-beam assumption becomes less appropriate at larger scanning angles and may introduce systematic errors into the retrieval results. Second, this study assumes that the particles are completely randomly oriented and considers only the first Stokes component I. Since the target sensors provide polarization observations, the random-orientation assumption is insufficient for simulating polarized observations at off-nadir viewing angles. Even for the first Stokes component I, the random-orientation assumption may introduce errors when ice cloud particles exhibit preferential horizontal alignment. Third, the OEM currently relies on perturbation methods to calculate the Jacobian, which is computationally expensive. Alternative approaches for more efficient Jacobian calculations are therefore being investigated. Future work will focus on addressing these limitations to better support the upcoming instrument missions.

In summary, compared to the BMCI-only algorithm, the additional OEM procedure significantly reduces TB discrepancies while keeping the a priori information largely consistent with that used in the BMCI step. Despite the significant improvements in TB agreements, enhancements in the retrieval accuracy of ice cloud microphysical quantities are moderate due to the limited information content in the TB observations. The posterior PDFs of the OEM retrievals are still largely constrained by the a priori information. The hybrid Bayesian framework also demonstrates its great extendibility to other remote sensing observations. As more information becomes available through techniques such as multi-sensor integration or hyperspectral measurements, the hybrid Bayesian algorithmic framework shows increasing potential to better constrain cloud microphysics.