Geostationary infrared satellite observations are spatially dense [>1/(20 km)2] and temporally frequent (>1 h−1). These suggest the possibility of using these observations to constrain subsynoptic features over data-sparse regions, such as tropical oceans. In this study, the potential impacts of assimilating water vapor channel brightness temperature (WV-BT) observations from the geostationary Meteorological Satellite 7 (Meteosat-7) on tropical convection analysis and prediction were systematically examined through a series of ensemble data assimilation experiments. WV-BT observations were assimilated hourly into convection-permitting ensembles using Penn State’s ensemble square root filter (EnSRF). Comparisons against the independently observed Meteosat-7 window channel brightness temperature (Window-BT) show that the assimilation of WV-BT generally improved the intensities and locations of large-scale cloud patterns at spatial scales larger than 100 km. However, comparisons against independent soundings indicate that the EnSRF analysis produced a much stronger dry bias than the no data assimilation experiment. This strong dry bias is associated with the use of the simulated WV-BT from the prior mean during the EnSRF analysis step. A stochastic variant of the ensemble Kalman filter (NoMeanSF) is proposed. The NoMeanSF algorithm was able to assimilate the WV-BT without causing such a strong dry bias and the quality of the analyses’ horizontal cloud pattern is similar to EnSRF’s analyses. Finally, deterministic forecasts initiated from the NoMeanSF analyses possess better horizontal cloud patterns above 500 km than those of the EnSRF. These results suggest that it might be better to assimilate all-sky WV-BT through the NoMeanSF algorithm than the EnSRF algorithm.
Forecast errors and uncertainties are generally sensitive to the errors and uncertainties in the initial conditions (Ehrendorfer and Errico 1995; Kleeman and Majda 2005; Zhang et al. 2003). Forecast errors refer to the literal difference between a forecast and true state of the atmosphere, and uncertainty refers to the imprecision of the forecast (e.g., the spread of an ensemble forecast). To constrain these initial errors and uncertainties, results from previous forecast cycles (also known as the first guess) and observations are combined through data assimilation (DA). The ability of DA to constrain initial errors and uncertainties depends on, among other things, the frequency and density of the assimilated observations (Benjamin et al. 2004; Liu and Rabier 2002; Ying and Zhang 2018; Ying et al. 2018).
Over the tropical oceans, in situ observations of the atmosphere, such as surface stations and radiosonde soundings, are often infrequent and sparse (Ingleby 2017). It is thus a challenge to constrain subsynoptic errors over the tropical oceans using in situ observations alone (Ying et al. 2018). In contrast, spatially dense remote infrared observations are frequently available from geostationary satellites. These high spatiotemporal resolution observations might have the potential to improve the constraining of subsynoptic features over tropical oceans.
Early direct uses of geostationary infrared satellite radiance observations in DA were mostly limited to clear-sky radiance observations. The introduction of these observations has improved analyzed temperature and humidity fields (Köpken et al. 2004; Munro et al. 2004; Yang et al. 2017; Burrows 2018). Nonetheless, about 75% of the globe is covered in clouds (Wylie et al. 2005), meaning that a large portion of infrared observations were not utilized. Some of the difficulties encountered in using cloud-affected geostationary infrared radiances are as follows. First, the nonlinear relationships between radiances, cloud fractions and cloud-top heights can cause variational DA minimization failure. (Bauer et al. 2011; McNally 2009). Second, when it comes to ensemble Kalman filters (EnKF), at pixels where the all ensemble members are cloud free, but the observations indicate cloudy skies, the EnKF cannot generate a cloud since the ensemble-derived covariances between simulated radiances and cloud parameters are zero.
The assimilation of overcast infrared observations (observations where cloud fraction is unity) is currently being explored in various global weather centers. This inclusion increased the amount of assimilated data by about 5%. The inclusion of overcast observations have improved the analyses of 500 hPa geopotential heights, mean sea level pressure, 850 and 250 hPa winds by up to 2% globally (Bauer et al. 2011). The European Centre for Medium-Range Weather Forecasts now operationally assimilates overcast infrared radiances and includes cloud parameters into its analyzed state at locations of overcast infrared radiance observations (ECMWF 2018).
In recent years, several studies have been conducted to estimate the impacts of assimilating infrared satellite observations under all-sky conditions (both clear- and cloudy-skies) on numerical weather prediction. Using idealized observing system simulation experiments (OSSEs), Otkin (2010, 2012) demonstrated that the assimilation of cloudy-sky infrared observations can potentially improve cloud and moisture analyses of midlatitude high impact weather systems. The potential in assimilating all-sky infrared satellite radiances into a convection-permitting ensemble was demonstrated by Zhang et al. (2016). Their OSSEs of Hurricane Karl (2010) showed that the assimilation of all-sky infrared satellite radiances can potentially improve the thermodynamic analyses of the hurricane. Subsequent real data experiments of tropical cyclones produced results that are consistent with their findings (Honda et al. 2018; Minamide and Zhang 2018; Zhang et al. 2019). Besides tropical cyclones, the assimilation of all-sky radiances also has the potential to improve the prediction of tornadic thunderstorms in North America (Zhang et al. 2018).
Compared to tropical cyclones and midlatitude severe weather systems, very little work has been done to estimate the potential impacts of assimilating all-sky infrared observations on the analysis and prediction of typical tropical convection. While typical tropical convection often appears to be organized over large scales (Wheeler and Kiladis 1999; Kiladis et al. 2009; Zhang 2005; Matsuno 1966), it manifests in a noisier spatial pattern than that of tropical cyclones and midlatitude severe weather. Given that, it is of interest to determine if the improvements seen in cases of tropical cyclone and midlatitude severe weather can be replicated on typical tropical convection.
Ying and Zhang (2018) examined the potential impact of assimilating various satellite-based observations into a gray-zone resolution ensemble through OSSEs. They demonstrated that the assimilation of all-sky infrared radiance observations can reduce errors in cloud hydrometeors during the active phase of the October 2011 Madden–Julian oscillation (MJO; Madden and Julian 1971, 1972). This improvement to hydrometeors is consistent with results from earlier case studies of tropical cyclone and midlatitude severe weather.
This study extends the work done by Ying and Zhang (2018) from OSSEs to the assimilation of real infrared radiance observations through the ensemble Kalman filter (EnKF). The EnKF experiments cover a 3-day period from 16 to 19 October 2011, which was during the active phase of the October MJO event observed by the 2011 Dynamics of the Madden–Julian Oscillation (DYNAMO) field campaign (Zhang et al. 2013). In this paper, we will focus on the impact of assimilating all-sky satellite infrared radiance observations on cloud patterns and moisture analyses. A stochastic EnKF scheme will be proposed to mitigate the erroneous drying impact introduced during the EnKF analysis step.
The methods utilized in this study will be described in section 2. The impacts of assimilating infrared radiance observations on horizontal cloud patterns across different horizontal scales are explored in section 3. These impacts will be broken down by intensity and position of the patterns. Section 3 will also show that assimilating infrared radiance observations can generate dry biases, discuss the mechanism that generated the bias, and, propose and evaluate a new stochastic EnKF scheme to mitigate this issue. The evolution of errors in the simulated infrared radiance fields of deterministic forecasts initialized from the analysis means will also be explored in section 3. Conclusions are provided in section 4.
2. Experimental design
a. Forecast model, domain of interest, and ensemble setup
The Weather Research and Forecasting (WRF; Skamarock et al. 2008) Model, version 3.8.1, was utilized in this study with a configuration similar to that of Ying and Zhang (2018). The model domain is set up over the equatorial Indian Ocean (Fig. 1), with 432 × 243 horizontal grid points at 9-km grid spacing. There are 45 model levels with the bottommost 9 levels within the lowest 1 km. The top of the domain is set to 20 hPa.
To simulate radiative processes, the updated Goddard shortwave scheme (Chou and Suarez 1999) and the Global Circulation Model version of the Rapid Radiative Transfer Model (RRTMG) longwave radiation scheme were used. The unified Noah land surface physics scheme (Chen and Dudhia 2001) was used to simulate surface processes, except that the surface skin temperature is diagnosed separately using the scheme of Zeng and Beljaars (2005). Subgrid-scale turbulent mixing is handled by the Yonsei University (YSU) boundary layer scheme (Hong et al. 2006). The WRF double-moment 6-class scheme (Lim and Hong 2010) was used to represent cloud microphysical processes.
No cumulus parameterizations were employed in this study, as the 9-km horizontal grid spacing (also called gray-zone resolution) used in this study is fine enough to resolve mesoscale convective complexes, although it’s still too coarse to resolve convective-scale updrafts. Wang et al. (2015) compared similar gray-zone simulations against the DYNAMO observations. Their results show that the gray-zone resolution simulation was able to realistically replicate the general precipitation, circulation, thermodynamic, and radiation features of the October and November MJO events observed during DYNAMO. Subsequent work by Chen et al. (2018b) showed that the gray-zone resolution simulation can capture the general features of atmospheric overturning across different scales well. Similar model setups have also been used to study the physical mechanisms of tropical convection (Zhang et al. 2017; Chen et al. 2018a,c; Chen and Zhang 2019).
The WRF ensemble used in this study was initialized from the European Centre for Medium-Range Weather Forecasts (ECMWF) 50-member perturbed forecast system at 0000 UTC 15 October 2011. These perturbed forecasts were archived in The Observing System Research and Predictability Experiment (THORPEX) Interactive Grand Global Ensemble (TIGGE; Swinbank et al. 2016). Missing data above 200 hPa, sea surface temperatures, orography, and land sea mask were filled in using the ECMWF interim reanalysis (ERA-I) dataset (Dee et al. 2011). The WRF ensemble was integrated forward for 24 h, without DA, to generate flow-dependent ensemble statistics. All EnKF experiments in this study start at 0000 UTC 16 October 2011 and end at 0000 UTC 19 October 2011.
b. Observation systems
This study utilized brightness temperature observations from the infrared water vapor channel (6.4 μm central wavelength; WV-BT) and infrared window channel (11.5 μm central wavelength; Window-BT) of the Meteosat Visible and Infrared Imager (MVIRI) on board the Meteorological Satellite 7 (Meteosat-7). Meteosat-7 is a geostationary satellite situated over the equator at 58°E, from December 2006 to April 2017, with a field of view that includes this study’s model domain. Half-hourly full-disk scans are available online from the Earth Observing Laboratory portal maintained by the University Corporation for Atmospheric Research/National Center for Atmospheric Research (UCAR/NCAR–Earth Observing Laboratory 2012a,b,c). Over the study’s domain, Meteosat-7’s infrared channels have a horizontal resolution of 5 km. The WV-BT observations were assimilated in this study.
Aside from the infrared satellite observations, soundings from the Dynamics of the MJO/Cooperative Indian Ocean Experiment on Intraseasonal Variability (DYNAMO) field campaign were used for validation (Johnson and Ciesielski 2013). Here, we will use soundings from the Northern Sounding Array and Southern Sounding Array (green circles in Fig. 1), every 6 h.
The rainfall produced by deterministic forecasts initiated from the DA experiments were validated against Tropical Rainfall Measuring Mission (TRMM)-derived precipitation rates (3B42; Tropical Rainfall Measuring Mission 2011). This dataset was downloaded from the National Aeronautics and Space Administration’s (NASA’s) Earth Science Data System archive.
c. DA methods and WV-BT assimilation specifications
The DA system used in this study is the Pennsylvania State University WRF EnKF system (PSU-EnKF; Meng and Zhang 2008). In the EnKF, an ensemble of model realizations is used to estimate flow-dependent error covariances between simulated observable quantities and model states (Evensen 1994). These covariances are used to pass information from observations into the ensemble of model states.
The PSU-EnKF system first ingests the input ensemble (prior ensemble). The prior ensemble is separated into an ensemble mean and an ensemble of perturbations from the ensemble mean state. The ensemble mean and ensemble perturbations are updated separately in the PSU-EnKF. The resulting mean state and ensemble of perturbations are then combined to produce new ensemble analyses of model states (posterior ensemble). The square root ensemble filter algorithm from Whitaker and Hamill (2002) was used to update the ensemble of perturbations.
The PSU-EnKF system is a serial filter, meaning that observations are assimilated one at a time. The procedure is as follows. Supposing there are M observations, the PSU-EnKF assimilates the first observation into the prior ensemble to produce the first updated ensemble. Then, the PSU-EnKF assimilates the second observation into the first updated ensemble to produce the second updated ensemble. This process repeats until all M observations have been assimilated, and the Mth updated ensemble is the posterior ensemble. For an ensemble of N members, this recursive process for the conventional version of the PSU-EnKF (EnSRF) can be mathematically expressed as
where m = 1, 2, …, M, and, n indicates the ensemble member being updated (n = 1, 2, …, N). and are the mean and nth perturbation joint state-observation vectors (henceforth, extended state vector) after assimilating the mth observation (i.e., mth updated ensemble), and likewise for and . These extended state vectors contain both the model state and simulated observations. Note that the mean and perturbation extended state vectors that are used to start the recursive process (i.e., m = 0, or the prior ensemble) are
where is the nth prior ensemble member’s model state and h is an operator that generates simulated observations from a model state. It should be noted that h is likely to be a nonlinear operator when it comes to satellite radiance observations; is a linear operator that extracts the mth simulated observation from an extended state vector; and and ϕm are, respectively, the Kalman gain matrix (a vector in this case) and square root filter factor. For the serial ensemble square root filter, they are defined as
where is the ensemble variance of the simulated mth observation, after assimilating m − 1 observations, and is the mth observation error variance. To reduce the negative impacts of clear/cloudy-air representativeness errors, the adaptive observation error inflation scheme (AOEI; Minamide and Zhang 2017) was applied in the conventional form of the PSU-EnKF.
It will later be shown and explained in section 3 that the conventional PSU-EnKF introduces an unphysical water vapor removal from the domain at each update step. To mitigate that, a stochastic filter version of the PSU-EnKF (NoMeanSF) is proposed. Unlike the conventional PSU-EnKF, the NoMeanSF updates each ensemble member directly, without the separation into ensemble mean and ensemble perturbations. The recursive update equation of NoMeanSF is
where and ϵm,n is a random number drawn from a monovariate normal distribution with mean zero and variance . This choice of is very similar to the AOEI, except that varies across ensemble members, while the AOEI uses the same for all ensemble members.
The reason for this ensemble-varying is to prevent potentially unphysical updates to members that disagree with observations by dramatically inflating the observation errors for these members. Note that this ensemble-varying causes the ensemble-averaged impact of including observation error perturbations ϵm,n to be nonzero. This effect might potentially degrade the quality of analysis and introduce a bias. Furthermore, the use of adaptive observation error inflation methods can potentially violate an important assumption of data assimilation; that is, observation and background errors are independent (Geer et al. 2018). Nonetheless, the potential introduction of a bias or violation of an assumption is preferable over unphysical updates in the model, which may lead to numerical instability. Hence, we have opted to use this ensemble-varying .
The NoMeanSF is a variant of the perturbed observation ensemble Kalman filter (Evensen 2003). The main difference is that the NoMeanSF uses an adaptive observation error inflation scheme that rescales the Kalman gain matrix differently for each ensemble member. If this scheme was turned off and assuming perfect Gaussian statistics, the posterior ensemble mean and covariances from the NoMeanSF should be identical to that of the perturbed observation ensemble Kalman filter and an EnSRF that uses instead of .
For both the EnSRF and NoMeanSF, the adaptive background error inflation scheme (ABEI; Minamide and Zhang 2019) was applied. Note that the regression coefficients used in ABEI were taken directly from Minamide and Zhang (2019) in this study. Further improvement might be expected if the regression coefficients are specifically calculated using an OSSE over the Indian Ocean domain, which deserves future exploration. To prevent filter divergence, 80% relaxation to prior perturbation (Zhang et al. 2004) was applied. The Gaspari and Cohn (1999) fifth-order polynomial tapering function was utilized to localize covariances to mitigate sampling noises in the ensemble error covariance.
d. Experiments conducted
Table 1 shows a list of experiments that are examined in this piece. The no DA experiment (NoDA) is a free ensemble run that is used to gauge the upper limit of errors. Experiments using the conventional PSU-EnKF system (EnSRF) and the stochastic filter experiment (NoMeanSF) to assimilate WV-BT were also performed.
To assimilate the infrared observations, the Community Radiative Transfer Model (CRTM), release 2.1.3, was used to calculate the infrared brightness temperatures (the h operator). Both vertical and horizontal localization were used in assimilating WV-BT. A vertical radius of influence (VROI) of 45 model layers (equivalent to the vertical extent of the model domain) was used and the center of the vertical localization depends on whether the satellite-observed pixel was under clear (400 hPa) or cloudy (250 hPa) sky conditions. Similar to Ying and Zhang (2018), Window-BT observations were first used to determine if an observed satellite pixel is clear (Window-BT ≥ 285 K) or cloudy (Window-BT < 285 K). The WV-BT observations were thinned to 10 km grid spacing, and assimilated hourly with a 100 km horizontal radius of influence (HROI). Several different cycling periods (1 and 3 h), observation thinning distances (10 and 30 km) and HROI’s (300 and 100 km) were also tested (not shown). The configuration that yielded the best fit to observed infrared radiances utilized hourly EnKF cycling with a thinning distance of 10 km and HROI of 100 km. This configuration is used in this paper. No bias correction was performed and none of the WV-BT observations were rejected. The variables updated by the Kalman filter include the three-dimensional winds, all hydrometeor mixing ratios, water vapor mixing ratios, temperature and pressure. Better results might be possible with more tuning, which deserves to be explored in future studies.
3. Results and discussion
a. Assimilation of WV-BT in EnSRF
Results from other studies using the conventional form (EnSRF) of the PSU-EnKF (Zhang et al. 2016; Minamide and Zhang 2018; Zhang et al. 2018) indicate that the most immediate impact of assimilating infrared satellite radiances is the improvement of cloud fields. This impact is also seen in the EnSRF experiment. Figure 2 shows the typical impact of assimilating WV-BT at 0600 UTC 16 October (the seventh DA cycle), in the EnSRF experiment. As can be seen from Fig. 2e, the prior mean has excessively large patches of cold Window-BT (blue patches) and generally warmer Window-BT minima than observed. These suggest excessively large clouds and lower-than-observed cloud tops. In fact, the NoDA experiment exhibits similar disagreements with the observations, albeit more strongly (not shown). The assimilation of WV-BT partially mitigated these issues, as can be seen in Fig. 2g.
The impact of assimilating WV-BT on cloud features of the prior ensemble over time was also examined. Figures 3a and 3e shows that the assimilation of WV-BT in the EnSRF was able to reduce the root-mean-square difference (RMSD) between the observed radiance fields and the radiance fields of the prior mean. Toward the end of the 3-day period, Figs. 3a and 3e indicates the EnSRF has less than 60% of NoDA’s RMSD in both WV-BT and Window-BT. Furthermore, while the assimilation process reduced the spread of the ensemble, as measured by the root-mean-square spread (RMSS; which also includes the observation error), the ensemble’s RMSS and RMSD values were similar throughout the experiment. This similarity indicates that the ensemble maintained an appropriate amount of dispersion with respect to cloud features on the overall. In other words, the assimilation of WV-BT in the EnSRF improved horizontal cloud patterns throughout the study’s time frame and constrained the uncertainty in cloud features appropriately.
It should be pointed out here that the EnSRF was able to maintain appropriate amounts of RMSS relative to RMSD, despite the NoDA experiment showing signs of underdispersion. The relaxation to prior perturbations (RTPP) applied in the EnKF is likely responsible for maintaining the spread in the EnSRF. Mixing the posterior ensemble perturbations with the prior ensemble perturbations generally result in perturbations larger than the posterior ensemble’s. Since model integration typically increases spread, the resulting prior perturbations are usually larger than if RTPP was not applied, helping to prevent underdispersion.
The implied improvement in horizontal cloud patterns can be broken down by ranges of spatial scales. This decomposition can be achieved by bandpass filtering the radiance fields in 2D Fourier wavenumber space. Three spatial scale ranges will be examined here: L-scale (features larger than 500 km), M-scale (features between 100 and 500 km), and S-scale (features smaller than 100 km). It can be shown that the decomposition is linear in terms of mean-squared differences:
According to Fig. 3e, after 24 h of DA, the mean-squared difference (MSD) of EnSRF’s Window-BT (WV-BT) is roughly 300 K2 (36 K2), which is 600 K2 (64 K2) less than the roughly 900 K2 (100 K2) MSD seen for NoDA. However, Fig. 3f indicates that the S-scale EnSRF Window-BT (WV-BT) MSD is slightly worse (only slightly better) than that of NoDA. Combining these facts with Eq. (4) reveal that most of the improvements in Window-BT (WV-BT) come from scales above 100 km. Furthermore, the improvement to the L-scale MSD seen in EnSRF’s Window-BT (WV-BT) is about 500 K2 (60 K2), which accounts for majority of the improvement seen in the overall RMSD. The assimilation of WV-BT also resulted in the three spatial scales having similar contributions to overall MSD (about 30% each). The forecast improvement at the L-scale and M-scale could indicate that errors at these spatial scales grow linearly on the time scale of the DA update frequency. Also, the RMSS at M- and L-scales are close to the corresponding RMSD values, meaning that the ensemble maintained an appropriate amount of dispersion with respect to horizontal cloud patterns at these scales. The maintenance of ensemble dispersion is likely due to the use of RTPP in this study, as explained earlier.
In contrast to the M- and L-scales, the assimilation of WV-BT only improved WV-BT’s S-scale RMSD by 1 K and slightly worsened Window-BT’s S-scale RMSD. There are two possible reasons for this weak update. First, the HROI in this study is set to 100 km, meaning that the analysis increment is smoothed with a similar length scale, possibly reducing corrections to sub-100 km scales. Another possibility is that the model’s effective resolution, which is approximately 4-6 times the 9 km horizontal grid spacing (Skamarock et al. 2014), might be too coarse to resolve sub-100 km features. In other words, the prior ensemble covariance and mean state might not realistically represent the sub-100 km features, potentially causing erroneous and/or weak updates at these scales. These two possibilities are not mutually exclusive.
The biases in the ensemble means’ Window-BT and WV-BT were also reduced by assimilating WV-BT observations (Fig. 4). Bias here is defined as the differences between the spatial average of prior mean’s BT field and observed BT field. Toward the end of the 3-day period, EnSRF’s prior Window-BT and WV-BT cold biases are more than 50% smaller than NoDA. These improvements imply that the assimilation of WV-BT helped to reduce the excessive cloud cover seen in NoDA.
To further understand how the EnSRF improved horizontal cloud patterns, the power spectra of radiance fields in terms of horizontal Fourier wavenumbers was examined. The spectral power corresponding to horizontal wavenumber kH,j is computed by first applying 2D fast Fourier transforms on the radiance fields, and then applying the following definition:
where ψ is a placeholder for BT fields or BT difference fields (e.g., WV-BT of NoDA ensemble mean minus observed WV-BT), is the Fourier coefficient corresponding to zonal and meridional wavenumbers k and l, and Sj is a set of (k, l) wavenumber pairs produced by the 2D-FFT such that . The horizontal power spectra of Window-BT and WV-BT reveal the intensity of cloud pattern variations at different length scales.
The horizontal spectral power of BT difference fields ΔBT is connected to overall RMSD range through
where NH is the number of horizontal wavenumbers. It should also be mentioned here that the spectral power of wavenumber 0 is equal to the square of the domain-average bias in Fig. 4.
Figure 5 shows various horizontal power spectra for the observed radiance fields and the simulated radiance fields from the prior ensemble means. For scales under 500 km, the horizontal power spectra of observed WV-BT and Window-BT fields have slopes that are similar to that of a −5/3 power law (e.g., Nastrom and Gage 1985). The origin of the −5/3 slope in radiance fields is not known and might be an area of future research. Nonetheless, the NoDA experiment was able to produce a sub 500 km slope that is similar to the observations. This indicates that even without DA, the intensity distribution of cloud patterns by horizontal length scales is generally correct, which reflects that the model at gray-zone resolution can capture the basic dynamics over the tropical ocean region.
A closer look at the power spectra from the NoDA experiment reveals that, in the absence of DA, features between 1000 and 2000 km are excessively intense and features between 100 and 1000 km are excessively weak (green line in Figs. 5a and 5e). The assimilation of WV-BT in EnSRF mitigated these effects and reduced the overall horizontal difference power for scales above 100 km. Furthermore, without DA, the prior ensemble spread spectra are lower than the difference power spectra for scales above 100 km for both Window-BT and WV-BT. This suggests that the ensemble is generally underdispersive for scales above 100 km, which is consistent with Figs. 3e–h. With the assimilation of WV-BT, the underdispersion is mitigated. This mitigation is consistent with those of the RMSD and RMSS time series.
The difference power spectra can be linearly decomposed into components relating to the horizontal wave pattern amplitude differences and phase differences. The phase difference term can be interpreted as the contribution of horizontal cloud pattern dislocation to the difference power spectra. To derive these terms, suppose that the Fourier coefficients of the two fields for wavenumber pair (k, l) are and . Then,
where ϕo and ϕf are the phase angles associated with and . The first term on the last line of Eq. (7) can be interpreted as the contribution of magnitude difference to the difference power, and the second term can be interpreted as the contribution of phase difference to the difference power. As a side-note, Eq. (7) reveals that if the magnitudes of and are identical, the maximum difference power is 4 times the amplitude squared of .
The difference power contributions from magnitude and phase differences were quantified through Eq. (7) and converted into horizontal wavenumber space via Eq. (5). The resulting magnitude and phase difference power spectra are shown in Fig. 5. Phase differences are generally the dominant contribution to the difference power spectra for horizontal wavelengths under 1000 km. In other words, the dominant source the disagreement between the observed and prior mean’s radiance fields, for scales under 1000 km, is the dislocation of cloud patterns.
Assimilating WV-BT mitigates both the magnitude and phase difference power spectra for scales above 100 km and the improvements generally diminish with smaller length scales. Furthermore, assimilating WV-BT generally reduced the phase difference power spectra more than that of the magnitude difference power spectra. Taken together, the assimilation of WV-BT generally improved the intensity and position of cloud patterns, with stronger improvements to cloud pattern position than cloud depth.
It is also interesting to examine the spectral performance of NoDA and EnSRF in terms of the average of prior ensemble members’ radiance fields (Fig. 6). The results are qualitatively similar to that of Fig. 5: EnSRF’s ensemble-averaged power spectra are generally closer to the observed power spectra than that of NoDA. Also, the dominant source the disagreement between the observed and prior ensemble’s radiance fields, for scales smaller than 1000 km, is the dislocation of cloud patterns.
b. Validation of EnSRF against sounding observations
Aside from comparisons against the infrared radiance fields, the experiments were also compared against the sounding observations obtained from the DYNAMO field campaign. The 6 sounding sites used in this paper are as indicated in Fig. 1. Without DA, the ensemble mean is generally dryer and cooler than observed (Figs. 7b,d). Furthermore, the average observed westward component wind is generally stronger (weaker) than NoDA’s at levels under (above) 200 hPa (Fig. 7f). When it comes to the meridional component of the wind, NoDA exhibits a more southward flow between 600 and 300 hPa, and above 200 hPa (Fig. 7h).
The assimilation of WV-BT in the EnSRF worsened the dry bias seen in NoDA in the mid- to lower troposphere. The worst dry bias is at 850 hPa, where EnSRF’s dry bias is roughly three times that of NoDA. In terms of temperature, the EnSRF is slightly colder than NoDA in general. The 1000 hPa temperature of the EnSRF, NoDA, and ECMWF-PF experiments are much colder than that of the observations (~3 K for EnSRF and NoDA, ~1.75 K for ECMWF-PF). It is not clear if this large difference is induced by the sounding observation errors near surface, or because the model has issues resolving near-surface temperatures. The exact reasons behind these temperature biases are beyond the scope of this paper and deserve further study.
When it comes to wind variables, EnSRF shows a similar performance with NoDA. In other words, it can be difficult to constrain wind variables using satellite radiance in this study’s case. Future work will assimilate conventional observations (e.g., atmospheric motion vectors and other soundings) to check if these variables can be better constrained. Among the four sounding variables, the impact of assimilating WV-BT through EnSRF on moisture is the most noticeable.
The fact that EnSRF’s moisture profile is much dryer than NoDA’s is contrary to expectations from earlier work (e.g., Minamide and Zhang 2018). To understand why, the domain-averaged prior and posterior ensemble mean water vapor mixing ratios were plotted as a time series (Fig. 8). Figure 8 shows that the domain-averaged moisture consistently falls during the update steps of the EnSRF (green lines in Fig. 8). On the other hand, the forecast steps generally moisten the domain, albeit insufficiently to counter the drying effects of the update steps. In other words, EnSRF’s update steps are responsible for unexpectedly the dry profile.
The analysis increments show that most of the moisture removals in cloudy pixels are associated with positive WV-BT analysis increments (not shown). This association implies negative correlations between water vapor mixing ratio and WV-BT under cloudy conditions. Since the WV-BT of EnSRF’s prior ensemble means have consistently more cloud cover than the observations (Fig. 4), the negative correlations would cause moisture removal, in a domain-averaged sense, during update steps.
Understanding the source of the excessive cloud cover can lead to ways to mitigate the dry bias. Model integration generally causes cloud dislocation among ensemble members (not shown). Considering that most water clouds are optically thick in the infrared regime, the WV-BT of any ensemble mean’s grid box approximately reflects the blackbody radiation from the tallest water cloud top among the ensemble members, as well as the temperature and water vapor content above the said cloud top. Figure 9 illustrates this line of reasoning conceptually. Thus, the excessive cloud cover seen in the WV-BT of the prior ensemble mean is due to the dislocation of clouds among ensemble members during integration. Even though the assimilation process does reduce the dislocation of clouds, model integration regenerates the dislocation effects. Following the reasoning in the previous paragraph, the dry bias introduced by EnSRF’s update step might be because the WV-BT of the prior mean was used to calculate the analysis increment.
To test the claim discussed in the previous paragraph, a different formulation of the ensemble Kalman filter was employed. Instead of updating the mean state using the WV-BT of the ensemble mean, the individual members are updated directly. This direct update utilizes the difference between the simulated WV-BT field of individual members and the observed WV-BT field, as indicated in Eq. (3) (NoMeanSF experiment). Compared to the EnSRF, the NoMeanSF’s update steps generally show much smaller impacts on the domain-averaged ensemble mean water vapor mixing ratio (Fig. 8). In fact, at the model level with an average pressure of 640 hPa, NoMeanSF’s update steps caused the domain-averaged ensemble mean water vapor mixing ratio to exceed that of NoDA. When compared against the DYNAMO soundings, NoMeanSF shows a significantly smaller dry bias than EnSRF, though with a slightly larger dry bias than that of NoDA (Fig. 7b). This reduced bias indicates that the drying introduced by EnSRF’s update step is likely due to the use of the WV-BT of the prior mean to calculate the analysis increment.
When compared to the observed radiance fields, while the performance of the radiance fields of NoMeanSF’s prior ensemble mean is slightly weaker than that of EnSRF, their performances are comparable (Figs. 3–5). This comparable performance happened even though NoMeanSF does not explicitly constrain the radiance fields of the ensemble mean, whereas the EnSRF does [see Eqs. (1) and (3)]. Taking the reduced dry bias into account, it would thus seem that the NoMeanSF’s prior ensemble means outperform those of EnSRF.
d. Comparisons of the deterministic forecasts of EnSRF and NoMeanSF
The deterministic forecasts initiated from the posterior mean of the EnSRF and NoMeanSF are compared here. Figure 10 shows the forecast-cycle-averaged radiance field RMSD’s of EnSRF and NoMeanSF’s deterministic forecasts. These deterministic forecasts are initiated every 6 h and run for 24 h. Similar plots for the radiance field biases of the two experiments’ deterministic forecasts are in Fig. 11. A total of 8 forecast cycles (between 0600 UTC 17 October and 0000 UTC 19 October) were used to produce Figs. 10 and 11.
While the overall and L-scale WV-BT RMSD performances of the deterministic forecasts of NoMeanSF and EnSRF are generally similar, there are some large differences when it comes to the Window-BT channel. NoMeanSF’s deterministic forecasts have much smaller Window-BT L-scale RMSDs than that of EnSRF (except at initialization). Furthermore, the biases of the radiance fields of NoMeanSF’s deterministic forecasts are generally better than that of EnSRF (Fig. 11). The difference between the squared biases of Window-BT for the EnSRF’s deterministic forecasts (roughly 300 K2) and NoMeanSF’s deterministic forecasts (roughly 200 K2) accounts for most of the L-scale deterministic forecast MSD differences (Fig. 10h). It can thus be said that NoMeanSF generally produces better cloud pattern deterministic forecasts than EnSRF for scales above 500 km. When it comes to the M-scale and S-scale RMSD, the deterministic forecasts of NoMeanSF do not perform as well as the EnSRF deterministic forecasts.
The distributions of cloud patterns by wavelengths in both sets of deterministic forecasts are also examined. Figure 12 shows the power spectra of the deterministic forecasts for both WV-BT and Window-BT, at different lead times. For both radiance fields and for both sets of deterministic forecasts, the power spectra first weakened, and then strengthened toward the end of the 24-h forecasts. This evolution pattern suggests that the cloud features at the start of the deterministic forecasts were not sustained, resulting in a general loss of spatial variability in the radiance fields. Insufficiently strong convective updrafts and/or less favorable thermodynamic conditions in the posterior mean might explain this lack of sustenance. It is also possible that the model needed time to spin up since the deterministic forecasts were initialized from an ensemble average and averaging tends to smooth out dislocated small-scale features. The second phase in the power spectra evolution suggests that after the first few hours, the model generated new convection that increased the spatial variability of the simulated radiance fields.
Comparing the deterministic forecasts’ power spectra is also another way to evaluate the performances of the two assimilation schemes. While NoMeanSF’s power spectra are usually weaker than that of the observations, the structures observed in NoMeanSF’s power spectra are generally similar to those observed, especially for the WV-BT channel. In contrast, the EnSRF’s deterministic forecasts’ power spectra changed from the observed shape over time. After 12 h, the proportions of large cloud patterns (>2000 km) to sub 2000 km cloud patterns are larger than that of the observations. In other words, beyond 12 h, NoMeanSF shows better cloud pattern distribution by horizontal scale.
It should also be pointed out that between 12 and 18 h lead time, the EnSRF generally has better cloud pattern intensities than the NoMeanSF at scales above 2000 km. However, at these lead time and spatial ranges, the difference power spectra of EnSRF exceed that of NoMeanSF (Fig. 12, second row). These stronger difference power spectra are due to a much stronger phase difference contribution spectra (not shown) in EnSRF than NoMeanSF, which overwhelmed the effect of the weaker amplitude performance in NoMeanSF. Thus, even these for these scales and lead times, NoMeanSF performed better than that of EnSRF.
The deterministic forecasts of NoMeanSF and EnSRF were also compared against TRMM rainfall estimates (Fig. 13). The deterministic forecasts from EnSRF and NoMeanSF both underestimate the observed rainfall intensity. A similar result is also noted for deterministic forecasts initiated from NoDA’s ensemble mean (Fig. 13b). The weaker rainfall in the deterministic forecasts might be due to the dry bias present in the initial and boundary conditions. NoMeanSF’s deterministic forecasts have the smallest rainfall bias compared to NoDA and EnSRF’s. This might be due to NoMeanSF generating a moister midtroposphere (Fig. 8a) than both NoDA and EnSRF. In addition, the deterministic forecasts of NoMeanSF better capture the westward propagation of rainfall than that of EnSRF and NoDA. Thus, NoMeanSF is better than EnSRF in terms of the rainfall forecast.
Overall, it seems that NoMeanSF generally outperforms the EnSRF in terms of the deterministic forecasts of radiance fields and rainfall rates. Furthermore, it has been found that NoMeanSF’s prior ensemble mean has a much smaller dry bias than that of EnSRF’s. Putting these together with the fact that the radiance fields of NoMeanSF’s and EnSRF’s prior ensemble means have comparable performance, it appears that NoMeanSF outperforms the EnSRF.
e. Error behavior of NoMeanSF’s deterministic forecast
The saturation times of the RMSD at various scales indicate the lifespans of the improvements introduced by DA at the corresponding scales. According to Fig. 10, the M-scale RMSD of NoMeanSF’s deterministic forecasts, for both infrared channels, are close to saturation within 6 h. Interestingly, the L-scale WV-BT RMSD appears to be close to saturation within 19 h, but that of Window-BT does not saturate within 24 h. Since Window-BT does reflect cloud tops below 800 hPa, whereas WV-BT should not, this difference in saturation times indicates that L-scale errors related to cloud and moisture above 800 hPa saturate faster than those beneath.
Figure 14 shows the evolution of the difference spectra, amplitude difference spectra and phase difference spectra of NoMeanSF’s deterministic forecasts. The total and phase difference spectra of both channels generally increase with time and saturate earlier at the smaller scales than at the larger scales. An exception to this pattern is the evolution of the Window-BT total and phase difference spectra during the first hour for scales above 1000 km—they decreased instead. The loss in phase difference spectra is partially due to the loss of spectral power above 1000 km in the lower troposphere during the first hour (Fig. 12c). Since the total difference spectra are dominated by the phase difference spectra at initialization, this reason might also explain the loss of total difference spectra above 1000 km at the first hour in the Window-BT channel.
The amplitude difference spectra for wavelengths above 1000 km generally rise for the first 16 h, and then fall. The increase in amplitude difference spectra are associated with the loss of power spectra at these scales for the first 16 h seen in Fig. 12. As discussed earlier, this loss in power spectra might be because the dynamic/thermodynamic fields in the initial conditions could not sustain the clouds, and/or the model requiring some time to spin up due to the loss of small-scale features during averaging. The later amplitude difference spectra reduction phase is associated with the power spectra in Fig. 12c rising toward the observed spectra. This might be the result of the model having spun up after 16 h, and/or the development of circulations and thermodynamic environments that can support clouds. Finally, the evolution of the amplitude difference power spectra at the smaller scales tends to slow down earlier than that of the larger scales.
Despite the nonmonotonic evolution seen in the amplitude difference spectra above 1000 km, the difference spectra evolve in a generally monotonic fashion. The reason is that like the various difference power spectra of the prior ensemble means (Fig. 5), phase differences generally contribute to the majority of the difference power spectra in the deterministic forecast. Furthermore, at lead times after 12 h, that contribution is responsible for most of the difference power spectra and the phase difference power spectra grows much faster than the amplitude difference power spectra. As such, the nonmonotonic evolution of the amplitude difference power spectra is overwhelmed by the monotonic evolution of the phase difference power spectra. This would explain the monotonic evolution of the difference power spectra.
4. Summary and future work
This study examined the impact of assimilating geostationary infrared satellite radiance observations, every hour, into a gray-zone resolution ensemble over the equatorial Indian Ocean, outside of a tropical cyclone. The resulting short-term ensemble mean forecasts (i.e., priors) show improved matches with the observed WV-BT and Window-BT, with most of the improvement coming from scales above 500 km. For scales under 100 km, DA slightly improved (slightly worsened) the match to the WV-BT (Window-BT) observations. Subsequent examination of the horizontal power spectra of the ensemble mean and ensemble members revealed that DA also improved the general distribution of features by length scales. It was also found that much of the difference between the ensemble means/members and the observed infrared fields by scale are due to mismatches in horizontal cloud pattern positions. DA constrained the cloud position and intensity mismatches in a generally proportionate sense.
When it comes to dynamical variables, the use of the WV-BT of the ensemble mean in the EnKF algorithm (EnSRF) appeared to consistently introduce dry biases into the model domain. Replacing the algorithm with a stochastic filter that did not utilize the WV-BT of the ensemble mean to calculate the mean increment (NoMeanSF) mitigated much of this drying effect. Furthermore, the simulated infrared brightness temperature fields of NoMeanSF’s and EnSRF’s prior ensemble means had similar performances. The simulated infrared brightness temperature fields of NoMeanSF’s deterministic forecasts also showed better root-mean-square performance at scales above 500 km, compared to that of the EnSRF. Furthermore, NoMeanSF’s deterministic forecasts show better intensity distribution by wavelength than that of the EnSRF’s. However, the EnSRF’s deterministic forecasts have better sub 100 km matches to observations than that of NoMeanSF. Considering the performance of the NoMeanSF above 100 km and that it did not introduce a strong dry bias, the NoMeanSF might be a better algorithm to assimilate WV-BT in this region.
The error growth behaviors of the simulated infrared fields of NoMeanSF’s deterministic forecasts were also examined. The 100 to 500 km errors appear to saturate within 6 h of time integration for both the WV-BT and Window-BT. For errors above 500 km, the errors in the WV-BT saturated within 24 h, whereas that of the Window-BT did not saturate within the period of the deterministic forecast (24 h). It was also found that the generally monotonic error growth across horizontal scales is because of the monotonic growth of the phase difference contribution, which is the dominant contribution.
There are some caveats that should be noted about this study. First, this study assimilated satellite infrared radiances only. In other words, this study does not quantitatively indicate the benefit of assimilating all-sky infrared radiances on top of conventional observations. Future work will assimilate conventional observations as well to estimate that benefit. Nonetheless, this study might indicate that assimilating all-sky infrared radiances on top of conventional observations should generally improve the horizontal cloud patterns. Second, given the domain size of this study, it is difficult to validate the dynamical fields through soundings at any given instant. As such, only the 2-day average field was validated. Subsequent work might employ a larger domain or utilize unassimilated retrievals for validation. Last, the impact of the all-sky satellite DA on interior cloud structures and precipitation has not been explored. This can be explored in the future using the CINDY/DYNAMO radar observations.
It should also be pointed out that the NoDA experiment in this study was dry biased. The dryness of the NoDA experiment could be caused by dry biases in the ECMWF-PF data, its low vertical resolution in the TIGGE archive (only 9 pressure levels), and/or dry biases of the WRF Model itself. We found that the dry biases are reduced when the ERA-Interim reanalysis is used to provide boundary conditions (not shown here), suggesting that the ECMWF-DF data at least partly contribute to the dry biases. Aside from that, while the NoMeanSF had a better dry bias than the EnSRF, the NoMeanSF still introduced an additional dry bias. This additional dry bias could be related to sampling errors, errors in the CRTM, or even shifts in the ensemble mean state introduced by the AOEI in NoMeanSF. Future investigation is needed to unravel the source of these dry biases.
It should also be pointed out that if the adaptive observation error inflation methods are turned off, the NoMeanSF formulation is mathematically equivalent to using the EnSRF algorithm with the mean simulated brightness temperature defined as the mean of the ensemble of brightness temperature fields. This alternative form of the EnSRF can be investigated in future work.
Aside from assimilating conventional observations, future work can also investigate the impact of combining all-sky microwave observations with all-sky infrared observations. While all-sky microwave observations are used in operational DA, the impacts of combining both types of all-sky radiances are not well explored in the context of a convection-permitting ensemble over the tropical ocean.
This algorithm can potentially be used to improve the operational forecasts of cloud systems over the tropics, where in situ observations are sparse and infrequent. Also, the algorithm can be used to generate high resolution reanalysis data for tropical convection studies. Future work is needed to test this algorithm on more cases. It is also worthwhile to explore assimilation of all-sky infrared radiances together with routine observations and microwave radiances in future studies.
It is interesting to note that other studies of all-sky infrared DA using the EnSRF formulation in the PSU-EnKF system did not result in an obvious, strong dry bias found in this study (e.g., Zhang et al. 2016; Minamide and Zhang 2018). One reason for the difference might be due to the fact that previous studies often focus on tropical cyclones. In those studies, the central pressures of tropical cyclones are typically assimilated, limiting the dislocations of tropical cyclones across ensemble members. Together with the fact that convection in tropical cyclones is much more organized than tropical convection, clouds are less likely to be dislocated among ensemble members for tropical cyclones, compared to tropical convection.
This work is a first attempt at assimilating actual all-sky infrared radiance into a gray-zone ensemble over a tropical ocean, outside of tropical cyclones. The resulting impact on horizontal cloud impacts by scale are encouraging and this study was able to largely mitigate the undesirable dry bias introduced by the EnSRF via the NoMeanSF. Considering the impact of assimilating all-sky infrared radiances in tropical cyclone cases, it seems conceivable that all-sky infrared radiances might improve analyses and operational forecasts over the tropical oceans in the near future.
This research is supported by the U.S. Department of Energy Office of Science Biological and Environmental Research as part of the Regional and Global Climate Modeling program, ONR Grant N000141812517, the Graduate School Fellowship from The Pennsylvania State University Graduate School, as well as the National Science Foundation (Award 1712290). The Pacific Northwest National Laboratory is operated for the Department of Energy by Battelle Memorial Institute under Contract DE-AC05-76RL01830. The Meteosat-7 data used in this study were obtained from https://data.eol.ucar.edu/dataset/347.025, https://data.eol.ucar.edu/dataset/347.027, and https://data.eol.ucar.edu/dataset/347.029 on 1 May 2018. The sounding array data were taken from https://data.eol.ucar.edu/dataset/347.145 on 1 October 2018. The initial and boundary conditions used to set up the WRF ensemble were created from the ERA Interim (https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/) and ECMWF perturbed forecasts (https://apps.ecmwf.int/datasets/data/tigge/levtype=sfc/type=pf/), which were accessed on 1 May 2018. The outputs of the PSU-EnKF system and assimilated observation data files are archived and accessible at the Penn State Data Commons (https://doi.org/10.26208/9wg1-g948).