Optimal estimation of cloud properties from thermal infrared observations with a combination of deep learning and radiative transfer simulation

While traditional thermal infrared retrieval algorithms based on radiative transfer models (RTMs) could not effectively retrieve the cloud optical thickness of thick clouds, machine-learning-based algorithms were found to be able to provide reasonable estimations for both daytime and nighttime. Nevertheless, stand-alone machine learning algorithms are occasionally criticized for the lack of explicit physical processes. In this study, RTM simulations and a machine learning algorithm are synergistically utilized using the optimal estimation (OE) method to retrieve cloud properties from thermal infrared radiometry measured by the Moderate Resolution Imaging Spectroradiometer (MODIS). In the new algorithm, retrievals from a machine learning algorithm are used to provide a priori states for the iterative process of the OE method, and an RTM is used to create radiance lookup tables that are used in the iteration processes. Compared with stand-alone OE, the cloud properties retrieved by the new algorithm show an overall better performance by using statistical a priori information obtained by the machine learning algorithm. Compared with the stand-alone machine-learning-based algorithm, the radiances simulated based on retrievals from the new method align more closely with observations, and physical radiative processes are handled explicitly in the new algorithm. Therefore, the new method combines the advantages of RTM-based cloud retrieval methods and machine learning models. These findings highlight the potential for machine-learning-based algorithms to enhance the efficacy of conventional remote sensing techniques.

1 Introduction

Clouds play an important role in the Earth's energy budget by altering radiation patterns at both the surface and the top of the atmosphere (TOA) (Liou and Davies, 1993; Stubenrauch et al., 2006). Cloud properties change in response to variations in greenhouse gases, aerosol concentrations, and global surface temperature, leading to large uncertainties in climate change projections (Forster et al., 2021; Sassen et al., 2007). Grasping the variations in cloud properties is crucial for a comprehensive understanding of cloud dynamics and their radiative impacts on global climate change. The advancement of science and technology has positioned satellite remote sensing as a pivotal tool for monitoring cloud behaviors across diverse spatial and temporal scales. Active satellites like CloudSat and CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations) offer unparalleled cloud-profiling capabilities (Marchand et al., 2008; Sassen et al., 2009). Conversely passive satellites, known for their extensive swath observations, are widely utilized in a range of atmospheric research.

In recent decades, numerous efforts have been made to retrieve cloud properties using passive satellite instruments (Lai et al., 2019; Li et al., 2023; Min et al., 2020; Minnis et al., 2011; Poulsen et al., 2012; Tan et al., 2022; Zhao et al., 2012). A common method involves combining data from visible (VIS) and near-infrared (NIR) channels to construct lookup tables (LUTs) for daytime cloud microphysical properties, such as cloud optical thickness (COT) and cloud effective radius (CER) (Painemal and Zuidema, 2011; Twomey and Seton, 1980; Nakajima and King, 1990). This approach is grounded in the principle that cloud reflectance at non-absorbing VIS wavelengths predominantly depends on COT, while reflectance at absorbing NIR wavelengths is closely related to cloud effective radius (Arking and Childs, 1985; Rossow et al., 1989). Additionally, distinguishing liquid water from ice clouds using NIR channels (e.g., 1.65 µm) has also proven beneficial for deriving cloud-top height (CTH) (Harshvardhan et al., 2009; Menzel et al., 2008; Pilewskie and Twomey, 1987). Nonetheless, these VIS/NIR-based methodologies are confined to daytime operations, owing to their reliance on incident solar radiation, which is absent during nighttime hours.

Alternatively, nighttime cloud properties can be retrieved using thermal infrared (TIR) radiometry from passive satellites. Inoue (1985) employed the split-window method, leveraging brightness temperature (BT) and BT differences across various window channels, to derive COT and CER. Subsequently, numerous improvements and enhancements have been made to this method (Hamada and Nishi, 2010; Iwabuchi et al., 2018; Yang et al., 2005). Wang et al. (2016a) implemented an algorithm based on optimal estimation (OE) with Moderate Resolution Imaging Spectroradiometer (MODIS) infrared (IR) observations for cloud property retrieval, demonstrating the suitability of IR channels for thin-ice-cloud properties during both daytime and nighttime (Wang et al., 2016b). In addition, the CO₂-slicing method, which utilizes channels with a wavelength between 13–15 µm, is able to retrieve CTH effectively (Smith et al., 1974; Menzel et al., 1983). The atmospheric IR window measurements are also useful for CTH determination by comparing with the ambient atmospheric temperature profile (Garrett et al., 2009; Hong et al., 2007). However, IR window methods are less effective for optically thick clouds as their BT approaches asymptotic values (Garrett et al., 2009; Iwabuchi et al., 2016). While far-infrared channels (> 15 µm) are useful for ice clouds with substantial optical thickness (Libois and Blanchet, 2017), these far-infrared channels have been rarely measured by satellites during the last decade. Moreover, the retrieval methods based on plane-parallel cloud radiative transfer (RT) models face global application challenges due to their high computational demands (Wang et al., 2013).

