A preliminary evaluation of FY-4A visible radiance data assimilation by the WRF (ARW v4.1.1)/DART (Manhattan release v9.8.0)-RTTOV (v12.3) system for a tropical storm case

Satellite visible radiance data that contain rich cloud and precipitation information are increasingly assimilated to improve the forecasts of numerical weather prediction models. This study evaluates the Data Assimilation Research Testbed (DART, Manhattan release v9.8.0), coupled with the Weather Research and Forecasting (WRF) model (ARW v4.1.1) and the Radiative Transfer for TOVS (RTTOV, v12.3) package, for assimilating the simulated visible imagery of the FY-4A geostationary satellite located over Asia in an Observing System Simulation Experiment (OSSE) framework. The OSSE was performed for the tropical storm Higos that occurred in 2020 and contains multi-layer mixed-phase cloud and precipitation processes. The advantages and limitations of DART for assimilating FY-4A visible imagery were evaluated. Both single-observation experiments and cycled data assimilation (DA) experiments were performed to study the impact of different filter algorithms available in DART, variables being cycled, observation outlier thresholds, observation errors, and observation thinning.

The results show that assimilating visible radiance data significantly improves the analysis of the cloud water path (CWP) and cloud coverage (CFC) from first-guess forecasts. The rank histogram filter (RHF) allows WRF to more accurately simulate CWP and CFC compared with the ensemble adjustment Kalman filter (EAKF) although it took roughly twice as long as the latter. By cycling both cloud and non-cloud variables, specifying large outlier threshold values, or setting smaller observation errors without thinning of observations, WRF achieved a better simulation of CWP and CFC. With model integration, DA of the visible radiance data also generated a slightly positive impact on non-cloud variables as they were adjusted through the model dynamics and physics related to cloud processes. In addition, the DA improved the representation of precipitation. However, the impact on the rain rate is limited by the inabilities of the DA to improve cloud vertical structures and cloud phases. A negative impact of the DA on cloud variables was found due to the nature of the non-linear forward operator and the non-Gaussian distribution of the prior. Future works should explore faster and more accurate forward operators suitable for assimilating FY-4A visible imagery, techniques to reduce the non-linear and non-Gaussian errors, and methods to correct the location errors which correspond to the clouds underestimated by the first guess.

1 Introduction

All-sky satellite data assimilation (DA) has shown great potential to improve weather forecasts (Bauer et al., 2011). Many satellite DA-related studies were conducted for microwave (MW) and infrared (IR) radiance data. DA of MW radiance data adjusts the atmospheric state variables such as humidity and temperature (Geer et al., 2019; Migliorini and Candy, 2019) as well as cloud-related parameters such as liquid and ice water content and cloud coverage (Zhang et al., 2013; Yang et al., 2016), exhibiting positive effects on cloud and precipitation forecasting. All-sky MW data have been operationally assimilated at some numerical weather prediction (NWP) centers (Bauer et al., 2010; Zhu et al., 2016). However, operational DA of all-sky MW data is limited to humidity- and temperature-sounding channels (Carminati and Migliorini, 2021) because MW radiance at these channels is insensitive to surface emissivity and skin temperature, which are difficult to be estimated accurately under cloudy-sky conditions (J. Hu et al., 2021​​​​​​​). DA of MW data is also challenging in terms of separating the radiance contribution of clouds from the non-cloud variables (especially temperature and humidity) (Geer et al., 2017). In addition, several studies showed positive effects on water vapor and temperature by assimilating the IR data in clear-sky conditions (McCarty et al., 2009; Ma et al., 2017). DA of IR radiance data in cloudy regions also improved the analysis of column-integrated water and forecasting skills in the mid- and upper troposphere (Stengel et al., 2013; Geer et al., 2019). However, DA of IR radiance data in cloudy regions is still complicated by the non-linear relationship between the observation and state variables, the non-Gaussian problems (Li et al., 2022), the difficulty to separate cloud signals and non-cloud signals (Geer et al., 2017), and the inability to constrain the layered structures in the case of multi-layer clouds (Prates et al., 2014).

Several studies suggested that there is great potential in assimilating the visible (VIS) and shortwave IR (collectively referred to as “shortwave”, SW) data (Vukicevic et al., 2004; Polkinghorne and Vukicevic, 2011; Scheck et al., 2020; Schröttle et al., 2020) because these measurements contain unique and supplementary cloud information to the IR and MW radiance data (Kostka et al., 2014; Schröttle et al., 2020). For example, SW radiation can penetrate a certain depth of cloud fields and connotate cloud microphysical properties such as effective particle radius (Nakajima and King, 1990). In comparison, satellite IR data only provide information on the cloud-top microphysics (Xue, 2009). SW data also complement precipitation radars that mainly measure large hydrometeors or precipitation particles (Keat et al., 2019) that usually do not occur during the initial growth stage of convective systems (Zhang and Fu, 2018), measuring mainly small cloud droplets. Furthermore, SW data usually have higher spatial resolution than MW data (Yang et al., 2017). Therefore, high-resolution satellite SW radiance data provide cloud properties that are of great significance for cloud-resolving models. Unlike MW and IR data, the radiance contributed by clouds could be easily extracted from the VIS observations because the VIS radiance data are much more sensitive to cloud variables than to non-cloud variables.

Several studies have attempted to assimilate the SW radiance data directly (i.e., direct DA). Unlike indirect DA, which assimilates the retrieved cloud parameters from SW radiance data, direct DA critically depends on observation operators. Several observation operators and relevant algorithms have been developed for the satellite VIS radiance DA. For example, Vukicevic et al. (2004) mapped the model state variables to the equivalent radiance by an observation operator for the VIS and IR radiance measurements (VISIROO). Later, Evans (2007) developed the spherical harmonics discrete ordinate method for plane-parallel data assimilation (SHDOMPPDA), which solves radiative transfer processes by discrete ordinate method (DOM) in Cartesian space while computing source functions using spherical harmonic series in spherical space. Compared with observation operators which solve source functions in Cartesian space, SHDOMPPDA has an advantage of high computation efficiency. To further speed up the computation, Scheck et al. (2016a) developed a method for fast satellite image synthesis (MFASIS), which is 2–4​​​​​​​ orders of magnitude faster than the other DOM-based observation operators (Scheck et al., 2016b). Scheck et al. (2018) further improved MFASIS and reduced its errors caused by three-dimensional (3D) radiative effects. MFASIS is one of the observation operators in the radiative transfer for TOVS (RTTOV). RTTOV contains several observation operators for satellite radiance DA (Saunders et al., 2018), including DOM and the single-scattering method for SW radiative processes. These solvers can tackle cloud fraction, parallax correction, and many other critical aspects with respect to molecular absorption and scattering, underlying surface reflection, etc. (Saunders et al., 2018). Apart from these aforementioned observation operators, several machine learning-based observation operators and relevant methods (Scheck, 2021; Zhou et al., 2021) were developed and achieved high computation efficiency and accuracy for the VIS radiance simulations.

Another critical aspect of assimilating satellite VIS radiance data is in the selection of DA algorithms. There are two groups of widely used DA approaches. The first is based on variational (VAR) methods. Vukicevic et al. (2004) assimilated GOES-9 VIS radiance data to the Regional Atmospheric Modeling System (RAMS) with a four-dimensional VAR (4DVAR) method and exhibited positive effects on the short-term forecasting of a stratus cloud field. Similarly, Polkinghorne and Vukicevic (2011) assimilated the GOES-8 VIS and IR radiance data to RAMS with the 4DVAR system and also achieved positive results. The second system is based on ensemble-based methods, which are remarkably stable for non-linear systems and were used for cloud and precipitation studies by several researchers (Lei et al., 2015; Kurzrock et al., 2019). Schröttle et al. (2020) assimilated VIS and IR radiance data in an idealized Observing System Simulation Experiment (OSSE) framework using a local ensemble transform Kalman filter (LETKF). Their results indicated that assimilating VIS radiance data alone could improve the forecasting skills of the regional model, COnsortium for Small-scale MOdeling (COSMO), and assimilating the VIS and IR radiance data collaboratively could further improve the forecasting skills. Their findings were further validated by Scheck et al. (2020), who concluded that assimilating the VIS radiance data of Spinning Enhanced Visible and Infrared Imager (SEVIRI) on METEOSAT could improve cloud and precipitation forecasts and, moreover, reduce the temperature and relative humidity forecast errors in most conditions.

