Optimizing cloud optical parameterizations in Radiative Transfer for TOVS (RTTOV v12.3) for data assimilation of satellite visible reflectance data: an assessment using observed and synthetic images

Radiative Transfer for TOVS (RTTOV) is a commonly used forward-operator software package for the data assimilation (DA) of satellite visible reflectance data. However, the wide choice of cloud optical parameterizations (COPs) in RTTOV poses challenges in discerning the optimal configuration. In this study, the performance of different COPs was evaluated by comparing the observed and synthetic visible satellite images. Observed images (O) were provided by Fengyun-4B (FY-4B) and Himawari-9, two operational geostationary meteorological satellites covering East Asia. Synthetic images (B) were generated by RTTOV (v12.3) with the discrete ordinate method (DOM) and the Method for FAst Satellite Image Simulation (MFASIS). The inputs to RTTOV were provided by the 3 h forecasts of the China Meteorological Administration Mesoscale (CMA-MESO) model and the fifth-generation European Centre for Medium-Range Weather Forecasts reanalysis (ERA5) data. On average for the domain, B was smaller than O, especially in cloudy situations. The minimum O−B bias was revealed for the COP of liquid water clouds in terms of effective diameter (D_eff) in combination with the COP of ice clouds developed by the Space Science and Engineering Center (SSEC), with the D_eff for ice clouds parameterized in terms of ice water content and temperature. Compared with the O−B biases, the standard deviations of the O−B departure were less sensitive to COPs. In addition, histogram analysis of reflectance indicated that the synthetic images with the minimum O−B bias resembled the observed images best. Therefore, the optimal cloud optical parameterization was proposed to be the “D_eff + SSEC” suite.

1 Introduction

Fengyun-4A (FY-4A) and Fengyun-4B (FY-4B) are China's two new-generation geostationary meteorological satellites that form a dual-satellite observation constellation. The two satellites carry the Advanced Geostationary Radiation Imager (AGRI), which provides radiance observations ranging from visible to infrared bands over East Asia and Western Pacific areas. From 1 February to 5 March 2024, FY-4B replaced FY-4A and started its operational observations from 00:00 UTC on 5 March 2024. Its spectral observations contain vital information on cloud, precipitation, aerosols, underlying surfaces, temperature, and humidity (Yang et al., 2017). Compared with infrared radiation, visible radiation is capable of penetrating deeper into the cloud and thereby provides more information about them. In addition, visible radiation is less sensitive to humidity and temperature, which facilitates the remote sensing of lower-layer clouds. Therefore, visible radiance data are receiving attention from the data assimilation (DA) community (Scheck et al., 2020; Schröttle et al., 2020; Zhou et al., 2022, 2023; Kugler et al., 2023).

For the DA of satellite radiance data, it is necessary to convert the model state into simulated observations to compute innovations. The accuracy of the forward operators influences the analysis fields and the subsequent forecasting results in two aspects. On one hand, forward operators influence the observation increments, which are converted into the analysis increments in the state variable space (e.g., Anderson, 2001). On the other hand, forward operators influence the bias correction of satellite radiance data because the bias correction is usually based on observation (O) minus background (B) analyses, where B is simulated by a forward operator based on the forecasts of a numerical weather prediction (NWP) model (Harnisch et al., 2016; Noh et al., 2023; Zhou et al., 2024). DA of FY-4A visible radiance data has been explored under the framework of observing system simulation experiments (OSSEs) (Zhou et al., 2022, 2023). For the experiment designs of OSSEs, the model equivalents of observation were simulated by the Radiative Transfer for TOVS (RTTOV) software using the discrete ordinate method (DOM) (RTTOV-DOM hereafter). Meanwhile, the background model state was converted into the simulated observations by the same forward operator. Therefore, errors due to the forward operator were neglected. However, to extend the DA of FY-4A and FY-4B visible reflectance to real-world cases, the performance of the forward operator must be thoroughly evaluated.

Recently, the O−B statistics of FY-4A visible band (0.55–0.75 µm) were explored using RTTOV-DOM based on the short-term forecasts of the China Meteorological Administration Mesoscale (CMA-MESO) (Zhou et al., 2024). One of the findings of Zhou et al. (2024) is that the synthetic images were more reliable in cloud-free situations than cloudy situations. In cloudy situations, an important factor constraining the performance of the RTTOV is the cloud optical parameterizations. In Zhou et al. (2024), cloud optical properties were estimated by the ice cloud optical parameterization developed by Baran et al. (2014) (the “Baran 2014” parameterization) and the liquid water cloud optical parameterization in terms of effective diameter (D_eff) (the “D_eff” parameterization) (Hocking et al., 2019). Considering that a wide selection of liquid water cloud and ice cloud optical parameterizations were available in RTTOV, the configuration of cloud optical parameterizations in Zhou et al. (2024) exhibited a certain degree of blindness.

The cloud optical parameterizations in RTTOV were developed based on theoretical computations and in situ measurements of cloud microphysical properties such as the particle size distribution (PSD) or the shapes and mixing ratios of different ice habits (for ice clouds only). (Baum et al., 2011). It is possible that there are discrepancies in cloud microphysical properties between the built-in parameterizations and NWP models (or real cases). As a result, the performance of RTTOV in synthesizing visible satellite images should vary with the configuration of different cloud optical parameterizations.

In order to provide guidance for optimizing the configuration of cloud optical parameterizations in RTTOV for the DA of the satellite visible reflectance data, the performance of RTTOV configured with all available liquid water cloud and ice cloud optical parameterizations was evaluated by comparing the observed and the synthetic visible satellite images. The remaining part of this paper is organized as follows. Cloud optical parameterizations are introduced in Sect. 2. Experiment designs, data, and the method are presented in Sect. 3. Results for the reference experiment are shown and discussed in Sect. 4. Discussions on the generalization of the results for different experiment designs and comparison with a previous study are presented in Sect. 5. Conclusions are summarized in Sect. 6.

2 Cloud optical parameterizations in RTTOV