Recently, machine learning techniques such as random forests, artificial neural networks, and deep learning have gained significant traction in remote sensing (Bai et al., 2021; Guo et al., 2022; Shi et al., 2020; Tan et al., 2023; Yuan et al., 2020; Zhao et al., 2023). Häkansson et al. (2018) used a neural network algorithm to retrieve cloud-top properties from several passive polar orbit sensors, greatly improving CTH retrievals. Advanced machine learning algorithms have particularly enhanced CTH retrievals for high and thin clouds (Min et al., 2020). Wang et al. (2022) developed a convolutional neural network (CNN)-based framework (TIR-CNN), which utilizes TIR radiometry from MODIS to retrieve COT, CER, and CTH. This framework demonstrates consistency with both passive and active cloud products and is effective during both daytime and nighttime (Wang et al., 2022, 2023). Tana et al. (2023) obtained cloud detection and cloud microphysical properties with high spatial and temporal resolutions from TIR spectral channels of Himawari-8 using a machine learning algorithm. Zhao et al. (2023) applied a deep learning ResUnet model for retrieving cloud phase (CLP), COT, CER, and CTH using Fengyun-4A satellite observations.

However, the reliance of these machine learning methods on mathematical and statistical approaches typically leads to an implicit assimilation of the relationships between cloud properties and radiance observations, making it difficult to trace how input data are transformed into output predictions. A great number of cloud property users favor remote sensing products that offer explicit physical interpretations. Therefore, enhancing traditional inversion algorithms with machine learning algorithms can be beneficial.

In this study, we combine traditional radiative transfer simulations with TIR-CNN retrievals using the OE method (OE-CNN-IR) to retrieve COT, CER, and CTH from MODIS, which is effective under both daytime and nighttime conditions. The Community Radiative Transfer Model (CRTM) is utilized to simulate MODIS IR observations and generate LUTs for cloud properties. The TIR-CNN retrievals are employed as a priori states, and an iterative process based on gradient descent is performed to get an optimal estimation. The performance of the proposed OE-CNN-IR model is subsequently compared with a stand-alone OE method utilizing fixed a priori states. Details of the data and the enhanced OE-CNN-IR method are presented in Sect. 2. Section 3 outlines the retrieval results and their evaluation against cloud products from passive and active sensors. Conclusions are summarized in Sect. 4.

2 Data and methodology

2.1 Data

2.1.1 MODIS data

This study utilizes global data observed by the MODIS instrument aboard the Aqua satellite. Aqua-MODIS continuously monitors the Earth–atmosphere system with 36 spectral bands ranging from 0.4 to 14.4 µm. The Aqua-MODIS official Collection 6.1 (C6.1) products (MYD021KM, MYD03, MYD35, MYD06, and MCD12C1), available at https://ladsweb.modaps.eosdis.nasa.gov/search/ (last access: 13 May 2024), have been selected for this study. These products, with spatial resolutions of 1 and 5 km, are chosen for their widely accepted quality (Wang and Christopher, 2003). In this study, the TIR radiation from the Aqua-MODIS Collection 6.1 (C6.1) Level 1B calibrated radiance products (MYD021KM) are converted to BTs using the Planck function. All channels with wavelengths greater than 6.5 µm are used, except that the 30th channel (primarily used for ozone retrievals) is not used to reduce uncertainties induced by ozone. Additionally, atmospheric parameters from MYD03 and MYD06, including surface temperature, land surface type, and cloud phase, are used as ancillary data for LUT construction and forward radiative simulations. The product in Table 1, reported in “Cloud Phase Optical Properties” is the daytime-only phase used in the MYD06 cloud optical retrievals, and “Cloud Phase Infrared” is a daytime and nighttime product derived from three IR window channel pairs. “Cloud Phase Optical Properties” is used in daytime to determine cloud types, while “Cloud Phase Infrared” is used in nighttime only in our paper. Cloud optical and physical parameters such as COT, CER, and CTH from MYD06 serve to verify the accuracy of daytime retrievals. All parameters are aligned to a 5 km spatial resolution grid, ensuring uniformity in data and variables. To validate the applicability of the inversion algorithm, retrievals are compared to MYD06 data from a single day each month in 2009 (1 January, 10 February, 12 March, 11 April, 11 May, 10 June, 10 July, 9 August, 8 September, 8 October, 7 November, and 7 December. The default spacing between adjacent days is 30, and the spacing is set to 40 if a date lies in the same month as the previous date), capturing the variability of atmospheric conditions throughout a seasonal cycle and facilitating a comprehensive evaluation across different scenarios. The total sample size of MYD06 for comparison is ∼ 4.7 million. The retrievals are performed over the whole globe, and the data between 60° S–60° N are used in the validation. Table 1 summarizes the data and parameters used in our retrieval model.