The VAR and ensemble-based approaches are complementary to each other. The ensemble-based approaches generate the flow-dependent background error covariance matrices which can be used to leverage the VAR approaches. Therefore, several hybrid approaches have been developed. Buehner et al. (2013) evaluated an ensemble-VAR DA approach by assimilating the observations that were operationally assimilated in Environment Canada and found it more skillful than the VAR method for the short- and mid-range forecasts over tropical and extra-tropical regions. Gao et al. (2013) developed a hybrid ensemble Kalman filter (EnKF)-3DVAR method to assimilate radar data and found that it outperforms 3DVAR or EnKF in shortening the spin-up time of a supercell storm. In addition, hybrid methods are increasingly applied in the DA of satellite radiance data. Xu et al. (2016) assimilated the FY-3B satellite MW radiance data with the WRF hybrid ensemble/3DVAR and improved the forecasts of typhoon tracks, intensities, and precipitation from 3DVAR. Similar results were also reported by Shen et al. (2020).

Today, there are many community DA resources of ensemble-based methods, such as the Data Assimilation Research Testbed (DART, Anderson et al., 2009). DART supports several numerical weather prediction (NWP) models including the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008). Recently, WRF/DART incorporated the RTTOV radiative transfer package, facilitating the DA of satellite VIS to MW wavelength radiance and enabling the DA of all-sky satellite SW radiance data. The Advanced Geostationary Radiation Imager (AGRI) on the geostationary FY-4A satellite launched in 2016, located over Asia, excels at high sampling frequency (5 min for intensive observation and 15 min for operational observation) and high spatial resolution (0.5–2 km, depending on channels). Zhang et al. (2019) reported that AGRI provides rich short-wave radiance measurements describing rapidly evolving and small- to meso-scale atmospheric systems. However, these data have not been assimilated in the operational NWP centers.

In this study, FY-4A/AGRI VIS radiance data were simulated in an OSSE framework and experiments of the VIS radiance DA were conducted using the WRF/DART-RTTOV system. We intended to answer the following three questions: (1) What are the advantages and limitations of assimilating the FY-4A VIS radiance for forecasting tropical storms? (2) How to specify the WRF/DART-RTTOV model settings and observations? (3) What needs to be considered for real DA applications of FY-4A VIS radiance data? The results of this study are preliminary and contribute toward evaluating the WRF/DART-RTTOV system in assimilating all-sky FY-4A/AGRI VIS radiance data. They should also be applicable to the DA of the upcoming FY-4B VIS radiance data because the designs of the AGRI payload for the FY-4A and -4B are similar. The paper is organized as follows. The models and experiment designs are introduced in Sect. 2. The impact of assimilating the FY-4A VIS radiance data on the analysis and first-guess forecasts of the tropical storm is discussed in Sect. 3. Finally, the conclusion and an outlook of future works are summarized in Sect. 4.

2 Models and experiment designs

The OSSE framework in this study consists of a nature run, a control run, a cluster of single-observation experiments, and a group of cycled DA experiments. The nature run was performed to generate a proxy true atmosphere state. The DA experiments that assimilated the simulated FY-4A VIS radiance data were conducted to explore the impact of the FY-4A VIS radiance DA on the tropical storm forecast. A control run that excluded DA was carried out for comparison. The OSSE was performed based on the tropical storm Higos in 2020.

On 16 August 2020, a tropical disturbance occurred over the north of Luzon, the Philippines. The system tracked northwest toward the South China Sea, and intensified into a tropical storm at 19:00 UTC on 17 August. The tropical storm developed into a typhoon system named “Higos” at 12:00 UTC, 18 August. Higos landed on Zhuhai, Guangdong province at 22:00 UTC, 18 August, and weakened into a tropical depression at 12:00 UTC, 19 August . This study focuses on the pre-landfall stage of the tropical storm (00:00–12:00 UTC, 18 August) considering that FY-4A VIS imagery is only available at daytime. During this period, the tropical storm had multi-layer and mixed-phase cloud structures, which facilitates the evaluation of the abilities of assimilating VIS radiance data for these cloud structures.

2.1 Configurations of the WRF model

The WRF model domain settings were kept the same for the nature run, control run, and DA experiments to avoid errors caused by displacement of grids between the observation and simulations. The WRF model domain covers parts of East Asia and the Western Pacific (Fig. 1). The domain contains 151 × 177 horizontal grid boxes with a grid spacing of 15 km in the horizontal directions and 40 vertical levels, and the model top was set to 50 hPa. To avoid the disturbances over the regions close to the model domain boundaries, simulations within the inner rectangle of 131 × 157 horizontal grids are analyzed.

Figure 1 The WRF model domain with 15 km horizontal grid spacing. Only the observations within the inner white rectangle were assimilated and those close to the model domain boundaries were discarded. The white dots denote the locations that were used for single-observation experiments (discussed in Sect. 2.3.2). Detailed discussions are given for points 1 (green), 2 (blue), and 3 (red) in the text. The two white lines, AB and CD, are for cross-section analyses.