All available cloud optical parameterizations built into RTTOV (version 12.3) (Saunders et al., 2018) are summarized in Table 1. For the liquid water cloud and ice cloud, cloud optical parameterizations were further divided into several subtypes, which are described in Sect. 2.2.1 and 2.2.2, respectively.

Table 1 The liquid water cloud parameterization (clw_scheme) and ice cloud parameterization (ice_scheme) in version 12.3 of the RTTOV software package (Hocking et al., 2019).

2.1 Liquid water cloud optical parameterizations

2.1.1 The OPAC parameterization

The Optical Properties of Aerosols and Clouds (OPAC) parameterization in RTTOV provides optical properties for five liquid water cloud types, comprising the continental stratus (STCO), the maritime cumulus (STMA), the clean continental stratus (CUCC), the polluted continental stratus (CUCP), and the maritime cumulus (CUMA) (Hess et al., 1998). For each of the five cloud types, the PSD was described by a pre-defined modified gamma function, and the D_eff determined by the pre-defined PSD was 8.0, 14.6, 11.5, 22.6, and 25.3 µm for CUCP, STCO, CUCC, STMA, and CUMA, respectively. The volumetric optical properties for unit cloud concentration (1 g m⁻³) were calculated by integrating the optical properties derived from Mie theory over the PSDs. In the RTTOV simulations, the extinction and scattering coefficients of liquid water cloud are computed as the product of the respective per-unit-concentration coefficients and the total liquid hydrometeor concentration (comprising cloud droplets and raindrops for the CMA-MESO forecasts).

To generate synthetic visible satellite images, the atmospheric fields were processed into the format of the RTTOV input files for each of the five liquid cloud types following the flowchart shown by Fig. 1. A pre-trail experiment indicates that the differences in simulated reflectance between CUCC and CUCP were on the order of 10⁻³. Therefore, liquid water cloud over land was set to CUCC.

Figure 1 The flowchart for assigning one of the five OPAC liquid water cloud subtypes from the atmospheric fields.

2.1.2 The “D_eff” parameterization

Cloud optical properties for the “D_eff” parameterization were parameterized in terms of D_eff and the mixing ratios of liquid water cloud hydrometeors. Since the atmospheric fields, introduced in Sect. 3.2.1 and 3.2.2, did not include D_eff, D_eff was estimated by Eq. (1) explicitly (Thompson et al., 2004):