Table 1 Summary of MODIS data sources and preprocessing parameters.

2.1.2 Lidar–radar detection cloud products

Cloud-Aerosol LIdar with Orthogonal Polarization (CALIOP), a space-based lidar instrument aboard the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite, provides vertical profiles of clouds and aerosols in Earth's atmosphere. CALIOP can perform observations at both daytime and nighttime, overcoming the limitations of passive optical instruments, but it can not penetrate thick clouds. The cloud-profiling radar (CPR) aboard CloudSat is a radar system that sends out microwave pulses and measures the reflected energy from clouds. This technique is particularly suited to determining the structure and ice content within clouds but fails to detect thin clouds. The DARDAR product (Delanoë and Hogan, 2010), integrating data from both CALIOP and CPR, offers a comprehensive atmospheric column view that neither instrument can achieve independently. This extensive dataset includes information on cloud-top and base heights, optical thickness, ice content, and aerosol layers. In our study, the ice cloud product of DARDAR in 2009 is used to evaluate the inversion results during both daytime and nighttime conditions, and ∼ 0.54 million pixels are colocated in the comparison processes.

2.2 Development of the retrieval algorithm

The core algorithm of our inversion method is the optimal estimation method, which utilizes the CRTM as the forward model and incorporates CNN results as a prior information. Figure 1 illustrates the architecture of our retrieval models. Initially, temperature, humidity, and ozone from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA5) (Hersbach et al., 2020) are used to construct lookup tables for each 0.25° × 0.25°spatial grid box. These LUTs enumerate the BT for each channel, corresponding to varying COT, CTH, and CER. Subsequently, the OE method is performed to retrieve cloud properties. The OE method can get the optimal solution by accounting for all spectral information. However, the iteration may have started a long way from the solution in a nonlinear problem, and the cost function decrease is much slower. Starting with a better first guess rather than climatology value can make the process converges much more quickly (Rodgers, 2000). The deep learning methods can achieve high accuracy, and once trained, they offer very fast prediction speeds. However, due to the nature of neural networks, deep learning results often lack interpretability, leading to the perception of deep learning as a black box model. In the OE-CNN-IR approach, the TIR-CNN-derived cloud properties provide a priori states for iterative processes, which is subsequently refined through iterative minimization of the objective cost function, while the climatology values were used as starting points in OE-IR. This method iteratively adjusts parameters to reconcile radiative transfer simulations with observed data. Further details are presented below.

Figure 1 The architecture of the retrieval model. (a) The establishment of look-up table. (b) The iteration steps in the optimal estimation process.

2.2.1 Forward model

The CRTM, developed by the U.S. Joint Center for Satellite Data Assimilation (JCSDA), spans a broad spectrum of channels from visible to microwave. It is widely used in simulating radiances at the top of the atmosphere for various satellite sensors, owing to its flexible interface, sophisticated radiative transfer processes, and efficient numerical computation (Han et al., 2006). The model divides the atmosphere into a series of vertical layers, and the temperature, pressure, and composition of each layer is assumed to be homogenous. CRTM solves the radiative transfer equations throughout the atmosphere. Its precision and reliability have been extensively validated by ground-based and satellite observations (Zou et al., 2016). Considering that the optical properties of ice cloud crystals in the CRTM and MODIS products differ due to variations in particle habit assumptions (Yi et al., 2016; Yao et al., 2018), the volumetric extinction cross section in the CRTM is adjusted by a scaling factor of 0.4, resulting in simulated brightness temperatures that are consistent with observations (Fig. 6a, b, c). For each grid cell, the CRTM simulates TIR radiances corresponding to a set of COT, CER, and CTH values at each location, from which a LUT is subsequently constructed. Table 2 provides a detailed list of the cloud properties and ancillary parameters used in these calculations.

Table 2 Geometries and cloud properties selected to calculate the cloud lookup tables.

The outputs of the forward model (F) can be expressed as a function of cloud properties and ancillary parameters:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{Y}& =[{\mathrm{BT}}_{\mathrm{1}},{\mathrm{BT}}_{\mathrm{2}},\mathrm{\dots },{\mathrm{BT}}_m{]}^T\\ & =\mathit{F}\left[\mathit{X}\right(\mathrm{COT},\mathrm{CER},\mathrm{CTH}),\mathit{P}]+\mathit{e},\end{array}\end{array}$$

where Y is a vector consisting of m MODIS IR observations of BT, and P is a vector encompassing various ancillary variables, including air temperature, water vapor concentration, ozone concentration profiles, surface emissivity spectrum, and surface temperature. The state vector (X) includes COT, CER, and CTH, and e is an error term.