For the nature run, the initial conditions (ICs) and lateral boundary conditions (LBCs) were extracted from the National Centers for Environmental Prediction (NCEP) Final (FNL) Operational Global Analysis data (with 1° × 1° resolution, available at https://doi.org/10.5065/D6M043C6, National Centers for Environmental Prediction et al., 2000). The WRF model configurations include the Thompson microphysical scheme (Thompson et al., 2008), the Tiedtke cumulus parameterization option (Tiedtke, 1989; Zhang et al., 2011), and the University of Washington (UW) planetary boundary layer scheme (Bretherton and Park, 2009). These are the optimal schemes for typhoon simulations over the Northwest Pacific Ocean as suggested by Di et al. (2019). The other model configurations include the revised MM5 Monin–Obukhov surface layer scheme (Jiménez et al., 2012), the five-layer thermal diffusion land surface scheme (Dudhia, 1996), and the rapid radiative transfer model for global climate models (RRTMG) longwave and shortwave radiation schemes (Iacono et al., 2008). The Thompson microphysical scheme provides prognostic variables of liquid water particles, including cloud water droplets and rain drops, and ice particles, including ice crystals, snow, and graupel. The nature run was initialized with a cold start at 12:00 UTC, 17 August 2020. To exclude a spin-up time of 14 h, the WRF model simulations between 02:00 and 12:00 UTC, on 18 August 2020 were used as a proxy true atmosphere of the tropical storm. The nature run captured the track and general properties of Higos. Synthetic observations of FY-4A VIS radiance were simulated with RTTOV, which is described in detail in Sect. 2.2. The simulated VIS imagery (15 km × 15 km resolution) was approximately equivalent to the superobbed 2 km resolution imagery (like those provided by FY-4A real observations) by averaging the 2 km × 2 km imagery for every block of about 7 × 7 pixels. Because the observation locations and model grid points are overlapped, the locations of the synthetic observations are directly assigned to the model grid points without interpolation during the DA processes.

For the cycled DA experiment, the ensemble size is set to 40 and the ICs and LBCs were extracted from the ERA5 hourly data at 0.25° × 0.25° resolution (the ERA5 hourly data (Hersbach et al., 2020) on pressure levels and on single levels are available at https://doi.org/10.24381/cds.bd0915c6 (Hersbach et al., 2018a) and https://doi.org/10.24381/cds.adbb2d47 (Hersbach et al., 2018b) respectively). Perturbations, which were extracted based on the WRF 3DVAR system using a generic background error option with proper scaling, were added to the ICs. The scaling factors for the variance, horizontal length scale, and vertical length scale are set to 0.25, 1.0, and 1.5, respectively. To avoid discontinuities and poor results at the boundary, LBCs at each analysis time were updated based on the analysis and WRF lateral boundary conditions using an approach built in the DART pert_wrf_bc module. This is why we choose the higher-resolution LBCs, as will be done in real DA applications, for the DA experiments than for the nature run. The WRF model microphysics configurations are the same as the nature run. The ensemble members were initialized by a cold start at 00:00 UTC on 18 August 2020. After a spin-up of 2 h, synthetic visible radiance observations were assimilated to the ensemble members from 02:00 to 09:00 UTC, 18 August 2020. The time span corresponds to the daytime when VIS imagery is available. The ensemble members were advanced to 12:00 UTC. The updating frequency of the first-guess state variables was set according to different experiment designs summarized in Table 1. With these set-ups, the effects of assimilating VIS imagery on the spin-up of the WRF model and on the analysis and first-guess forecasts of the state variables including cloud and precipitation were explored. The model settings for the control run are the same as the cycled DA experiments, except that no observations were assimilated.

2.2 Configurations of the RTTOV radiative transfer package

Synthetic AGRI channel 2 radiance was simulated based on WRF outputs using the RTTOV radiative transfer package. The input parameters of RTTOV include cloud-related parameters (the vertical structures of liquid water mixing ratio, ice water mixing ratio, cloud water effective radius, cloud fraction, etc.), atmosphere profiles (the water vapor mixing ratio profile, temperature profile, etc.), surface properties (elevation, surface type, etc.), and sun-satellite viewing geometries, etc. In the partly cloudy regions, the top-of-atmosphere (TOA) radiance is a weighted average of the radiance for a clear sky and a cloudy sky. The weight attached to the cloudy sky (i.e., the effective cloud fraction) was calculated as the vertically averaged cloud fraction weighted by the mixing ratio of hydrometeors (Geer et al., 2009). It is noted that cloud fraction (CFC) parameterization in the WRF model depends on relative humidity (RH), the saturation water vapor mixing ratio (q^∗), and cloud water + ice mixing ratios (q_l+i) (Xu and Randall, 1996),

$$\begin{array}{}\text{(1)}& \mathrm{CFC}=\left\{\begin{array}{l}{\mathrm{RH}}^{\mathrm{p}}\left[\mathrm{1}-\exp \left(\frac{-\mathit{\alpha }{q}_{\mathrm{l}+\mathrm{i}}}{{\left[\left(\mathrm{1}-\mathrm{RH}\right)q^{\ast }\right]}^{\mathit{\gamma }}}\right)\right],\text{if}\mathrm{RH}<\mathrm{1}\\ \mathrm{1},\text{if}\mathrm{RH}\ge \mathrm{1}\end{array}\right.,\end{array}$$

where p, α, and γ are suggested to be 0.25, 100, and 0.49, respectively.

The solar zenith angle, solar azimuth angle, satellite viewing zenith angle, and satellite azimuth angle were calculated using the Python astropy library according to the UTC time and the FY-4A satellite position (104.7° E). In addition, errors ranging from 1 to 4 mW m⁻² sr⁻¹ were assigned to the observations (Table 1). RTTOV includes the schemes considering different predefined cloud optical properties. For liquid water clouds, the “Deff” scheme where cloud optical properties are parameterized in terms of r_e (Mayer and Kylling, 2005) was used. The cirrus scheme developed by Baren et al. (2014) was used to calculate ice cloud optical properties, which has no explicit dependence on ice particle size. Therefore, analyses of the results were simplified given that cloud variables were adjusted collectively, but we do not have to analyze the effective radius of ice particles.

The radiative transfer processes were simulated by the DOM solver in RTTOV. The surface was treated as a specular reflector for downwelling emitted radiance. For land surface, the surface bidirectional reflectance distribution function (BRDF) was drawn from the land surface atlases (Vidot and Borbás, 2014; Vidot et al., 2018). For sea surface, BRDF was calculated with the JONSWAP (Hasselmann et al., 1973) solar sea BRDF model. The lay-to-space transmittance was computed by the v9 predictor on 54 levels (Matricardi, 2008). The downwelling atmospheric emission was computed using the linear-in-tau approximation for the Planck source term. Water vapor profiles were drawn from the WRF outputs. Other parameters not explicitly mentioned are set to the default values in DART/RTTOV.

Based on the above model configurations, the dependence of AGRI channel 2 radiance on the cloud water path (CWP) and effective radius of cloud water droplet (r_e) is presented in Fig. 2. CWP denotes the vertically integrated cloud liquid and ice water mixing ratio in an atmospheric column, which is calculated by

$$\begin{array}{}\text{(2)}& \mathrm{CWP}={\int }_{P_{\mathrm{s}}}^{P_{\mathrm{t}}}{\displaystyle \frac{\mathrm{1}}{g}}(Q_{\mathrm{c}}+Q_{\mathrm{i}})\mathrm{d}P,\end{array}$$

where P_s and P_t denote the surface and model top pressures. Q_c and Q_i are the liquid water mixing ratio (the sum of the mixing ratio of cloud droplet and rain) and ice water mixing ratio (the sum of the mixing ratio of ice, snow, and graupel), and g is the gravitational acceleration (9.8 m s⁻²).

Figure 2 Dependence of AGRI channel 2 radiance on (a) cloud water path (CWP) and effective radius (r_e​​​​​​​), and (b) cloud fraction (CFC) for CWP of 10 kg m⁻² and r_e of 15 µm. The simulation was performed with the “Deff” scheme for liquid water cloud optical properties and the Baren et al. (2014) scheme for cirrus optical properties. The solar zenith angle, viewing zenith angle, and relative azimuth angle were set to 25, 40, and 135°, respectively.

The curvature properties in Fig. 2 indicate a non-linear relationship between the observation (radiance) and cloud parameters (CWP and r_e). The variation of the radiance-CWP functions with different effective radii becomes smaller as r_e increases. For r_e larger than 30 µm, the radiance-CWP functions almost do not change with effective radii. Because raindrops are several orders larger than cloud droplets, the r_e of cloud droplets is sufficient to describe the radiative transfer processes for the clouds where cloud droplets and raindrops coexist. As a result, r_e in the following discussion explicitly denotes the effective radius of cloud droplets, which corresponds to the WRF state variable “RE_CLOUD”.

2.3 DA experiment design and DART configurations

2.3.1 DART filters

DART was configured to employ the ensemble adjustment Kalman filter (EAKF, Anderson, 2001) and the rank histogram filter (RHF, Anderson, 2010) for this study. EAKF and RHF are two variants of the deterministic filters. Therefore, no perturbations were added to the observations. EAKF is a serial ensemble DA algorithm and the observations are assimilated as scalars. The model state variable x_m is updated by Eq. (3) (Anderson, 2001),

$$\begin{array}{}\text{(3)}& {\mathit{x}}_m^{\prime }={\mathit{x}}_m+\mathrm{\Delta }{\mathit{x}}_{m,n},m=\mathrm{1},\mathrm{\dots },M,n=\mathrm{1},\mathrm{\dots },N,\end{array}$$

where x_m denotes the mth state variable, ${\mathit{x}}_m^{\prime }$ the updated value of x_m, and Δx_m,n the state variable increment for the mth state variable due to the nth observation. Δx_m,n is calculated by Eq. (4),

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{\mathit{x}}_{m,n}=\left({\mathit{\sigma }}_{\mathrm{p},m}/{\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}}\right)\mathrm{\Delta }{y}_n,m=\mathrm{1},\mathrm{\dots },M,n=\mathrm{1},\mathrm{\dots },N,\end{array}$$

where the subscript “p” stands for the prior estimate (i.e., the first guess), σ_p,m is the first-guess sample error covariance between the observation and the mth state variable x_m, and ${\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}}$ is the first-guess sample error variance of the observed variable. Δy_n is the observation increment for the nth ensemble, which is calculated by the following equation:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{y}_n=\left(y_n^{\mathrm{p}}-{\stackrel{\mathrm{‾}}{y}}_{\mathrm{p}}\right)\left({\mathit{\sigma }}_{\mathrm{u}}/{\mathit{\sigma }}_{\mathrm{p}}\right)+{\stackrel{\mathrm{‾}}{y}}_{\mathrm{u}}-y_n^{\mathrm{p}},n=\mathrm{1},\mathrm{\dots },N,\end{array}$$

where $y_n^{\mathrm{p}}$ denotes the first guess of the observed variable for the nth ensemble, ${\stackrel{\mathrm{‾}}{y}}_{\mathrm{p}}$ the ensemble mean of the first guess of the observed variable, ${\stackrel{\mathrm{‾}}{y}}_{\mathrm{u}}$ the ensemble mean of the posterior estimate (i.e., the analysis) of the observed variable, and σ_u the updated standard deviation of σ_p. ${\stackrel{\mathrm{‾}}{y}}_{\mathrm{u}}$ and σ_u are calculated by Eqs. (6) and (7):

$$\begin{array}{}\text{(6)}& {\displaystyle }\stackrel{\mathrm{‾}}{y_{\mathrm{u}}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{y}}_{\mathrm{p}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}}}}{y}_{\mathrm{o}},\text{(7)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{u}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{o}}{\mathit{\sigma }}_{\mathrm{p}}}{\sqrt{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}}}}},\end{array}$$