$$\begin{array}{}\text{(1)}& D_{\text{eff}}={\left({\displaystyle \frac{\mathrm{1.91}{\mathit{\rho }}_{\mathrm{a}}{q}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}{N}_{\mathrm{w}}}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where ρ_a and ρ_w denote the density of air and cloud water droplets, respectively. q_w denotes the mixing ratio of cloud droplets, and N_w denotes the number concentration of cloud droplets. ρ_a and q_w were derived from the CMA-MESO forecasts. In this study, ρ_w was set to 1000 kg m⁻³. In addition, N_w was set to 300 cm⁻³, which is the default number concentration of cloud condensation nuclei (CCN) in the CMA-MESO model configured with the single-moment six-class microphysical scheme (the default physics scheme for operational forecasting). The presumed CCN number concentration may influence the overall results, and this potential impact is discussed in Sect. 5.5.

Following the parameterization scheme of Thompson et al. (2004), raindrops were formed from cloud droplets through an autoconversion process involving collision and coalescence. Therefore, Eq. (1) does not explicitly account for the impact of raindrops on D_eff. Nevertheless, the mixing ratios of cloud droplets and raindrops were summed (representing total liquid hydrometeor concentration) to capture the influence of liquid cloud hydrometeors on the radiative transfer simulations.

2.2 Ice cloud optical parameterizations

2.2.1 The “SSEC” parameterization

For the ice cloud optical parameterization developed by the Space Science and Engineering Center (SSEC) (the “SSEC” parameterization), the ice cloud optical properties were parameterized in terms of ice water content (IWC) and D_eff. D_eff was not provided by the RTTOV users explicitly. Instead, D_eff was estimated by the following four built-in parameterizations.

Using the parameterization of Ou and Liou (1995) (OL95 hereafter), D_eff was parameterized in terms of ambient temperature T_c in degrees Celsius:

$$\begin{array}{}\text{(2)}& D_{\text{eff}}=\mathrm{326.3}+\mathrm{12.42}{T}_{\mathrm{c}}+\mathrm{0.197}{T}_{\mathrm{c}}^{\mathrm{2}}+\mathrm{0.0012}{T}_{\mathrm{c}}^{\mathrm{3}}.\end{array}$$

Using the parameterization of Wyser (1998) (W98 hereafter), D_eff was parameterized in terms of ambient temperature T in kelvins and IWC (unit: g m⁻³):

$$\begin{array}{}\text{(3)}& D_{\text{eff}}=\mathrm{377.4}+\mathrm{203.3}B+\mathrm{37.91}{B}^{\mathrm{2}}+\mathrm{2.3696}{B}^{\mathrm{3}}.\end{array}$$

B is described by Eq. (4):

$$\begin{array}{}\text{(4)}& B=-\mathrm{2}+{\mathrm{10}}^{-\mathrm{3}}(\mathrm{271}-T{)}^{\mathrm{1.5}}{\log }_{\mathrm{10}}(\text{IWC}/\mathrm{50}).\end{array}$$

Using the parameterization of Boudala et al. (2002) (B02 hereafter), D_eff was parameterized in terms of T_c and IWC:

$$\begin{array}{}\text{(5)}& D_{\text{eff}}=\mathrm{53.005}{\text{IWC}}^{\mathrm{0.06}}\exp \left(\mathrm{0.0013}{T}_{\mathrm{c}}\right).\end{array}$$

Using the parameterization of McFarquar et al. (2003) (MF03 hereafter), D_eff was parameterized in terms of IWC:

$$\begin{array}{}\text{(6)}& D_{\text{eff}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{1.78449}+\mathrm{0.281301}\log \left(\text{IWC}\right)+\mathrm{0.01777166}\left[\log \right(\text{IWC}){]}^{\mathrm{2}}}.\end{array}$$

Apart from the four built-in parameterizations of D_eff for ice cloud, we evaluated an extra parameterization in terms of the number concentration of ice nuclei (IN) (described by Eqs. 2–5 in Yao et al., 2018). However, the IN-related parameterization did not generate a better performance than the four built-in parameterizations. Therefore, the results for the IN-related parameterization are not shown for simplicity.

Similarly to the liquid water clouds, the extinction and scattering coefficients of ice cloud are computed as the product of the respective per-unit-concentration coefficients and the total ice hydrometeor concentration. The total ice hydrometeor concentration is the sum of the concentration for ice, snow, and graupel for the CMA-MESO forecasts.

2.2.2 The “Baran” parameterization

For the “Baran” parameterization, the ice cloud optical properties were parameterized in terms of IWC and temperature (Vidot et al., 2015). Unlike the “SSEC” parameterization, the “Baran” parameterization does not have dependence on D_eff. The old “Baran” parameterization (“Baran 2014”; see Baran et al., 2014) was updated (“Baran 2018”; Saunders et al., 2020) to smooth the discontinuous variation in absorption and scattering coefficients within the shortwave spectral range (Saunders et al., 2020). Therefore, Baran 2018 is more spectrally consistent than Baran 2014, and as such it is recommended over the latter by Hocking et al. (2019).

In summary, there are 12 combinations of cloud optical parameterizations (including the four built-in parameterizations of D_eff for ice cloud) for liquid water cloud and ice cloud (Table 2). The performance of different combinations of cloud optical parameterizations was evaluated by the experiment designs introduced in Sect. 3.1. Based on the configuration of cloud optical parameterizations, radiative transfer simulations were performed by the DOM solver and the Method for FAst Satellite Image Simulation (MFASIS). For the DOM solver, 16 streams were used. The general radiative transfer options account for atmospheric refraction and curvature. The surface bidirectional reflectance distribution function (BRDF) was either derived from monthly mean atlases for land surface (Vidot and Borbás, 2014; Vidot et al., 2018) or calculated by the JONSWAP (Hasselmann et al., 1973) solar BRDF model for the sea surface. Since aerosol variables were not provided by the atmospheric fields, the contribution of aerosols to visible reflectance was neglected during the radiative transfer simulations. Nevertheless, the aerosol impact on the evaluation results is discussed in Sect. 5.6.

Table 2 Different combinations of cloud optical parameterizations. clw_scheme and ice_scheme denote the liquid water cloud and ice cloud optical parameterizations, respectively. “idg” denotes the built-in parameterization of D_eff for ice cloud. idg = 1: OL95; idg = 2: W98; idg = 3: B02; idg = 4: MF03.

3 Experiment designs, data, and method

3.1 Experiment designs

The performance of different cloud optical parameterizations was evaluated by four experiments summarized in Table 3. The experiments cover different observing systems (the visible bands of FY-4B and Himawari-9, where the latter is the operational geostationary meteorological satellite operated by the Japan Meteorological Agency), different atmospheric fields (the CMA-MESO forecasts and the fifth-generation European Centre for Medium-Range Weather Forecasts reanalysis (ERA5) data), and different radiative solvers (DOM and MFASIS) to generalize the main findings.

Table 3 The four experiments covering different atmospheric fields, observing systems, and radiative solvers.

The comparison was performed for band 2 (centered at 0.65 µm) of AGRI on board FY-4B and band 3 (centering at 0.64 µm) of the Advanced Himawari Imager (AHI) on board Himawari-9. Himawari-9 band 3 was chosen due to the spectral and spatiotemporal matches with the FY-4B visible band. On one hand, the visible bands of the two satellites share considerable similarities with respect to spectral characteristics (Fig. 2). Therefore, results derived from the two instruments should have certain consistency. On the other hand, FY-4B and Himawari-9 are above the Equator at the longitudes of 104.7 and 140.7° E, respectively. In addition, the full-disk scanning cycles of AGRI and AHI are 15 and 10 min, respectively. Therefore, the observation areas and times of the two satellites exhibit substantial overlaps.

Figure 2 The spectral response functions for band 2 of AGRI on board FY-4B and for band 3 of AHI on board Himawari-9.

MFASIS is a lookup-table (LUT)-based emulator of a 1D solver (Scheck et al., 2018). When constructing the LUT for the MFASIS, radiative transfer simulations were also affected by the cloud optical parameterizations. Therefore, synthetic images generated by the DOM and MFASIS solvers should reveal some common characteristics in terms of the cloud optical parameterizations built into the RTTOV software package.

Each of the four experiments includes 12 configurations of cloud optical parameterizations (Table 2), except for the CM-HW-MF experiment, which only involved the “SSEC” ice cloud optical parameterization. In fact, MFASIS was initially designed for the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the satellites of the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT). For version 12.3 of RTTOV, MFASIS was extended to the Himawari-9 observing systems, but only the “SSEC” ice cloud optical parameterization was available.

3.2 Data

3.2.1 The 3 h forecasts of the CMA-MESO model

The CMA-MESO forecasts from 15 March to 25 April 2024 were used to provide the inputs to RTTOV. The CMA-MESO model is a cycled DA and forecasting system. The model was initialized eight times per day at 00:00, 03:00, 06:00, 09:00, 12:00, 15:00, 18:00, and 21:00 UTC. At 00:00 and 12:00 UTC, the model was cold-started, with the initial conditions (ICs) and lateral boundary conditions (LBCs) provided by the CMA Global Forecasting System (GFS). At other startup times, the model was warm-started, with the ICs and LBCs updated from the analysis fields that were generated by a cloud analysis technique and by assimilating synergic observations with a 3D variational (3DVar) DA method (Shen et al., 2020).

To avoid low solar elevation, only the 3 h forecasts at 06:00 UTC were selected. The 3 h forecasts were chosen based on an evaluation study of the CMA-MESO model (originally termed the GRAPES_3 km model) by Zhang et al. (2020). The study suggested that the CMA-MESO forecasts were much more reliable within the initial 3 forecasting hours, and the “spin-up” issue was essentially absent with the incorporation of a cloud analysis technique. The domain coverage of the CMA-MESO model includes 2501 × 1671 horizontal grids with a grid spacing of 0.03° and 50 vertical layers with a model top of 10 hPa. To avoid 3D radiative effects on high-resolution radiative transfer simulations, the CMA-MESO forecasts were superobbed to a horizontal resolution of 0.09° with 833 × 557 horizontal grids. The superobbing was performed by simply averaging the 0.03° × 0.03° products every three grids. In the following, the CMA-MESO forecasts refer to the 0.09° × 0.09° products unless otherwise specified.

3.2.2 The ERA5 data

The 0.25° × 0.25° gridded ERA5 data (Hersbach et al., 2020) on pressure levels (Hersbach et al., 2023a) and on single levels (Hersbach et al., 2023b) were used to provide the inputs to RTTOV. The ERA5 data on pressure levels were generated by the ECMWF Integrated Forecast System (IFS) based on worldwide observations using a 4DVar DA method. The pressure-level data provide temperature, the water vapor mixing ratio, and cloud-related parameters (the mixing ratio of cloud, ice, rain, and snow, as well as cloud cover) on 137 pressure levels. In addition, the ERA5 data at the surface level were generated by a coupled model of the atmospheric model in the IFS and a land-surface model. The surface-level data provide the 2 m (and 10 m) wind, 2 m (and 10 m) temperature, 2 m (and 10 m) humidity, surface pressure, etc.

3.2.3 FY-4B data

The FY-4B 4 km × 4 km visible reflectance data, cloud mask (CLM) product, and synchronous observation geometry (GEO) data were used. To generate spatially collocated observed and synthetic images, the FY-4B full-disk visible reflectance data were horizontally averaged to the locations of the CMA-MESO forecasts (or ERA5) grids. The horizontal averaging was performed by the following two steps: first, centering at a given CMA-MESO (or ERA5) grid and finding all the pixels (matched pixels hereafter) in the FY-4B/AGRI visible image within ± 0.045° for the 0.09° × 0.09° CMA-MESO forecasts (or within ± 0.125° for the 0.25° × 0.25° ERA5 data) in both the zonal and the meridional directions, and second, averaging the reflectances of all these matched pixels to generate a reflectance that is spatially matched to the selected CMA-MESO (or ERA5) grid. Repeating the two steps for all CMA-MESO (or ERA5) grid points generated an observed image gridded at 0.09° × 0.09° (or 0.25° × 0.25°). The full-disk scanning cycle of FY-4B/AGRI is 15 min, and the scanning starts at 00:00 UTC. In addition, the CMA-MESO forecasts and ERA5 data were produced at hourly intervals. Therefore, the maximum allowable time difference between the FY-4B observations and CMA-MESO forecasts and ERA5 data was within 15 min to ensure a temporal match.

To facilitate the radiative transfer simulations by RTTOV, the sun-viewing geometries (i.e., solar zenith and azimuth angles, satellite zenith and azimuth angles) were derived from FY-4B GEO data. The 4km × 4 km GEO data were horizontally averaged to the CMA-MESO (or ERA5) grids in the same way as described above. In addition, the FY-4B CLM product was used to provide a first-step estimate of cloud mask for an arbitrary CMA-MESO or ERA5 grid. The 4km × 4 km CLM product was matched to the CMA-MESO (or ERA5) grids by a maximum occurrence method; a i.e., cloud mask for a CMA-MESO (or ERA5) grid was set to the maximum occurrence frequency amongst the entire matched FY-4B pixels.

3.2.4 Himawari-9 data

The Himawari-9 5 km × 5 km products (Bessho et al., 2016) – including the reflectance at band 3, the GEO data, and cloud type product – were matched to the CMA-MESO and ERA5 grids in the same way as introduced in Sect. 3.2.3. The Himawari-9 cloud type product provided information not only on cloud or clear sky for a certain pixel, but also on the cloud types for cloudy scenarios. Clouds were divided into eight subtypes, comprising cirrus (Ci), cirrostratus (Cs), deep convection (Dc), altocumulus (Ac), altostratus (As), nimbostratus (Ns), cumulus (Cu), and stratus (Sc).

3.3 Equivalent criteria of cloud mask for the observed and synthetic images

Since statistical characteristics of the O−B departure are different for cloudy and clear pixels, it is critical to evaluate the results for the two scenarios separately. To ensure equivalent criteria of cloud mask for the observed and synthetic images, the observed and synthetic visible images were compared with the images simulated by ignoring cloud impacts (Zhou et al., 2024). For synthetic images, a pixel was designated cloudy if Eq. (7) was satisfied. Otherwise, the pixel was designated cloud-free.

$$\begin{array}{}\text{(7)}& r_{\text{sim}}>r_{\text{clr}},\end{array}$$

where r_sim denotes the reflectance for an arbitrary pixel in a synthetic image and r_clr denotes the spatiotemporally collocated reflectance simulated by ignoring cloud impacts.

The aerosol contributions were neglected by the RTTOV simulations. However, the observed reflectance inevitably included aerosol contributions. To account for the aerosol impacts on the cloud masking for the observed images, a pixel was designated cloudy if the observed reflectance r_obs satisfied Eq. (8):

$$\begin{array}{}\text{(8)}& r_{\text{obs}}>r_{\text{sim}}+r_{\text{aer}}^{\mathrm{75}},\end{array}$$

where $r_{\text{aer}}^{\mathrm{75}}$ denotes the aerosol contribution to the reflectance for cloudy pixels. $r_{\text{aer}}^{\mathrm{75}}$ was set to the upper quartile of r_obs−r_clr for the preliminarily estimated cloud-free pixels, which were designated by the cloud mask derived from the FY-4B CLM product. The second-step estimate of cloud-free pixels was determined as follows:

$$\begin{array}{}\text{(9)}& r_{\text{obs}}<r_{\text{clr}}+r_{\text{aer}}^{\mathrm{25}},\end{array}$$

where $r_{\text{aer}}^{\mathrm{25}}$ denotes the estimate of aerosol contributions to the cloud-free reflectance. Similarly, $r_{\text{aer}}^{\mathrm{25}}$ was set to the lower quartile of r_obs−r_clr for the preliminarily estimated cloud-free pixels.

Utilizing $r_{\text{aer}}^{\mathrm{75}}$ and $r_{\text{aer}}^{\mathrm{25}}$ to identify aerosol impact in cloudy and cloud-free conditions may contain inherent logical flaws. For instance, inaccurate surface albedo introduces confounding effects into reflectance comparisons. Even under idealized conditions where atmospheric parameters and other RTTOV components are precisely known, discrepancies between simulated and observed reflectances simultaneously reflect the combined influences of both aerosol effects and surface albedo. The aerosol impact on the evaluation of different cloud optical parameterizations (COPs) is further discussed in Sect. 5.6.

Figure 3 Time series of the O−B biases and standard deviations for (a, b) all, (c, d) cloudy, and (e, f) clear pixels. The results are for the CM-FY-DM experiment.

4 Results for the CM-FY-DM experiment

4.1 Temporal variation in the bias and standard deviation

The time series for the biases and standard deviations of the O−B departure from 15 March to 15 April 2024 is shown in Fig. 3. On average for the domain, B was underestimated compared with O, with the O−B biases ranging from 0.05 to 0.12. The O−B biases of cloudy conditions were several orders of magnitude larger than those of the clear conditions. In clear conditions, the unresolved aerosol processes, the errors due to measurement calibration processes, and the errors in BRDF over snow-covered areas or sun-glint areas could contribute to the O−B biases (Zhou et al., 2024). In cloudy conditions, the contributing factors to the errors in B include not only the above-mentioned error sources in clear conditions, but also the errors in the atmospheric fields and the deficiencies of cloud optical parameterizations. The minimum O−B bias was revealed for the C213 cloud optical parameterization, i.e., the “D_eff” liquid water cloud parameterization and the “SSEC” ice cloud parameterization with the effective diameter of ice crystals parameterized by the B02 parameterization.

The maximum difference in the standard deviations for different cloud optical parameterizations was within 0.03. The probability density distribution function (PDF) of reflectance for the synthetic images showed over- or underestimation (compared with the observed images) in the occurrence frequency at low (< 0.03) and moderate-to-large (> 0.2) reflectance levels (Fig. 4). In comparison, the PDF for C112, which was taken as an example for the sub-optimal cloud optical parameterization, is also shown in Fig. 4. The PDF for C213 resembles that for the observation better than other cloud optical parameterizations (not shown for simplicity). Since the PDFs were not sensitive to the location of the clouds (Geiss et al., 2021), the PDF analysis suggested an overall improvement compared with other cloud optical parameterizations.

Figure 4 The probability density distribution functions (PDFs) of reflectance for the observed and synthetic images from 15 March to 15 April 2024.

Figure 5 Dependence of FY-4B/AGRI band 2 reflectance on cloud water path (CWP) for different (a) liquid water cloud optical parameterizations and (b) ice cloud optical parameterizations. For this simulation, the solar zenith angle, viewing zenith angle, and relative azimuth angle are set to 25, 40, and 135°, respectively.

To illustrate the reasons to the performance of different cloud optical parameterizations, the sensitivity study of reflectance to different cloud optical parameterizations was demonstrated by Fig. 5. According to Mie theory, the larger the cloud particle size, the stronger the forward scattering, and vice versa. Therefore, the largest (smallest) reflectance was expected for the CUCP (CUMA) cloud subtype. The “D_eff” parameterization generated larger top-of-atmosphere (TOA) reflectance than the “OPAC” scheme, which implies that the backscattering effects of the “D_eff” parameterization are stronger than the “OPAC” parameterization.

Figure 6 (a) The effective diameter estimated by the four built-in parameterizations in RTTOV. (b) The vertical distribution of temperature and ice water content for the sensitivity study shown by Fig. 5.

For the ice cloud optical parameterizations, Baran 2014 outperformed Baran 2018 when measured by the O−B biases. The mean volumetric scattering coefficients of Baran 2014 were larger than those of Baran 2018 (Saunders et al., 2020). Therefore, more photons were backscattered for the Baran 2014 parameterization. The results in Figs. 3 and 5 suggest that the Baran 2018 ice scheme should be used with caution for the radiative transfer simulations in visible spectral ranges. In addition, the performance of the “SSEC” ice cloud optical parameterization was sensitive to the parameterization of D_eff for ice clouds. An inter-comparison between the four parameterizations of D_eff is illustrated by Fig. 6, which reveals that the minimum (maximum) effective diameter was found for the B02 (W98) parameterization. Therefore, the largest (smallest) reflectance was expected for the B02 (W98) parameterization.

The optical characteristics of the “Baran 2014 + B02” ice cloud parameterization may partially explain the observed phenomenon where the C213 parameterization successfully corrects the systematic underestimation of B while introducing greater standard deviation (or local variability) (Fig. 3a–d). Within the same cloud water path (CWP) range, C213 produces the largest variation amplitude in B. Notably, when transitioning from cloud-free to cloudy conditions, the “Baran 2014 + B02” scheme exhibits the most significant reflectance variation range (Fig. 5b). This enhanced variability in B could potentially broaden the distribution of O−B discrepancies, thereby increasing the standard deviation of O−B. However, it is important to note that the standard deviations of the O−B departure showed less sensitivity to cloud optical parameterizations compared to the biases.

Figure 7 (a) Spatial distribution of the 1-month bias of O−B departure for the CM-FY-DM experiment; (b) spatial distribution of the 1-month standard deviation of O−B departure; (c) spatial distribution of the 1-month correlation coefficient between O and B.

Figure 8 (a) Spatial distribution of the O−B biases for collocated cloud-free pixels from 15 March to 15 April 2024; (b) spatial distribution of the mean BRDF from 15 March to 15 April 2024.

4.2 Spatial distribution of the O−B departure

One-month O−B departure statistics over the study domain reveal systematic errors and aid in understanding of results. The spatial distribution of the O−B biases, the standard deviations of the O−B departure, and the correlation coefficients between O and B were derived from the 1-month observed and synthetic images (Fig. 7). Compared with O, B was underestimated over Siberia, the Mongolian Plateau, the southern foothills of the Himalayas, the Sichuan basin, and the Yunnan–Kweichow Plateau (Fig. 7a). In general, the spatial distribution of the correlation coefficients between O and B agreed well with the spatial distribution of the O−B biases (Fig. 7c). Namely, small O−B biases agreed well with large correlation coefficients. Siberia, the Mongolian Plateau, and the southern foothills of the Himalayas were covered with snow in March and April. Reflectance simulated in these areas should be less accurate compared with other places because the BRDF atlas is questionable in snow-covered areas (Ji et al., 2022). Since the areas with O−B bias larger than zero were in good agreement with the areas potentially being covered by snow (Fig. 8), a tentative conclusion could be drawn that the surface albedo (equivalent to the BRDF for Lambertian radiators) over the snow-covered areas was underestimated.

Figure 9 (a) Visible image observed by FY-4B/AGRI band 2 at 06:00 UTC on 28 March 2024. (b) Synthetic visible image of FY-4B/AGRI band 2 simulated by RTTOV (C213 configuration). (c) Difference between the observed and synthetic images (observed minus synthetic). (d) Spatiotemporally collocated cloud types derived from the Himawari-9 cloud type product.

In addition, the performance of the CMA-MESO model was reduced in some circumstances. An example for the observed and synthetic images at 06:00 UTC on 28 March 2024 is shown by Fig. 9. The selected case reported a typical comma-shaped cloud, which was caused by a frontal system. In general, B was undervalued in some parts of Eastern China and the Western Pacific areas (Fig. 9c), which were classified as featuring convective clouds (Fig. 9d). Similar results were reported by Zhou et al. (2024). A potential explanation for this is the deficiency of the CMA-MESO model in forecasting the strongly convective weather systems (Wan et al., 2015). In addition, some of the Ci and Cs clouds were missed or underestimated over the Himalayas, the Sichuan basin, and the Yunnan–Kweichow Plateau. This is most likely caused by the reduced performance of the CMA-MESO model over complex-terrain areas (Wan et al., 2015; Zhu et al., 2017).

Figure 10 Time series of the O−B biases and standard deviations for (a, b) all, (c, d) cloudy, and (e, f) clear pixels. The results are for the E5-FY-DM experiment.

5 Discussions

5.1 Results for E5-FY-DM: influences of atmospheric fields

The comparison between CM-FY-DM and E5-FY-DM experiments was designed to reveal the influences of different atmospheric fields on the statistical characteristics of O−B departure. The results are shown by Fig. 10.

In general, the results for E5-FY-DM are similar to those for CM-FY-DM. To be specific, the minimum O−B bias was revealed for the cloud optical parameterization of C213, and the standard deviations were less sensitive to cloud optical parameterizations compared with the O−B biases. However, the biases and standard deviations of the O−B departure were smaller than those for the CM-FY-DM experiment (Fig. 3). There are two potential explanations. On one hand, the coarser grids for the E5-FY-DM experiment meant more averaging of the atmosphere fields and the subsequent reflectance fields. The horizontal averaging tended to smooth the reflectance fields; i.e., the larger reflectance tended to be reduced and vice versa. As a result, the PDF of the O−B departure shrank, and the differences between O and B were reduced. On the other hand, it was possible that atmospheric fields, especially cloud variables, were better represented by the ERA5 data than by the CMA-MESO forecasts. Since errors in B were mainly determined by atmospheric fields and the forward operators, consistent results for different atmospheric fields increased the robustness of the findings in Sect. 4.

Figure 11 Time series of the O−B biases and standard deviations for (a, b) all, (c, d) cloudy, and (e, f) clear pixels. The results are for the CM-HW-DM experiment.

5.2 Results for CM-HW-DM: influences of observing systems

The comparison between CM-FY-DM and CM-HW-DM experiments was designed to reveal the influences of different observing systems on the statistical characteristics of the O−B departure. The results for the CM-HW-DM experiment are shown by Fig. 11.

In general, the results for the CM-HW-DM experiment were similar to those for CM-FY-DM. However, the O−B biases and standard deviations were smaller for the CM-HW-DM experiment than for CM-FY-DM. Although the centering wavelengths of the FY-4B and Himawari-9 visible bands are close, the spectral response function has a wider range for the FY-4B visible band than for the Himawari-9 visible band (Fig. 2). As a result, the convolution of monochromatic cloud optical properties over the spectral response function would involve a broader wavelength range. If cloud optical parameterizations contain errors, it is likely that including a broader wavelength range would amplify the errors in volumetric optical properties, leading to larger O−B biases and standard deviations. In addition, the differences in the O−B departure between the two observing systems could be related to the measurement calibration processes. The radiometric calibration techniques for the two visible bands were performed by different methods (Okuyama et al., 2018; Zhang et al., 2025). Since 29 May 2023, the National Satellite Meteorological Center (NSMC) has not updated the calibration coefficients of the FY-4B/AGRI solar reflection bands. It is possible that the radiometric performance of the instrument declined due to the influences of space particle erosion, device aging, etc.

Figure 12 Time series of the O−B biases and standard deviations for (a, b) all, (c, d) cloudy, and (e, f) clear pixels. The results are for the CM-HW-MF experiment.

5.3 Results for CM-HW-MF: influences of radiative solvers

The comparison between CM-HW-MF and CM-HW-DM experiments was designed to reveal the influences of different radiative solvers on the statistical characteristics of O−B departure. The results for the CM-HW-MF experiment are shown by Fig. 12.

In general, similar results were revealed between the two experiments. However, the O−B biases and standard deviations were smaller for CM-HW-MF (Fig. 12) than for CM-HW-DM. Since B was smaller than O on average for the domain, the reflectance simulated by MFASIS was slightly larger than that simulated by the DOM solver. This is most likely related to the differences between the two solvers in tackling the scattering interactions between clouds and gaseous molecules. The multiple-scattering processes in MFASIS were considered for all cloudy and clear layers (Scheck et al., 2016). However for the DOM solver, the multiple-scattering processes were only limited to cloudy layers, and the scattering interactions between clouds and gaseous molecules were simplified into single-scattering processes. As a result, the MFASIS-simulated reflectance was larger than the DOM-simulated reflectance, especially in dense cloud areas where the scattering interactions between clouds and gaseous molecules were non-negligible (Scheck et al., 2016).

Figure 13 (a) Spatial distribution of the mean BRDF in March 2024; (b) spatial distribution of the mean BRDF in April 2024.

5.4 Comparison with a previous study on FY-4A observing systems

The O−B biases were also explored for FY-4A/AGRI channel 2 based on the CMA-MESO forecasts and RTTOV-DOM forward operator in a previous study by Zhou et al. (2024). In Zhou et al. (2024), the “D_eff” liquid cloud optical parameterization and the “Baran 2014” ice cloud optical parameterization were used. It is obvious that the configuration of the cloud optical parameterization in Zhou et al. (2024) was sub-optimal. Nevertheless, the results in this study revealed many common characteristics with Zhou et al. (2024). For example, B consistently underestimated O on average for the domain. In addition, an abrupt change was reported from 8 to 9 September 2020, shown by Fig. 7 in Zhou et al. (2024), which was caused by the update of the calibration coefficients. In this study, an abrupt change in the O−B biases was also revealed from 31 March to 1 April for the four experiments (Figs. 3a, 10a, 11a, and 12a). The abrupt changes in this study were more likely related to the abrupt change in the BRDF from March to April (Fig. 13). In RTTOV, the BRDF of the underlying land surface was taken from a monthly mean atlas, meaning it was quasi-static and could not reflect the true temporal-variation characteristics.

One topic of Zhou et al. (2024) was to promote a bias-correction method based on the first-order approximation of the O−B bias. The findings in this study provide guidance for estimating more accurate B and the subsequent bias-correction coefficient. Therefore, this study can be regarded as an extension of Zhou et al. (2024).

Figure 14 Sensitivity of the reflectance, denoted by R, to the pre-assumed CCN number concentration, denoted by N. (a) Variation in reflectance with the cloud water path (CWP) for different N values. (b) Variation in reflectance difference with CWP for different N values. The simulation was conducted for FY-4B band 2 by RTTOV-DOM, with the satellite zenith angle, satellite azimuth angle, solar zenith angle, and solar azimuth angle set to 36, 0, 50, and 55°, respectively.

Figure 15 Time series of the O−B biases and standard deviations for (a, b) all, (c, d) cloudy, and (e, f) clear pixels. The settings of the experiment are similar to those of the CM-FY-DM experiment (CCN number concentration = 300 cm⁻³), except that the CCN number concentration was set to 500 cm⁻³ to calculate the effective diameter for liquid water cloud.

5.5 Influences of pre-assumed CCN number concentration

The pre-assumed CCN number concentration in Eq. (1) impacts D_eff for liquid water cloud. An increase in CCN number concentration distributes available water vapor among a greater number of aerosol particles, resulting in the formation of smaller cloud hydrometeors. Conversely, a reduction in CCN leads to fewer but larger hydrometeors and vice versa. The typical CCN number concentration varies from 50, 100, and 200 to 500 cm⁻³ (Thompson et al., 2004). Accordingly, we conducted two additional experiments by configuring the CCN number concentration in Eq. (1) as 100 and 500 cm⁻³. The sensitivity of TOA reflectance to the CCN number concentration is shown by Fig. 14. For thin to moderately thick clouds (CWP < 3.2 (i.e., 10^0.5) kg m⁻², Fig. 14), higher CCN concentrations yield reduced TOA reflectance, as smaller cloud hydrometeors enhance backscattering. In contrast, for optically thick clouds (CWP > 3.2 kg m⁻²), increased CCN number concentrations result in greater TOA reflectance, likely due to more efficient multiple-scattering processes. Nevertheless, the pre-assumed CCN number concentration demonstrates minimal impact on the evaluation results. The simulation results obtained with a CCN number concentration of 500 cm⁻³ (Fig. 15) show remarkable consistency with those generated using 300 cm⁻³ (Fig. 3). The results suggest that CCN's influence becomes less pronounced in synthetic visible satellite images. This occurs primarily because liquid water clouds are often overlain by upper-level ice clouds. Since visible-wavelength radiation has limited penetration depth through cloud layers, the observed reflectance is predominantly sensitive to the properties of the upper ice cloud layer rather than the underlying liquid water clouds affected by CCN. Consistent results were observed for the 100 cm⁻³ CCN case (omitted for brevity), further demonstrating the robustness of our findings to variations in the prescribed CCN concentration.

5.6 Influences of aerosol impact derived from the FY-4B aerosol products

In Eqs. (8) and (9), the aerosol impact was crudely identified by $r_{\text{aer}}^{\mathrm{75}}$ and $r_{\text{aer}}^{\mathrm{25}}$ for cloudy and cloud-free scenarios, respectively. In this part, the aerosol impact is further identified based on the FY-4B aerosol product. We conducted an extra experiment to account for the aerosol impact on the evaluation results. The experiment settings were similar to the CM-FY-DM experiment, except that the aerosol optical properties were included in the RTTOV inputs when synthesizing visible satellite images. For a certain atmosphere column, the aerosol-related inputs to RTTOV include the extinction coefficient profile, the scattering coefficient profile, and the phase function at specified scattering angles. The aerosol extinction coefficient profile was derived from FY-4B aerosol optical depth (AOD) products at 650 nm (data quality flag = 3), assuming an exponential decrease in aerosol concentration with a scale height of 3 km. The assumption was based on the aerosol climatology over China (Zhou et al., 2017), which revealed that the vertical distribution of aerosols commonly conforms to an exponential function and the typical value of the scale height is around 3 km. In addition, sand dust is the most predominant aerosol type. Therefore, we used the phase function of the sand dust aerosol type. The per-layer scattering coefficient was calculated by multiplying the extinction coefficient by the single-scattering albedo of sand dust aerosol. The optical properties of the dust aerosol at the central wavelength of FY-4B AGRI band 2 were calculated by a logarithmic interpolation of the optical properties at 0.532 and 1.064 µm provided by Zhou et al. (2017).

Figure 16 Time series of the O−B biases and standard deviations for (a, b) all, (c, d) cloudy, and (e, f) clear pixels. The settings of the experiment are similar to those of the CM-FY-DM experiment (aerosols are ignored in RTTOV simulations), except that the aerosol impact was included in the RTTOV simulations to synthesize visible satellite images.

The FY-4B AOD products only provide valid retrievals under clear-sky conditions. Due to the cloud displacement errors in the CMA-MESO forecasts, the FY-4B AOD product provides aerosol optical properties over a pixel which is either cloudy (r_sim>r_clr) or cloud-free for the CMA-MESO forecasts. Therefore, $r_{\text{aer}}^{\mathrm{75}}$ in Eq. (8) was replaced by the mean differences in reflectance simulated with and without aerosol impact over cloudy pixels. Similarly, $r_{\text{aer}}^{\mathrm{25}}$ in Eq. (9) was replaced by the mean differences in reflectance simulated with and without aerosol impact over cloud-free pixels. Based on the experiment design, we re-calculated the biases and standard deviations of the O−B departure. The results indicate general consistency in cloudy conditions with the results for the CM-FY-DM experiment (Figs. 3 and 16). This aligns with established understanding, as aerosol influences become negligible under cloudy conditions, particularly with optically thick clouds. However, O was generally smaller than B in cloud-free conditions (Fig. 16e), which is contrary to the results shown by Fig. 3e. Including aerosols in the RTTOV inputs for cloud-free conditions tends to increase the reflectance. Although the systematic underestimation of O relative to B in Fig. 16e suggests a potential overestimation of the aerosol impact in cloud-free conditions, the optimized cloud optical parameterizations derived from Fig. 16 remain consistent with those shown in Figs. 3, 10–12, and 15.

6 Conclusions

The performance of RTTOV, a commonly used forward-operator software package, is critical to the DA of satellite visible reflectance data in many aspects. During the radiative transfer simulations of RTTOV, cloud optical properties were determined by several built-in parameterizations in terms of the cloud water content, cloud effective diameter, and ambient temperature. The radiative transfer is influenced by the cloud optical parameterizations. However, it is unclear which combination of liquid water cloud and ice cloud optical parameterizations reproduces the observed reflectance best.

In view of this problem, the performance of RTTOV under different cloud optical parameterizations was evaluated based on observed and synthetic visible satellite images. To generalize the main findings, four experiments were performed. The experiments covered two observing systems, the FY-4B and Himawari-9 visible bands; two radiative solvers, DOM and MFASIS; and two atmospheric fields, the CMA-MESO forecasts and ERA5 data. Statistical characteristics of the O−B departure varied with the observing systems, the representativeness of the atmospheric fields, and the accuracy of radiative transfer modeling (e.g., scattering interactions between clouds and gaseous molecules). Nevertheless, consistent findings were revealed for different experiment designs.

In general, the O−B biases were sensitive to the cloud optical parameterizations. An analysis of 1-month O−B biases revealed that the simulated reflectance was lower than the observed reflectance. The smallest O−B bias was revealed for the liquid water cloud optical parameterization in terms of effective diameter (the “D_eff” parameterization). This was in combination with the ice cloud optical parameterization developed by SSEC (the “SSEC” parameterization). For the “SSEC” parameterization, the effective diameter for ice clouds was parameterized in terms of ambient temperature and IWC, i.e., the B02 parameterization. Compared with the worst-performing cloud optical parameterization, the optimized configuration achieved a bias reduction of 0.02–0.04 on average, with only a marginal increase (< 0.01) in standard deviation. Despite the largest standard deviation of the O−B departure being revealed for the “D_eff + SSEC + B02” parameterization, the standard deviations were less sensitive to the cloud optical parameterizations. In addition, the PDF of the reflectance for the synthetic images simulated by the “D_eff + SSEC + B02” parameterization resembled the best one for the observed images. Therefore, the optimal cloud optical parameterization was suggested to be the “D_eff + SSEC + B02” parameterization suite.

It is noted that a skewed PDF for the O−B departure would be expected even for the optimal cloud optical parameterization. Since conventional DA methods assume that the PDF of observation errors conforms to an unbiased Gaussian function, it is necessary to correct the systematic biases of the FY-4B or Himawari-9 visible reflectance data for DA applications. A potential bias-correction method was promoted based on the first-order approximation of O−B biases (Zhou et al., 2024). Since O is statistically larger than B, the optimal cloud optical parameterization would generate the smallest bias-correction coefficient (defined as $\stackrel{\mathrm{‾}}{O-B}/\stackrel{\mathrm{‾}}{O}$, where the upper horizontal lines denote the domain averaging). As a result, fewer extra errors in the background fields would be introduced into the observations for the optical cloud optical parameterization.

In addition, the “D_eff” liquid water cloud parameterization explicitly depends on the effective diameter of liquid water cloud, which is not a state variable for most of the NWP models. The effective diameter was usually parameterized in terms of the mixing ratio of cloud droplets, the ambient temperature, etc. For the DA of satellite visible reflectance data in real-world cases, the parameterization of the effective diameter for liquid clouds should be assessed. In addition, the performance of RTTOV with the optimal cloud parameterization should be tested by DA experiments in real-world cases and the results should be evaluated by synergic observations. Extending the optimal configuration of cloud optical parameterization in RTTOV to the bias correction and DA applications of FY-4B and Himawari-9 visible reflectance data is ongoing.