Figure 2 depicts the variation in CRTM output (F) (expressed in BT) as a function of ice cloud properties, derived from a simulation using the atmospheric profile dated 10 June 2009, at 00:00 UTC, at coordinates 175.87° E and 60.55° N. In Fig. 2a, with fixed CER and CTH, the TOA BTs in MODIS IR bands generally decrease with increasing COT. Notably, for COT > 10, the slopes approach zero, causing challenges in inversion accuracy. In the case of fixed COT and CTH (Fig. 2b), TOA BTs decrease with increasing CER values when CER is below 10 µm across all channels, followed by minor oscillations in most channels, except that band 29 shows significant variations. In Fig. 2c, for CTH values under 11 km, TOA BTs are negatively correlated with CTH noticeably.

Figure 3 shows the relationship between TOA BTs and liquid cloud properties, which reveals a weaker response to changes in COT and CER compared to ice clouds. Nevertheless, TOA BTs decrease noticeably with increasing liquid cloud CTH. In summary, CTH is the most accurately determinable variable for both ice and liquid clouds due to the high sensitivity of TOA BTs to CTH. For ice clouds, COT values below 10 generally allow for more accurate retrieval of cloud properties in theory. However, retrieving CER for ice clouds poses greater challenges due to the complexity of ice particle size distribution and shape. Liquid clouds, conversely, show no strong sensitivity of TOA BTs to both COT and CER, and it is more difficult to accurately retrieve these cloud properties solely based on TOA BT observations.

Figure 2 Radiative transfer model simulations for ice clouds. The atmospheric profile is from the coordinates with a longitude of 175.87° E and a latitude of 60.55° N, on 10 June 2009, at 00:00 UTC. (a) TOA BTs as a function of COT, when CER and CTH are set to 20 µm and 10 km, respectively. (b) BT as a function of CER, when COT and CTH are set to 5 and 10 km, respectively. (c) BT as a function of CTH, when COT and CER are set to 5 and 20 µm, respectively.

Figure 3 Same as Fig. 2 but for liquid clouds.

2.2.2 The convolutional neural network infrared method

The convolutional neural network using thermal infrared (TIR-CNN) model is trained with solar-independent variables (thermal infrared radiances, viewing zenith angles, and altitude) as inputs and uses standard MYD06 products (COT, CER, and CTH in the daytime) as targets. To capture a comprehensive range of the Earth's surface and viewing geometries while accounting for seasonal variations, Wang et al. (2022) collected full-year granules from 2010 to create the training dataset. Products with a 10 d interval from 2011 were selected as the validation dataset during the training phase. Additionally, the 10 d interval data from 2009, which are independent of both the training and validation datasets, served as the testing dataset. The granules were divided into samples sized 256 × 256 km. After preprocessing, there are 1 888 680 samples in the training dataset, 191 520 in the validation dataset, and 382 760 in the testing dataset. This TIR-CNN model is an asymmetric architecture, featuring an equal number of encoding and decoding layers arranged sequentially. The basic convolutional block consists of two 2D convolutional layers with 3 × 3 kernels. Each convolutional layer is followed by a batch normalization layer and a leaky rectified linear unit (Wang et al., 2022). Through training, the model can capture context and learn the complex nonlinear relationship between the input variables and targets, which can be applied in the cloud property retrievals during both daytime and nighttime. The convolutions in the TIR-CNN model are beneficial in considering statistic information from neighboring fields in training. Theoretically, spatial distributions and optical and microphysical properties of clouds are all determined by the meteorological backgrounds, so cloud properties are statistically connected to their horizontal distributions. In addition, the effective radius of ice cloud particles are functions of cloud temperature. The CNN-based deep learning architecture is able to capture the statistical features among adjacent pixels of satellite observations as a solution for retrieving cloud optical and micro-physical properties (Wang et al., 2022, 2023), so it is able to provide more information than traditional algorithms that retrieve cloud properties from infrared radiances of individual pixels. The benefits of machine learning in IR cloud retrievals have also been demonstrated independently by the results of Tana et al. (2023) and Zhao et al. (2023).

2.2.3 Optimal-estimation-based retrieval method

The OE-based retrieval method, as introduced by Rodgers (2000), is designed to derive the best estimates of atmospheric quantities (such as temperature, humidity, aerosol concentration, or trace gas concentrations) by minimizing the discrepancy between observed measurements and the model simulations. This method combines information from both the measurement data and a priori knowledge, typically obtained from atmospheric models or ancillary data sources. A key strength of the OE method is its capability in addressing complex atmospheric retrieval challenges, enabling simultaneous retrieval of multiple parameters in contexts where physical processes are nonlinear and highly coupled. It provides a rigorous and statistically robust method to estimate atmospheric parameters, along with quantifying the associated uncertainties.

The OE method aims to identify the most probable state variables by minimizing a cost function J (Wang et al., 2016a):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathit{J}& ={\left[\mathit{F}\left(\mathit{X},\mathit{P}\right)-\mathit{Y}\right]}^T{\mathbf{S}}_y^{-\mathrm{1}}\left[\mathit{F}\left(\mathit{X},\mathit{P}\right)-\mathit{Y}\right]\\ & +{\left[\mathit{X}-{\mathit{X}}_a\right]}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left[\mathit{X}-{\mathit{X}}_{\mathrm{a}}\right],\end{array}\end{array}$$