where y_o and σ_o denote the observation and its corresponding observational error standard deviation.

Anderson (2007, 2009) promoted a spatially varying state–space adaptive covariance inflation to the first-guess state to increase the prior ensemble spread. The same option was adopted in this study and in several other papers (Lei et al., 2015; Kurzrock et al., 2019). The adaptive inflation uses 1.0, 0.6, and 0.9 as the initial value, fixed standard deviation, and damping settings, respectively. The sampling error due to the use of the limited ensemble size was corrected with the method developed by Anderson (2012). Since observations such as satellite VIS radiance data do not have a specific single vertical location, no vertical localization was used in this study.

The RHF produces a posterior ensemble based on a continuous approximation of the prior probability density function (PDF) and a given likelihood function. The prior PDF is approximated by a rank histogram which has a piecewise constant between two ensemble members and follows Gaussian distributions beyond the lower and upper bounds of the ensemble members. The posterior distribution is calculated by the Bayes theorem, and the state variable is updated by searching for the appropriate position in the state variable space which partitions the posterior distribution to unity probability for each ensemble member. The prior PDF does not have to respect the Gaussian form for RHF. Therefore, the method is declared to be more suitable for non-Gaussian problems. Details on this algorithm can be found in Anderson (2010).

2.3.2 Single-observation experiments

With the OSSE set-ups, a set of single-observation experiments were performed by employing EAKF to assimilate the VIS radiance data. The single-observation experiments assimilate an observation at a given pixel, and the adjustment of state variables in the column at a targeting pixel is only caused by assimilating the one observation. This is convenient for evaluating the basic functionality of assimilating VIS radiance data. The potentials, inabilities, and ambiguities of assimilating VIS radiance data are discussed herein. The single-observation experiments were performed at 02:00 UTC on 18 August 2020, and no forecast was carried out. The affected cloud variables include Q_c, Q_i, r_e, CFC, and the non-cloud variables including water vapor mixing ratio (Q), perturbation potential temperature (T), and the x- and y-wind components (U and V).

The single-observation experiments were performed for the most inner parts of the satellite imagery to avoid disturbances near boundaries. The observations at 02:00 UTC were thinned by selecting every six pixels to ensure that the selected observations are far from each other. This resulted in 176 points shown in Fig. 1. By setting a localization distance of 15 km, assimilating the VIS radiance at a pixel would not influence the state variables at the surrounding pixels. Therefore, like Scheck et al. (2020), we performed a cluster of 176 single-observation experiments in one DA cycle to save computational cost. Among the 176 selected pixels, special focus was placed on the three colored points, which were designed to illustrate the ambiguities related to cloud layered structures and cloud phases and to show the limitations caused by the non-Gaussian and non-linear problems.

2.3.3 Cycled DA experiments

In total, 14 cycled DA experiments were performed to evaluate the influence of different model settings and observation pre-processing on the analysis and forecast with VIS radiance DA. The purpose of the cycled DA experiments is to reveal the forecast quality and growth of the forecasting errors during assimilating satellite VIS radiance data, and to provide guidance on the use of WRF/DART-RTTOV with an outlook to future applications. The experiment set-ups cover the tests of different filter algorithms, cycling intervals, cycling variables, outlier threshold values, observation errors, and observations with or without data-thinning. The outlier threshold value for the observation is a predefined threshold for rejecting an observation depending on its distance from the ensemble mean of the first guess. If the distance is more than N (the predefined outlier threshold value) standard deviations from the square root of the sum of the first-guess ensemble and observation error variance, the observation is rejected. A description of the experiment designs is summarized in Table 1.

The comparison between Exp-01–Exp-03 and Exp-4–Exp-06 groups was designed to reveal the pros and cons of EAKF and RHF in the WRF analysis and forecast. The comparison between Exp-01 and Exp-02 and between Exp-07 and Exp-08 groups was designed to reveal the influence of updating the thermal and dynamic variables. The comparison between Exp-03, Exp-09, and Exp-10 was designed to reveal the influence of the observation errors. The comparison between Exp-01–Exp-03 and Exp-11–Exp-13 groups was designed to reveal the influence of the outlier threshold values. Finally, the comparison between Exp-10 and Exp-14 was designed to reveal the influence of observation thinning.

Table 1 Parameter settings for the cycled data assimilation experiments. x_cloud denotes the WRF cloud variables including cloud fraction (CLDFRA), the mixing ratio of cloud droplet (QCLOUD), rain (QRAIN), ice (QICE), snow (QSNOW), graupel (QGRAUP), the effective radius of cloud water droplet (RE_CLOUD), and the effective radius of cloud ice droplet (RE_ICE). x_atmos denotes the WRF non-cloud variables including water vapor mixing ratio (QVAPOR), water vapour mixing ratio at 2 m height (Q2​​​​​​​), x-, y-, and z-wind components (U, V, W), x- and y-wind components at 10 m height (U10 and V10), temperature at 2 m height (T2), perturbation geopotential (PH), perturbation potential temperature (T), perturbation dry air mass in column (MU), and surface pressure (PSFC).

2.4 Metrics of simulation errors

Root mean square error (RMSE) and mean absolute error (MAE) are two of the most widely used metrics in assessing weather simulation errors (Kurzrock et al., 2019). RMSE is much more sensitive to extremely large errors than MAE. For satellite VIS radiance DA, some exceptionally large analysis increments of CWP were rarely expected (details provided in Sect. 3), implying that the difference in RMSE between the first guess and the analysis was not as distinct as MAE. Thus, MAE is used to measure the difference between the simulated CWP and the theoretical true CWP (derived from the nature run). MAE is calculated by the following formula:

$$\begin{array}{}\text{(8)}& \mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{n_x{n}_y}}{\sum }_{i,j}\mathrm{abs}(x_{i,j}^{\mathrm{sim}}-x_{i,j}^{\mathrm{obs}}),\end{array}$$

where $x_{i,j}^{\mathrm{sim}}$ ($x_{i,j}^{\mathrm{obs}}$) denotes the simulated (true) CWP at the ith (in the zonal direction) and jth (in the meridional direction) model grid. n_x and n_y denote the number of pixels in zonal and meridional directions, respectively, of the relevant model domains.

A fraction skill score (FSS) was developed to measure the accuracy of spatially inhomogeneous variables at a specific spatial scale. Therefore, it can mitigate the “double-penalty” problem (Mittermaier et al., 2013) for small spatial shifts in features of interest. FSS is calculated by the following formula:

$$\begin{array}{}\text{(9)}& \mathrm{FSS}=\mathrm{1}-{\displaystyle \frac{\frac{\mathrm{1}}{m_x{m}_y}{\sum }_{i,j}{\left(p_{i,j}^{\mathrm{obs}}-p_{i,j}^{\mathrm{sim}}\right)}^{\mathrm{2}}}{\frac{\mathrm{1}}{m_x{m}_y}({\sum }_{i,j}{p}_{i,j}^{\mathrm{obs}}+{\sum }_{i,j}{p}_{i,j}^{\mathrm{sim}})}},\end{array}$$

where $p_{i,j}^{\mathrm{obs}}$ denotes the cloud fraction within a subdomain covering 3 × 3 model grids. m_x and m_y denote the dimensions of subdomains in the zonal and meridional directions, respectively.

For evaluation of the precipitation simulation, the threat score (TS) is used. TS is computed as

$$\begin{array}{}\text{(10)}& \mathrm{TS}={\displaystyle \frac{H}{F+M+H}},\end{array}$$

where H denotes the number of pixels with the correct representation of precipitation (hits), F denotes the number of pixels where simulation indicates precipitation while the true state indicates non-precipitation (false alarms), and M denotes the number of pixels where simulation indicates non-precipitation while the true state indicates precipitation (under predictions).

Following Scheck et al. (2020), we also use the mean profile error (MPE, ε) to assess the error of the model state with respect to the nature run. If the difference between the MAE of the analysis (ε_pos) and that of the first guess (ε_pri), calculated as $\mathit{\delta }\mathit{\epsilon }={\mathit{\epsilon }}_{\mathrm{pos}}-{\mathit{\epsilon }}_{\mathrm{pri}}$, is negative, a positive impact is generated by the DA procedure, and vice versa.

3 Results

3.1 Single-observation experiments

The results discussed herein correspond to the OSSE set-ups described in Sect. 2.3.2. Only the ensemble mean of the first guess and analysis of state variables are discussed. We focus on three cases. (1) In both state variable (or a diagnosed parameter such as CWP) and observation spaces, the analysis is within the range bounded by the first guess and the truth; (2) the analysis is within the first guess and the truth in the observation space but not in the state variable space. This case is associated with the spurious covariance and non-linear properties of the forward operator; (3) The analysis is beyond the first guess and the truth in both observation space and state variable space, which is closely related to the non-Gaussian properties of the prior PDF. The results including the three cases are shown in Fig. 3.

Assimilating VIS radiance data generated a nonsignificant impact on the non-cloud variables including the x- and y-wind components (U and V in Fig. 3a) and the temperature and water vapor mixing ratio (T and Q in Fig. 3b). From the perspective of radiative transfer, the VIS radiance is insensitive to U and V at the analysis time. Therefore, the adjustment in observation space should not influence U and V. In addition, the VIS radiance is closely related to CFC (Fig. 2b), and thus an implicit relationship between the VIS radiance and RH could be expected because the parameterization of CFC involves water vapor (Eq. 1). The VIS radiance is positively related to CFC. However, given that RH depends not only on Q, but also on T and pressure, spurious covariance between the VIS radiance and $Q/T$​​​​​​​ may be generated due to the ensemble spread of $Q/T$. The ensemble spread of $Q/T$ for the ensemble members would blur the relationship between VIS radiance and Q or T. Therefore, only a neutral impact on Q and T was revealed for the single-observation experiments.

Figure 3 Differences in the mean profile error (MPE) between the ensemble mean of the first guess and the analysis, denoted as $\mathit{\delta }{\mathit{\epsilon }}_X=X_{\mathrm{MPE},\mathrm{pos}}-X_{\mathrm{MPE},\mathrm{pri}}$, where X denotes a state variable or a diagnosed parameter, and pos and pri the analysis and first guess, respectively. The variables include (a) the x- and y-wind components (U and V), (b) the perturbation potential temperature (T) and water vapor mixing ratio (Q), (c) the cloud fraction (CFC) and effective radius (r_e), and the (d) cloud water path (CWP) and radiance (RAD). The plus signs in green, blue, and red correspond to points 1, 2, and 3, respectively, in Fig. 1.

Assimilating VIS radiance data improved the ensemble mean of r_e, CFC, CWP, and VIS radiance for most points in Fig. 3c and d. Let us take point 1 marked in Fig. 1 as an example for case 1. The profiles of the ensemble mean of cloud and non-cloud variables for the first guess and analysis are shown in Fig. 4. Point 1 corresponds to a single ice cloud layer between 400 and 200 hPa with a CWP of 0.01 kg m⁻² and a TOA radiance of 3.63 mW m⁻² Sr⁻¹ for the nature run. The first guess features a two-layer mixed-phase cloud with a false-alarm liquid water cloud layer below 500 hPa. The ensemble mean first-guess CWP and the equivalent VIS radiance are 1.33 kg m⁻² and 7.29 mW m⁻² Sr⁻¹, respectively. After assimilating the satellite VIS radiance data, the first guess was drawn toward the nature run both in the observation space, with a decreased ensemble mean radiance of 6.40 mW m⁻² Sr⁻¹, and in CWP, with a decreased ensemble mean CWP of 0.85 kg m⁻². As a result, Q_c, Q_i, CFC, and r_e were adjusted collaboratively toward the nature run.

Figure 4 The vertical profiles of state variables for the nature run (theoretical truth), the first guess and analysis (ensemble mean). Q_c denotes the liquid water mixing ratio, Q_i the ice water mixing ratio, CFC the cloud fraction, r_e the effective radius of liquid water droplets, QVAPOR the water vapour mixing ratio, δT the perturbation potential temperature, U and V the x- and y-wind components.

Since the simulated VIS radiance observation is not sensitive to cloud vertical structures but to the accumulated cloud water/ice mass, assimilating the VIS radiance observation could not reduce cloud vertical location errors. In addition, it could not remove the false-alarm liquid clouds produced by the spurious covariance between the VIS radiances and liquid water clouds in the background. According to Eq. (4), the analysis increment of each state variable is linearly related to its covariance with observations. Therefore, the vertical structures and hydrometeor phases of the analysis are mainly determined by those of the first guess. A larger first guess of the state variable would generate larger covariance, and a larger adjustment to the first guess would be expected. Because the ensemble mean of the first-guess Q_c and Q_i was larger in the lower layer (≥400 hPa), the adjustment of Q_c or Q_i was much more distinct for the lower layer than the upper layer (≤300 hPa). Similar results were also found for r_e except that larger liquid water particles occurred in the middle layer (∼ 600 hPa) and smaller liquid water particles occurred in the lower layer (∼ 800 hPa). The covariance between CFC and the synthetic observation is 0 in the upper layer (≤200 hPa) because CFC is almost a constant of 1 for all ensemble members (the ensemble spread of CFC is 0). Compared with the cloud variables, the non-cloud variables remain almost unchanged after the DA.

Case 2 is used to illustrate that positive impact on VIS radiance does not ensure a positive impact on CWP. Let us take point 2 in Fig. 1 as an example. The theoretical truth, the ensemble mean of the first guess and the analysis of radiance/CWP, are 8.59 mW m⁻² Sr${}^{-\mathrm{1}}/$0.17 kg m⁻², 7.94 mW m⁻² Sr${}^{-\mathrm{1}}/$2.50 kg m⁻², and 8.00 mW m⁻² Sr${}^{-\mathrm{1}}/$2.65 kg m⁻², respectively; that is, the analysis of CWP is beyond the range bounded by the first guess and the truth. This is partly due to the non-linear relationship between CWP and VIS radiance. To elaborate this problem, we calculated the ensemble mean with Eq. (4), and substituting $\stackrel{\mathrm{‾}}{y_{\mathrm{u}}}$ in Eq. (6) would give the following formula:

$$\begin{array}{}\text{(11)}& \stackrel{\mathrm{‾}}{\mathrm{\Delta }{x}_m}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{p},m}}{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{p}}^{\mathrm{2}}}}{R}_{\mathrm{inc}},\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathrm{\Delta }{x}_m}$ denotes the ensemble mean of the mth state variable increment and R_inc denotes the ensemble mean radiance increment, which is calculated by the following formula:

$$\begin{array}{}\text{(12)}& R_{\mathrm{inc}}=y_{\mathrm{o}}-{\stackrel{\mathrm{‾}}{y}}_{\mathrm{p}}.\end{array}$$

When considering a simplified case with just two ensemble members, the ensemble mean observation increment is calculated by the following formula:

$$\begin{array}{}\text{(13)}& R_{\mathrm{inc}}=F\left(W_{\mathrm{o}}\right)-{\displaystyle \frac{F\left(W_{\mathrm{1}}\right)+F\left(W_{\mathrm{2}}\right)}{\mathrm{2}}},\end{array}$$

where W_o​​​​​​​ denotes the observed CWP and F denotes the forward operator. W₁ and W₂ represent CWP of the two ensemble members. However, considering the relationship between CWP and the VIS radiance, the theoretical true observation increment should be

$$\begin{array}{}\text{(14)}& R_{\mathrm{inc}}^{\mathrm{t}}=F\left(W_{\mathrm{o}}\right)-F\left(\stackrel{\mathrm{‾}}{\mathit{W}}\right),\stackrel{\mathrm{‾}}{\mathit{W}}={\displaystyle \frac{W_{\mathrm{1}}+W_{\mathrm{2}}}{\mathrm{2}}}.\end{array}$$

As indicated by Fig. 4, R_inc is larger than $R_{\mathrm{inc}}^{\mathrm{t}}$; that is, the ensemble mean observation increment was overestimated by Eq. (14), leading to an overestimated ensemble mean of the analysis of CWP.

Figure 5 Illustration of the non-linear effects of the observation operator on the calculation of radiance increments with two ensemble members. F denotes the observation operator, W₁ and W₂ denote the cloud water path (CWP) for the first and second ensemble member, respectively, and W_O denotes the observed CWP.

Case 3 is to show that some negative impact could be generated in the observation space in certain conditions. Let us take point 3 in Fig. 1 as an example. The ensemble mean radiance (2.51 mW m⁻² Sr⁻¹) of the analysis is beyond the range bounded by the ensemble mean of the first guess (2.56 mW m⁻² Sr⁻¹) and the true radiance (3.41 mW m⁻² Sr⁻¹). The EAKF algorithm assumes that the prior PDF, p(x), of the equivalent observations and model state variables (or the diagnosed variable such as CWP in this study) conforms to Gaussian functions. To see how well the assumption was respected, p(x) of radiance and CWP is presented in Fig. 6; it shows non-Gaussian prior PDFs. Several studies concluded that the non-Gaussian properties negatively affect the performance of ensemble methods (Lawson and Hansen, 2004; Lei et al., 2010). Therefore, we tentatively ascribe the negative impact in the observation space partly to the non-Gaussian properties of p(x). Accordingly, DA experiments using RHF, which is not bound by the Gaussian assumption, were added to the cycled DA experiments in comparison with EAKF.

Figure 6 The first-guess probability density functions (PDFs) of (a) equivalent visible radiance and (b) the cloud water path (CWP). The PDFs were estimated from the 40 prior ensembles for point 3 (red dot) in Fig. 1.

3.2 Cycled DA experiments

The results in this section correspond to the OSSE set-ups described in Sect. 2.3.3. The main focus is the impact of the VIS radiance DA on the analysis and first-guess forecast of CWP, cloud coverage, non-cloud state variables, and precipitation.

3.2.1 Impact on CWP and cloud coverage

The time evolution of the CWP for the nature run, control run, and the first-guess forecast and the analysis of CWP for Exp-01 are presented in Fig. 7.

Figure 7 Time evolution of the cloud water path (CWP) for the nature run (column 1), control run (column 2), first-guess forecast (column 3), and analysis (column 4) of Exp-01. From top to bottom, the row panels correspond to 02:00, 04:00, 06:00, 08:00, and 10:00 UTC on 18 August 2020.

The results indicate distinct differences between the first guess and the analysis of CWP on 02:00 UTC, 18 August 2020 (Fig. 7a3–a4). After assimilating the VIS radiance data, the horizontal distribution of the ensemble mean first-guess CWP is quite similar to that of the analysis. An extremely large analysis increment of CWP was rarely expected, as mentioned in Sect. 2.4. The similarities between the cycled DA experiments and the nature run also indicate the improvements in the analysis and the first-guess forecast of CWP and cloud coverage. Compared with the control run, assimilating the VIS radiance data clearly suppressed the false-alarm clouds. However, DA of VIS radiance could not generate clouds which were underpredicted. The inability to correct the underprediction is illustrated with a cross-section analysis shown in Fig. 8.

Figure 8 The cross section of the liquid and ice water mixing ratio of Exp-01 at 02:00 UTC, 18 August 2020 for the AB (a1–a3)​​​​​​​ and CD (b1–b3) lines shown in Fig. 1. The panels correspond to the nature run (a1, b1), the ensemble mean of the first guess (a2, b2), and the ensemble mean of the analysis (a3, b3).

Assimilating VIS radiance data did not improve the underprediction of the cloud distributions in vertical and horizontal directions. In the vertical direction, a single-layer cloud was present between 4 and 12 km height in the nature run (Fig. 8a1 and b1). However, clouds were present between 0 and 12 km height in the first guess (Fig. 8a2 and b2). After the DA, $Q_{\mathrm{c}}/Q_{\mathrm{i}}$ was decreased but the clouds below 4 km height were not removed (Fig. 8a3–b3). In the horizontal direction, a thin-layer cloud was present at 14–16° N in the nature run (Fig. 8b1). This cloud fragment was not simulated by the first-guess atmosphere state (Fig. 8b2), nor was it regenerated after assimilating the VIS radiance data (Fig. 8b3). Therefore, assimilating the VIS radiance would not generate cloud hydrometeors for the region with clear sky in the first guess because there is only zero spread of cloud variables from the unrepresentative background error covariance.

Quantitative analyses of CWP and CFC indicated improved analysis and first-guess forecasts for the cycled DA experiments. The results varied according to the filter algorithms, cycling variables, outlier thresholds, observation errors, and observations with or without thinning, which are discussed in the following.

Influence of filter algorithms

To evaluate the performance of different filters on the analysis and first-guess forecast of cloud variables, quantitative analyses of FSS and MAE of the ensemble mean CWP for Exp-01–Exp-03 (EAKF) and Exp-04–Exp-06 (RHF) are presented in Fig. 9. In general, the performance of RHF is comparable to or slightly better than EAKF. At some analysis times before 03:00 UTC, the FSS of CWP for the analysis is larger for Exp-01 than for Exp-04, but the first-guess FSS at the next analysis time is larger for Exp-04 than Exp-01. Similar results were found for the 1 and 3 h first-guess forecasts. These results suggest that better analyses do not always ensure better forecasts.

Figure 9 Time evolution of (a) FSS and (b) MAE of the ensemble mean of the first-guess forecast and analysis of CWP for the cycled DA experiments which are designed to illustrate the sensitivity to filter algorithms.

Unlike EAKF that assumes a Gaussian prior PDF, RHF does not. However, the performance of RHF is also subject to sampling errors due to limited ensemble members and other factors, as indicated by Anderson (2010). Therefore, only comparable or slightly better analysis and first-guess estimates were revealed for RHF than EAKF. In addition, updating the state variables by RHF comes with more computational cost. For the DA of 20 567 observations in one DA cycle at a Linux cluster equipped with a 2.20 GHz Xeon Silver 4214 CPU with 12 cores, the elapsed CPU time is 775 and 440 s for the RHF and EAKF methods, respectively.

Influence of cycling variables

The single-observation experiments indicate that assimilating VIS radiance data only generated a nonsignificant impact on non-cloud variables at 02:00 UTC, 18 August 2020. However, the non-cloud variables are an important part of a cycling DA system. To explore the impact of including or excluding the updated non-cloud parameters in the ensemble cycles, FSS and MAE of the 10 min and 1 h first-guess forecast and analysis of the ensemble mean CWP for the Exp-01 to Exp-02 and Exp-07 to Exp-08 groups are analyzed. Figure 10 shows that Exp-01 (and Exp-02) outperforms Exp-07 (and Exp-08), indicating that including the cloud and non-cloud variables in the ensemble cycling makes the forecasting more skillful than cycling the cloud variables alone. The benefit was introduced to the non-cloud variables during the model integration.