where X_a and X are the prior and posterior state vectors, respectively. S_y and S_a are the covariance matrices of the observation-to-simulation differences and the uncertainty of the prior state vector, respectively. When the uncertainties of the a priori state are large (e. g., OE-IR), the cost function J is primarily influenced by the first term in OE-IR. If the uncertainties of the a priori state are small, then the second term is also important in the iteration process (e. g., OE-CNN-IR). Then we employ an iterative process to find an optimal solution based on observed data and a priori states. In this paper, the gradient descent method for the (i+1)th iteration is implemented by

$$\begin{array}{}\text{(3)}& {\mathit{X}}_{i+\mathrm{1},n}={\mathit{X}}_{i,n}-\mathit{\theta }{\displaystyle \frac{\partial {\mathit{J}}_{i,n}}{\partial X_{i,n}}},\end{array}$$

where

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial {\mathit{J}}_{i,n}}{\partial {\mathit{X}}_{i,n}}}={\displaystyle \frac{\mathit{J}\left({\mathit{X}}_{i,n}+\mathit{\delta }{\mathit{x}}_n\right)-\mathit{J}\left({\mathit{X}}_{i,n}\right)}{\mathit{\delta }{x}_n}},\end{array}$$

and where θ represents a learning rate, n represents the nth cloud parameter (COT, CER, and CTH), and θ is set to be the same for all three variables. In this paper, θ is initially set to 0.05 for the first 200 iterations, and after those initial 200 iterations, the learning rate is then reduced to 0.01 for the subsequent 100 iterations. δx_n represents a small increase in the n cloud parameter, and $\mathit{J}\left({\mathit{X}}_{i,n}+\mathit{\delta }{\mathit{x}}_n\right)$ is calculated using LUTs.

In this paper, the stand-alone OE-IR method relies on a fixed a priori value for its iterative process, whereas OE-CNN-IR utilizes results from TIR-CNN as its a priori states for further refinement. These methods illustrate the integration of traditional optical estimation techniques with advanced machine learning models to potentially enhance the accuracy and reliability of atmospheric measurements.

Figure 4 shows the iterative variations in cost function, COT, CER, and CTH for both OE-IR and OE-CNN-IR under various conditions. Both algorithms show a significant decrease with an increasing number of iterations, and all iterations can achieve successful convergence. However, the initial value for OE-CNN-IR is lower than OE-IR. For smaller COT values (COT < 10), OE-CNN-IR and OE-IR exhibit consistent effects on COT and CTH, converging with close values. For both methods, the CER depends on a priori states. A distinct difference is that OE-CNN-IR starts with a significantly lower cost than OE-IR and maintains more stability throughout the iteration process. For larger COT values (COT > 10), the CTH of these two methods converge to same value, despite differences in a priori states. The COT of OE-IR struggles to iterate towards the expected target during the iteration process, while OE-CNN-IR maintains stably around the a priori values. The iterative results indicate that both methods perform well on COT and CTH for small COT values. However, for large COT values, the OE-IR method is unable to produce accurate results under these conditions. In contrast, the OE-CNN-IR is able to retrieve COT for thick clouds effectively.

Figure 4 The change of cost function and cloud parameters in the iteration processes, with OE-IR in red and OE-CNN-IR in blue. Left panels are for an illustrative ice cloud layer with a small optical thickness case, and right panels are for a large optical thickness case.

2.3 Metrics for performance evaluation

In this study, the magnitude of calculation errors, systematic bias errors, and linear correlation between outputs and standard values are quantitatively assessed using three key statistical metrics: root mean squared error (RMSE), mean bias error (MBE), and the Pearson correlation coefficient (r).

$$\begin{array}{}\text{(5)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\left(y_i-f_i\right)}^{\mathrm{2}}}{N}}},\end{array}$$

$$\begin{array}{}\text{(6)}& \mathrm{MBE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{f}_i-y_i}{N}},\end{array}$$

$$\begin{array}{}\text{(7)}& r={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left(f_i-\stackrel{\mathrm{‾}}{f}\right)\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}{\sqrt{\left({\sum }_{i=\mathrm{1}}^N{\left(f_i-\stackrel{\mathrm{‾}}{f}\right)}^{\mathrm{2}}\right)\sqrt{\left({\sum }_{i=\mathrm{1}}^N{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}\right)}}}},\end{array}$$

where N is the total number of calculated points; y and f are the true and estimated values, respectively.

3 Results and discussion

3.1 Case studies of OE-CNN-IR and OE-IR retrievals

To illustrate the daytime efficacy of the proposed method, a granule from Aqua-MODIS, captured at 03:00 UTC on 10 June 2009, has been chosen. This particular granule spans the southwestern Pacific Ocean, encompassing the geographical region from 0 to 20° S latitude and from 150 to 175° E longitude, as depicted in Fig. 5.

Figure 5a–c shows the spatial distribution of BT for each respective channel. These measurements reveal variations in thermal radiation, which are correlated with cloud properties specific to the wavelengths of the channels used. Figure 5d–f shows the cloud physical properties as derived using standard MODIS retrieval algorithms. The COT, CER, and CTH from the MYD06 product provide a benchmark for comparison with other inversion methods. The analysis of BT from channels 29, 31, and 32 shows a clear negative correlation with both COT and CTH, and regions with higher BT typically correspond to clouds with lower optical thickness and lower cloud-top heights. This is in line with the principle that thinner clouds permit more infrared radiation to escape from the Earth's surface and atmosphere, leading to higher observed BT. Furthermore, the analysis indicates that clouds with higher BT generally have lower altitudes. The patterns in Fig. 5, which display cloud properties derived from various inversion techniques, corroborate the physical relationships illustrated in Figs. 2 and 3. Figure 5g–i presents the retrieval results from the deep learning algorithm TIR-CNN method. The CNN-derived retrievals are not only consistent with MYD06 products in spatial patterns, but also agree well with the magnitudes of results. Figure 5j–l presents the retrieval results from the OE-CNN-IR method, showing similar spatial distributions to the standard MYD06 products for COT and CTH. However, significant differences are noted in CER retrieved by OE-CNN-IR and MYD06 products. This finding aligns with the work of Wang et al. (2016), which highlighted substantial discrepancies in CER retrieval when using OE-IR versus MYD06 products. Specifically, the C6 MYD06 cloud particle size information presented here is inferred from the 2.1 µm reflectance, which may capture signals reflected from the lower parts of a cloud (Zhang et al., 2009).

Figure 5m–o shows the results retrieved using the traditional OE-IR method with climatological a priori states, which employs climatological values of COT, CER, and CTH as a priori states for OE iteration. In the case where COT values are below 10, the OE-IR COT closely matches MYD06 products, indicating that it is able to capture the COT of thinner clouds. However, the inability of OE-IR to retrieve COT values greater than 10 suggests a limitation in the technique's sensitivity to optically thicker clouds, aligning with findings from Wang et al. (2016). This threshold effect arises from the TIR BT independent of COT in thick clouds (as shown in Fig. 2a). The performance of CTH retrievals using the OE-IR method is comparable to that of the OE-CNN-IR method, while the inversion of CER is not very effective due to limitations in the physical mechanisms.

Figure 5 Comparison of cloud properties obtained from the TIR-CNN model, OE-CNN-IR model, OE-IR model, and standard MODIS products for an illustrative daytime granule on 10 June 2009 (03:00 UTC). Panels (a, b, c) are BT images of MODIS band 29, 31, and 32, respectively. Panels (d, e, f) are the COT, CER, and CTH from the MYD06 product, respectively. Panels (g, h, i) are the COT, CER, and CTH from the TIR-CNN model, respectively. Panels (j, k, l) are the COT, CER, and CTH from the OE-CNN-IR model, respectively. Panels (m, n, o) are the COT, CER, and CTH from the OE-IR model, respectively.

Figure 6 compares MODIS-observed BTs with simulated BTs derived from CRTM. Figure 6a–c shows CRTM-simulated radiances using baseline MODIS cloud products, serving as a control scenario for comparative analysis. The correlation coefficients for channels 29, 31, and 32 are 0.877, 0.905, and 0.891, respectively, indicating CRTM's proficiency in simulating MODIS cloud products. However, there is a persistent negative MBE across these channels. Figure 6d–f presents a comparison between observations and BT simulated by CRTM and TIR-CNN retrievals, with outcomes that are analogous to those depicted in Fig. 6a–c. The pronounced correlation indicates that CNN-based inputs proficiently replicate the spatial and radiometric features of clouds, showing high consistency with MODIS MYD06 products. When OE-CNN-IR or OE-IR cloud property retrievals are used to simulate BT, the correlation coefficients between the simulated BT and observations increase significantly, and the absolute values of MBE and RMSE decrease significantly (as shown in Fig. 6g–l). The OE-CNN-IR model incorporates the OE iterations, which reduce the discrepancy between simulated and observed BT. The results indicate that retrievals of the OE-CNN-IR method align more closely with BT observations compared to the stand-alone TIR-CNN model.

Figure 6 Comparison between MODIS BT observations and simulated BT based on MODIS cloud properties (a, b, c), TIR-CNN a priori inputs (d, e, f), OE-CNN-IR estimations (g, h, i), and OE-IR estimations (j, k, l), based on an illustrative granule of 10 June 2009 (Fig. 5). The first column is the comparison between the simulation and observation of band 29, the middle column is for band 31, and the right column is for band 32. The unit of BT is kelvin (K).