Figure 10 Time evolution of (a) FSS and (b) MAE of the ensemble mean of the first-guess forecast and analysis of CWP for the cycled DA experiments designed to illustrate the sensitivity to cycling variables.

To demonstrate the error growth for the non-cloud variables, the temporal evolution of the ratio of benefits , calculated by ​​​​​​​Eq. (15), is presented in Fig. 11:

$$\begin{array}{}\text{(15)}& \textcolor{red}{}=N_{\mathrm{bet}}/N_{\mathrm{eff}},\end{array}$$

where N_bet denotes the number of horizontal grid boxes with negative differences of MPE between the analysis and the first guess (refer to Sect. 2.4 and Scheck et al., 2020), and N_eff denotes the number of observations effectively assimilated by the DA system (see the next section).

Figure 11 Time evolution of the ratio of beneficial impact , calculated by ​​​​​​​Eq. (15) for the first-guess forecast and analysis of (a) the perturbation potential temperature, (b) the water vapor mixing ratio, (c) the x-wind component, and (d) the y-wind component.

Figure 11 indicates the positive impact of SW radiance DA on the major non-cloud variables, especially at the later cycling steps. We believe that the main reason for this is the positive feedback to the non-cloud variables due to the adjustment to cloud variables during the WRF ensemble integration.

A potential explanation for the slightly positive impact of the DA on the water vapor mixing ratio is given here. If the ensemble mean of the first-guess equivalent radiance is overestimated for a cloud containing liquid hydrometeors, Q_c tends to be decreased in the analysis field in order to match a negative analysis increment in the observation space. In the next ensemble forecast cycle, the grid boxes with decreased Q_c may generate more hydrometeors due to condensation. This is because the supersaturation of the surrounding atmosphere is increased due to a loss of some cloud hydrometeors (i.e., a loss of condensation surface). The increased supersaturation could activate new condensation nuclei and enhance the condensation process, which consumes more water vapor. As a result, the VIS radiance could be positively related to the water vapor mixing ratio, and vice versa. According to Eq. (4), the covariance between the VIS radiance and water vapor mixing ratio will adjust Q correctly.

The adjustment of temperature is more likely related to the interactions between clouds and radiation. For example, decreased CWP or CFC tends to enhance the direct radiation flux in the surface layer, increasing the low-level temperature toward the truth´, as indicated by Scheck et al. (2020), and vice versa. In addition, the interactions between clouds and longwave radiation tend to generate cooling effects at cloud top and heating effects at cloud bottom (Zhang et al., 2020). Therefore, a relationship between cloud radiance and temperature is expected, and the covariance between VIS radiance and temperature could adjust temperature correctly.

The impact of SW radiance DA on the x- and y-wind components is slightly negative before 04:00 UTC, and becomes positive thereafter. We believe that the positive impact is mainly caused by the convergence and divergence related to the thermal instability, which is closely related to cloud formation (increased radiance) and dissipation (decreased radiance) for the convective systems. As the cloud–radiation interactions modify the temperature profile, the change of clouds could strengthen or weaken the thermal instability and impact the z-wind component. The z-wind component is closely related to horizontal x- and y-wind components by adjusting the convergence and divergence (White et al., 2018). Therefore, an indirect “radiance–cloud–vertical velocity–convergence and divergence–horizontal wind” interaction could map the observation increment to the U and V increments.

Influence of outlier threshold values

Not all observations were effectively assimilated by the WRF/DART-RTTOV system. Some of the observations were rejected by the DA system due to two reasons: (1) existence of non-monotonic pressures, i.e., pressure increases with altitude. This situation was generated at some points by the improper interpolation of the perturbed first-guess model state to the RTTOV predefined layers. For this case study, non-monotonic pressures were mainly located in the Qinghai–Tibet Plateau, Tianshan Mountain, and Central Mountain Range in Taiwan (not shown for simplicity), where complex high terrain exists; (2) the differences between the observations and the first-guess ensemble mean equivalent observations are too large and assimilating these data may cause the WRF model to collapse. For the observations rejected due to the second reason, the ratio of observations to be assimilated is mainly determined by setting an outlier threshold value. Increasing the outlier threshold value could increase the ratio of observations effectively assimilated (Fig. 12), but it may introduce unstable adjustments to model state variables and may destroy the forecast. Therefore, a trade-off between large outlier threshold value and the potentially detrimental effects on forecasts should be assessed. The analysis and 10 min, 1 h, and 3 h first-guess forecasts indicate improved results when increasing the outlier threshold values (Fig. 13); an outlier threshold value of 6 does not cause the collapse of the WRF and leads to improvements in the analysis and first-guess forecasts of CWP and cloud coverage.

Figure 12 Time evolution of the ratio of observations that are assimilated or rejected by the WRF-DART system for different experiment designs with different outlier threshold values and observation errors.

Figure 13 Time evolution of (a) FSS and (b) MAE of the ensemble mean of the first-guess forecast and analysis of CWP for cycled DA experiments with different outlier threshold values.

Influence of observation errors and thinning

The ratio of observations effectively assimilated is also influenced by the observation errors. Enlarging the observation errors will retain a larger ratio of observations effectively assimilated (Fig. 12). However, increasing observation errors implies less weight assigned to the observations during the DA. Therefore, the analysis and first-guess forecast are influenced by the observation errors in both ways with a trade-off. To evaluate the influence, the FSS and MAE of the ensemble mean CWP for the cycled DA experiments with different observation errors are presented in Fig. 14. The analysis and first-guess forecasts of CWP and cloud coverage are negatively related to the observation errors. As complementary to the parameters controlling the number of observations, the influence by thinning observations (Exp-14) is presented. The results indicate that thinning observations led to slight improvements only.

Figure 14 Time evolution of (a) FSS and (b) MAE of the ensemble mean of the first-guess forecast and analysis of CWP for the cycled DA experiments with varying observation errors and observation pre-processing (with and without thinning).

3.2.2 Impact on precipitation

The first-guess forecast of the rain rate for the nature run, control run, Exp-01, and Exp-11 is shown in Fig. 15. On the domain average, the rain rates were overestimated by both control and cycled DA experiments. Compared with the control run, the precipitation for the cycled DA experiments was decreased in most cases, and the spatial distribution of the precipitation agreed better with the nature run (Fig. 16a).

Figure 15 Time evolution of the rain rate for the nature run (column 1), control run (column 2), first-guess forecast of Exp-01 (column 3), and Exp-11 (column 4). From top to bottom, the row panels correspond to the results for 02:00–02:10, 04:00–04:10, 06:00–06:10, 08:00–08:10, and 10:00–10:10 UTC on 18 August 2020.

Figure 16 Time evolution of (a) the TS score and (b) MAE of the ensemble mean of the first-guess forecast for rain rate.

The metrics of the rain rate forecast indicate that the forecasting skills were improved at most analysis times (Fig. 16b). However, the improvements in rain rate did not happen at all times. For example, at the initial cycling step (before 04:00 UTC, 18 August 2020), the control run seemed to outperform the cycled DA experiments. As time progressed, the advantages of assimilating VIS radiance data became apparent. The improved TS score was closely related to the improved CWP and cloud coverage. It is noted that the improvements in the rain rate were not as significant as CWP. Some studies indicated that precipitation was closely related to cloud vertical structure (Kubar and Hartmann, 2008; Yan and Liu, 2019​​​​​​​), the presence and distribution of liquid and ice hydrometeors (Field and Heymsfield, 2015; Mülmenstädt et al., 2015; Korolev et al., 2017), surrounding atmosphere, and dynamic state variables (Kanji et al., 2017), etc. Therefore, the impact of the SW radiance DA on the rain rate suffers from its inability to constrain cloud vertical structures, improve cloud phase simulations, and correct cloud location errors (underestimation), etc., as discussed in the previous sections.

4 Conclusions and future work

In this study, a series of single-observation experiments and cycled DA experiments were performed in an OSSE framework to evaluate the WRF/DART-RTTOV system for assimilating FY-4A/AGRI VIS (channel 2) radiance data. The single-observation experiments were designed to reveal the positive impact and limitations of assimilating satellite VIS radiance data on the cloud variables (liquid / ice water mixing ratio and CWP, effective radius of liquid water droplets, and cloud fraction) and non-cloud variables (water vapor mixing ratio, perturbation potential temperature, and wind). The cycled DA experiments were designed to explore the impact of assimilating VIS radiance data on the analysis and first-guess forecasts of a tropical storm case with different model settings and observation pre-processing, including varying the filter algorithms, cycling variables, updating frequencies, outlier thresholds, observation errors, and observation thinning. The main findings are as follows.

1. Single-observation experiments show that assimilating the satellite VIS radiance data generated an overall positive impact on cloud variable analysis for the first DA cycle, although in some rare cases, the DA increased the errors both in the observation space and in the cloud variable (or diagnosed parameter CWP) space. The negative impact was closely related to the non-linear properties of the forward operator and the non-Gaussian properties of the first-guess PDF. In addition, a neutral impact was revealed on non-cloud variables including water vapor mixing ratio, temperature, and horizontal winds.

2. Although a neutral impact was revealed for the non-cloud variables in the first DA cycle, both the non-cloud variables and cloud variables were improved in the following-on ensemble forecast and DA cycles. The beneficial impact on the non-cloud variables was closely related to the analysis increments of cloud variables by condensation/evaporation and deposition/sublimation processes, and to the cloud–radiation interactions of the model microphysics. The RHF filter slightly outperformed the EAKF filter thanks to its ability to deal with the non-Gaussian problems, but at a cost of ∼ 1.8 times more computational time.

3. The cycled DA experiments revealed that the first-guess forecast and analysis are positively related to increasing the outlier threshold value and the updating frequency, and negatively related to increasing observation errors and thinning length scale. Similar results were found for the precipitation forecasts. Assimilating the SW radiance data effectively improved the representation of locations with or without precipitation, but less for the quantitative metrics of the rain rate. The limited impact on the rain rate is due to the fact that the DA scheme is not effective for constraining cloud vertical structures, modifying cloud phases, and correcting cloud location errors in the case of underprediction that caused the zero-spread problem.

The findings from this study provide useful guidance on setting WRF/DART-RTTOV configurations and pre-processing of observations for assimilating the FY-4A and the upcoming FY-4B VIS radiance data in real-time models. This study complements the work by Scheck et al. (2020), who reported the potentials and limitations of assimilating the Meteosat SEVERI VIS imagery with the COSMO/KENDA system using LETKF. In general, the pros and cons of assimilating the VIS radiance data found in this study are in agreement with those reported by Scheck et al. (2020), except that a slightly positive impact on horizontal wind speeds was found in this study while Scheck et al. (2020) reported a neutral impact. We tentatively ascribe the positive impact on U and V to the convergence or divergence related to the “radiance–cloud–vertical velocity–convergence and divergence–horizontal winds” interaction, which should differ in weather systems. Therefore, the different impact on horizontal winds of the two studies could be caused by the differences in the weather systems studied. In addition, the two studies used different models/tools and corresponding configurations and, moreover, this study explored the properties unique to the WRF/DART system.

It is noteworthy that the present study was based on the low-resolution model simulations at 15 km × 15 km. Such grid spacing is large enough to avoid radiance simulation errors caused by 3D radiative effects, which are apparent for high-resolution simulations (Várnai and Marshak, 2001). Although the 3D radiative transfer effects could be properly corrected by some of the forward operators in RTTOV, the related parameters and datasets specific to FY-4A/AGRI are currently unavailable. Further studies should be extended to cloud-resolving model simulations, as in the study by Scheck et al. (2020), to fully take advantage of high-resolution satellite VIS radiance data. Moreover, some attention should be paid to the forward operator because the enhanced 3D radiative effects in high-resolution modeling may make the non-linearity of the forward operator more complicated. An outlook of future works is summarized below.

1. To tune up the forward operators and more accurately estimate the observation errors. The findings in this study suggest that decreasing observation errors improves the first-guess forecast and analysis. One of the factors determining the observation errors is the forward operator (Janjić et al., 2017). Therefore, it is necessary to tune up the RTTOV configurations in simulating synthetic FY-4A visible imagery from WRF model state variables and to estimate the observation errors under different modeling resolutions, weather conditions, sun-viewing geometries, etc. Scheck et al. (2018) assessed the performance of the forward operator MFASIS by comparing the synthetic visible imagery simulated on the basis of the state variables from COSMO with the SEVIRI visible image. These works could be referenced when assessing the performance of RTTOV simulations against FY-4A visible observations.

2. To improve the computational efficiency and accuracy of the forward operators. Assimilating visible radiance data is quite time-consuming for the current WRF/DART-RTTOV system (around 7 min in an EAKF cycle and 13 min in a RHF cycle). Increasing the updating frequency and outlier threshold value makes the computational cost even more expensive. We need an accurate and fast observation operator for assimilating the FY-4A (and the upcoming FY-4B) visible radiance data at both low- and high-resolution simulations. Scheck et al. (2016a) and Scheck (2021) developed a Look-up Table (LUT) and machine learning-based forward operator that is several orders faster than DOM-based methods. In addition, 3D radiative effects could be corrected for high-resolution modeling without adding expensive computation costs (Scheck et al., 2018; Albers et al., 2020; Zhou et al., 2021). These methods should be explored to improve forward operators that are suitable for assimilating the FY-4A visible radiance data.

3. To correct the errors due to the non-Gaussian and non-linear problems with the SW radiance DA. The performance of EAKF was limited in dealing with the non-linear and non-Gaussian problems. The particle filter (PF) is reported to have advantages in dealing with non-Gaussian and non-linear problems. With certain localized methods included, PF shows great potential when applied to high-dimension numerical prediction models such as the WRF model (Shen and Tang, 2015; Poterjoy, 2016; Pinheiro et al., 2019). Therefore, some newly developed PF methods could be a candidate to further improve the forecasting skills of the WRF model with the DA of satellite visible radiance.

4. To develop techniques to reduce cloud location errors. The performance of the WRF/DART system is limited in correcting the cloud location errors in the case that the first guess underpredicts clouds or a zero spread of prior ensemble members occurs. Dowell et al. (2011) proposed a method to tackle the location errors in assimilating radar data by EnKF. Their basic idea was to add some perturbations to the base state randomly and to add local perturbations in and near precipitation areas regularly so as to produce clouds in the precipitation areas. In addition, White et al. (2018) also promoted a method to produce clouds comparable to satellite observations, but this method needs to be supported by additional observations (such as brightness temperature at infrared bands). Improvements in horizontal cloud location errors could improve precipitation forecasts in some cases. For example, the vertical cloud structure is less important for some deep convection since the cloud may extend from the boundary layer up to the troposphere (X. Hu et al., 2021).

5. To evaluate the impact of assimilating FY-4A visible radiance data for long-term forecasts. The limited impact of the DA on rain rate is partly caused by the shortcomings of the DA procedure discussed above but also by the spin-up effects. The spin-up effects may introduce false-alarm precipitation due to the interactions between the model dynamics and microphysics when smaller scales are now well represented in the initial conditions and lateral boundary conditions (Short and Petch, 2022). In this study, we started to assimilate the synthetic FY-4A visible radiance data after 2 h cold start and ran the model for a 10 h forecast (02:00–12:00 UTC). This time period is too short to exclude the spin up. Improvements in precipitation should be expected for longer forecasts.

Finally, for the real FY-4A visible radiance DA, the NWP model errors in cloud and precipitation forecasts should be considered in the DA processes. On the one hand, the parameterization (such as the subgrid-scale cloud fraction parameterization in this study) is closely related to the calculation of synthetic visible radiance by a forward operator. On the other hand, the formation and dissipation of cloud and precipitation depend greatly on the model parameterization. Suboptimal parameterization may introduce large model bias in some cases due to unsolved scales and processes (Janjić et al., 2017). The model bias could introduce a negative impact on cloud and precipitation forecasts. Therefore, the NWP model should be tuned to properly represent the scale-dependent microphysical processes in order to fully realize the effects of the FY-4A visible radiance DA.