Figure 7 illustrates a nighttime case of cloud parameter retrievals using TIR-CNN, OE-CNN-IR, and OE-IR methods. The data for this analysis are sourced from a randomly selected granule captured on 10 February 2009, at 21:00 UTC. Figure 7a–c displays the BTs at channels 29, 31, and 32, and Fig. 7d–f shows the COT, CER, and CTH retrieved by the TIR-CNN algorithm. The relationship between COT and CTH with BT at night is generally consistent with that during the day. Figure 7g–i shows the COT, CER, and CTH retrieved by the OE-CNN-IR algorithm. The OE-CNN-IR retrievals align well with the high- and low-value areas in the BT images, indicating that OE-CNN-IR effectively discerns the intricate spatial variations in cloud properties during nighttime conditions. Figure 7j–l displays the cloud parameters retrieved using the OE-IR method. In this analysis, the predominance of values falls below 10, which signifies a more constrained retrieval scope when contrasted with the OE-CNN-IR method. Nevertheless, the distribution of CTH derived from OE-IR closely mirrors that obtained from OE-CNN-IR, affirming its dependability for estimating the top height of clouds. Additionally, both methods exhibit comparable distributions in CER.

Figure 7 Comparison of cloud properties obtained from the OE model and standard MODIS products for an illustrative nighttime granule on 10 February 2009 (21:00 UTC). (a) BT image of MODIS band 29. (b) BT image of MODIS band 31. (c) BT image of MODIS band 32. Panels (d, e, f) are the COT, CER, and CTH from the TIR-CNN, respectively. Panels (g, h, i) are the COT, CER, and CTH from the OE-CNN-IR model, respectively. Panels (j, k, l) are the COT, CER, and CTH from the OE-IR model, respectively.

Figure 8 Scatterplots of the pixel-level comparisons between the retrievals and MYD06 products for ice clouds over oceans. (a, d, g) Pixel-by-pixel comparisons of COT, CER, and CTH from OE-CNN-IR with the MYD06 ice cloud products over ocean in 2009. (b, e, h) Scatterplots of the pixel-level comparisons between the MYD06 cloud products and OE-IR comparable retrievals. (c, f, i) The probability density functions obtained from MYD06 products, OE-CNN-IR-, and OE-IR-derived results. Color shadings denote the number of observations in each respective pixel. All comparable retrievals are constrained to cases with SZA < 60°and latitudes between 60° S and 60° N.

3.2 Comparison between retrievals and MYD06 products in the daytime

Figure 8 presents scatterplots that provides pixel-level comparisons of cloud property retrievals from OE-CNN-IR and OE-IR against the MYD06 ice cloud products over ocean for 2009. The left column of Fig. 8 offers a detailed pixel-by-pixel comparison for COT, CER, and CTH between OE-CNN-IR and the MYD06 ice cloud products. The middle column displays comparisons between MYD06 cloud products and OE-IR retrievals. The right column displays the probability density functions obtained from MYD06 products, OE-CNN-IR-, and OE-IR-derived results. The color scale in these plots indicates the number of observations in each grid, visually representing data point density. Due to the large uncertainties of MODIS in retrieving COT in polar regions, retrieval constraints have been established. These include limiting the solar zenith angle (SZA) to less than 60° and restricting latitudes to between 60° S and 60° N, thereby ensuring consistency and reliability in these comparisons. In Fig. 8a, the correlation coefficient between OE-CNN-IR COT and MYD06 COT is 0.835, indicating a strong positive correlation. In comparison, OE-IR achieves a COT correlation coefficient of 0.667 against MODIS products, indicating a slightly weaker relationship than that reported in Wang et al. (2016). In Fig. 8c, the distributions provided by MYD06, OE-CNN-IR, and OE-IR are relatively similar for COT values less than 10. The OE-CNN-IR retrievals contain a lot of cases with COT > 15, which is consistent with MODIS, but OE-IR retrievals do not contain clouds with COT > 15. The underestimation of COT for thick clouds by OE-IR is consistent with Wang et al. (2016). Therefore, it is concluded that both OE-CNN-IR and OE-IR show consistent performance for COT below 10, but OE-CNN-IR performs much better for thicker clouds. With respect to CER, both algorithms demonstrate moderate to weak correlation coefficients, reflecting the inherent physical constraints of the retrieval process. Nonetheless, OE-CNN-IR outperforms OE-IR with a correlation coefficient of 0.794, suggesting enhanced performance. In Fig. 8f, the results from OE-IR appear to be concentrated around an a priori value of 30 µm, whereas the results from OE-CNN-IR maintain a distribution that is more similar to that of MYD06. For the CTH retrieval, both OE-CNN-IR and OE-IR demonstrate good performance, with correlation coefficients of 0.871 and 0.808, respectively. Overall, the statistical analysis in Fig. 8 underscores the retrieval capability of OE-CNN-IR, particularly for COT and CER, compared to stand-alone OE-IR.

Figure 9 expands the ice cloud analysis from Fig. 8 to encompass liquid and ice clouds over both land and ocean, offering a more comprehensive evaluation of the retrieval algorithms across varied cloud conditions. In the case of liquid clouds above 10, the BT is not sensitive to COT, leading to most OE-IR COT retrievals clustering around a value of 10. This indicates difficulties in effectively retrieving COT for liquid clouds, so the OE-IR method has been used to retrieve cloud properties of ice clouds only. In contrast, the performance of OE-CNN-IR is much better. This shows that OE-CNN-IR can be improved by using TIR-CNN outputs as a priori state, allowing for accurate retrievals even in situations of lower BT sensitivity, as observed in liquid clouds. Regarding CER, the gradient of CER with respect to BT of liquid clouds tends toward zero. These artifacts signal the limitations of the retrieval algorithm under minimal BT gradient conditions. Despite these challenges for CER, both OE-CNN-IR and OE-IR perform exceptionally well in retrieving CTH, with correlation coefficients of 0.913 and 0.931, respectively. These high correlations reflect the algorithms' effectiveness in estimating CTH.

Figure 9 Same as Fig. 8 but includes liquid clouds over ocean, as well as ice and liquid clouds over land.

3.3 Comparison with products from active sensors

Under nighttime conditions, where standard MYD06 cloud products do not offer cloud optical properties, the evaluation is supplemented by incorporating near-real-time data from active sensors (DARDAR, derived from CloudSat/CALIPSO observations). Spatially and temporally colocated samples from 2009 are employed to evaluate the performance and generalization capabilities of the OE model during night conditions. These criteria are applied to achieve the closest possible data correspondence between the two different instruments, facilitating a meaningful assessment of the OE model's nighttime performance.

Figure 10 presents a detailed comparison of COT retrievals for ice clouds using OE-CNN-IR and OE-IR methods, benchmarked against DARDAR cloud products. The comparisons are confined to latitudes between 60° N and 60° S to ensure a comprehensive assessment across both daytime and nighttime conditions. The daytime correlation coefficient for OE-CNN-IR versus DARDAR COT is 0.651, with a slightly lower nighttime correlation of 0.583. These values are similar to the correlation between MYD06 COT and DARDAR COT 0.647 (Wang et al., 2022). In contrast, OE-IR exhibits lower correlation coefficients, with 0.546 during the day and 0.503 at night. Nevertheless, the RMSE of OE-IR is lower than that of OE-CNN-IR. Notably, the OE-CNN-IR method demonstrates better performance for COT > 10. This suggests that OE-CNN-IR is more adept at capturing the variability of thicker ice clouds, which is important for understanding cloud radiative effects and their implications for weather and climate systems.

Figure 10 Comparisons of OE-CNN-IR COT, OE-IR COT and DARDAR COT for ice clouds over oceans. (a–c) are comparisons between OE-CNN-IR COT and DARDAR products. (d–f) are comparisons between OE-IR COT and DARDAR products. The left column is for daytime comparisons, the middle column is for nighttime comparisons, and the right column is for all-day comparisons.

4 Conclusions

This study introduces a cloud property retrieval method based on optimal estimation (OE-CNN-IR), which integrates traditional radiative transfer simulations with a machine learning method. Designed for retrieving COT, CER, and CTH, this method is applicable for passive satellite imagery under both daytime and nighttime conditions. Retrievals from a machine learning algorithm (TIR-CNN) are used to provide a priori values for OE iteration, and a radiative transfer model (RTM) is used to create radiance lookup tables that are used in the iteration processes. Subsequently, the retrievals are iteratively adjusted to minimize discrepancies between the IR observations and RTM simulations. The efficacy of OE-CNN-IR is validated against MYD06 products and active sensor cloud products, and the results are compared to a stand-alone optimal estimation model (OE-IR).

The validation results reveal that the OE-CNN-IR method outperforms the stand-alone OE-IR model, especially for cloud optical thickness of thick clouds, due to its better a priori values. Correlation coefficients with MYD06 products have exhibited marked improvements: correlation coefficients for COT increase from 0.667 to 0.835, correlation coefficients for CER increase from 0.348 to 0.794, and correlation coefficients for CTH increase from 0.808 to 0.871. In nighttime evaluations, the OE-CNN-IR method consistently outperforms the traditional OE model when compared with DARDAR COT. The consistency between OE-CNN-IR retrievals and MYD06 products is as good as that of the stand-alone machine learning retrieval algorithm (i.e., TIR-CNN), and the radiance simulations based on OE-CNN-IR retrievals exhibit greater consistency with actual observations, as depicted in Fig. 6. Furthermore, the algorithm explicitly addresses physical processes, aligning with the preferences of scientists who advocate for physically based methodologies. While the OE-CNN-IR method in this study is primarily applied to Aqua-MODIS imagery, it can be potentially applied to other sensors with similar infrared (IR) channels. For instance, it can be readily adapted to geostationary satellites, given their analogous wavelength ranges (Tana et al., 2023; Zhao et al., 2023).