HTML
-
In Z13, an experimental GSI-based EnKF DA system was established and tested at 40-km grid spacing with 40 ensemble members and three-hourly assimilation cycles. The system was intended to be a prototypical system for the operational RAP (Benjamin et al., 2016), and the forecast model was configured based on the operational RAP except for the assimilation interval and resolution. The operational RAP runs hourly and has an approximate 13-km grid spacing. The choice for the current EnKF tests was dictated by the availability of computational resources required by a large number of sensitivity experiments with continuous cycles. Another consideration is the potential use of a lower-resolution EnKF in combination with a higher (native) resolution EnVar in a dual-resolution EnKF-EnVar hybrid setup (Schwartz et al., 2015; Pan et al., 2018) to save computational cost of operational implementation.
The eventual goal is to provide hourly ensemble perturbations to the operational RAP hybrid DA system. As a pilot study, an EnKF-En3DVar hybrid DA system at 40-km grid spacing was established based on the RAP system in Pan et al. (2014), and more recently a dual-resolution version that uses ensemble covariance downscaled from the 40-km EnKF DA system has also been tested (Pan et al., 2018). Results show consistent advantages of EnKF and hybrid En3DVar over GSI 3DVar, and that the hybrid DA further improves upon EnKF in some aspects of forecasts. These earlier encouraging results prompted the implementation of hybrid En3DVar DA for the operational RAP, but it borrows ensemble perturbations from the NCEP GFS EnKF system instead of producing them using RAP’s own EnKF (Hu et al., 2017; Wu et al., 2017), partly due to computational constraints. The current operational 13-km RAP hybrid DA uses 25% and 75% static and ensemble background error covariances, respectively (Benjamin et al., 2016), with the GFS EnKF perturbations available only every six hours and with delayed availability. For optimal results, it is desirable to operationally implement RAP’s own EnKF DA cycles, and to couple these with the EnVar hybrid DA system. Our prior and current studies represent some of the efforts toward such a goal.
For the cycled ensemble DA system, the lateral boundary conditions are from three-hourly NCEP GFS forecasts but perturbed using the “randomCV” option of the WRF DA system (Barker et al., 2012) following Torn et al. (2006). The same method was used to generate perturbations for initializing the first EnKF DA cycle. For the GSI-based EnKF system, the observation innovations, y − H(x), are calculated within GSI and then used within the EnKF system. The analysis algorithm is the ensemble square root filter of Whitaker and Hamill (2002). Forty-member ensemble forecasts are run for three hours between the analysis times, while 18-h short-range deterministic forecasts are run from the ensemble mean analyses of each cycle. As in Z13, a modified digital filter launch procedure (Lynch and Huang, 1994) where the land surface fields are subject to the same filtering, is used before launching the forecasts. All members use the same physics options, including Grell-G3 cumulus parameterization, Thompson microphysics, RRTM longwave radiation, Goddard shortwave radiation, the MYJ planetary boundary layer scheme, and the RUC-Smirnova land-surface model. In short, the EnKF DA system used is the same as that described in Z13, and the specific configurations are the same as experiment EnKF_CtrHDL of Z13, except for the addition of satellite radiance data and global positioning system (GPS) PW (Smith et al., 2007).
Following Z13, height-dependent background error covariance localization scales are used, in which the horizontal covariance localization scale increases from 700 km at the surface to 1050 km at the model top. The vertical localization scale ln(pcut) depends on the analysis variable. For relative humidity (RH) and temperature, it is set to a quarter and one half of 1.1 at the surface and model top, respectively. For horizontal wind components, U and V, the vertical correlation length is twice as large as that for RH and temperature. A value of 1.6 is used for surface pressure observations. These settings were found to produce the best results in Z13 based on sensitivity experiments. Here, the vertical localization scale for radiance observations is 1.6. The covariance inflation used is the same as in Z13, containing a fixed and an adaptive part (Anderson, 2009; Whitaker and Hamill, 2010). For the experiments with satellite radiance assimilation, the coefficient of the fixed part, b in Eq. (5) of Z13, is increased from 0.1 (no satellite radiance assimilation) to 0.16 (AMSU-A radiance assimilation) and 0.18 (full radiance assimilation). This is because the ensemble spread tends to be reduced when more satellite observations are assimilated. More details on other settings can be found in Z13.
For radiance assimilation, GSI uses the Community Radiative Transfer Model (CRTM) developed by the Joint Center for Satellite Data Assimilation (Han et al., 2006; Weng, 2007) as its observation operator. CRTM is capable of simulating most geostationary/polar-orbiting satellite observations, covering microwave and infrared frequencies.
-
The radiance data assimilated in this study include AMSU-A, AMSU-B, AIRS, as well as the Microwave Humidity Sounder (MHS) and High-resolution Infrared Radiation Sounder (HIRS/3 and HIRS/4) instruments on polar-obitting satellites. These data were part of the operational RAP and GFS observation data streams. For a regional, frequently updated forecasting system, radiance data preprocessing is an important step towards obtaining a positive impact of radiance assimilation (Lin et al., 2017a). To be consistent with the three-hourly assimilation cycles tested in this paper, the radiance data are reprocessed into three-hourly data batches, including data within the ±1.5-h windows centered on the analysis times. In this study, clear-sky radiance data are assimilated. GSI uses the minimum residual method to detect clouds (Eyre and Menzel, 1989). Additionally, for the quality control step, GSI uses outliers to reject observations along the scan edges. Table 1 gives the GSI default values for each sensor. Note that some satellite sensors have higher spatial resolution, and hence more scan points than others. In general, about 10% of the scan points along both left and right edges were removed.
Sensor name Number of scan points Left bound Right bound AMSU-A 30 4 27 AMSU-B 90 10 81 MHS 90 10 81 HIRS 56 7 50 AIRS 90 10 81 Table 1. Lists of scan points used by GSI for each satellite sensor. Points outside the left and right bounds are removed.
Figure 1 shows the number of assimilated data for each channel used and the corresponding height level of the maximum response function. AMSU-A and MHS are at microwave frequency, while AIRS and HIRS are at infrared frequency. For infrared sensors, the number of assimilated data is smaller near the surface compared to higher-level channels (Figs. 1a, e and f). This is especially true for infrared sensors such as HIRS4 (Figs. 1e and f). AIRS behaves similarly but benefits from more channels available near the surface (Fig. 1b), and the total number of near-surface data is still large (Fig. 1a). In fact, for this study, the total number of assimilated AIRS data is the largest. For microwave sensors, the number of surface channel data is not as many as for upper-level channels, but is still comparable (Figs. 1c and g). AMSU-A mainly provides atmospheric temperature profiling with data from several satellites and provides by far the best coverage. Among all instruments, the amount of assimilated AMSU-A data is the second largest. In all, the satellite radiance data from multiple sensors and satellites provide a reasonably good coverage in the 3D model domain over a period of time.
-
Online bias correction methods can be classified into variational versus non-variational approaches. Variational bias correction estimates the correction coefficients simultaneously with state estimation within a variational framework. Specifically, in a 3DVar framework, the cost function, following Dee and Uppala (2009), is given by
where
is the correction to the departure of simulated radiance
${H_2}({{x}})$ from radiance observations${{{y}}_2}$ . Here, the modified 3DVar cost function$J$ includes the estimation of bias correction coefficients in coefficient vector$\,{\beta }$ , which are analyzed together with state vector x. We write the observation term in two parts: one for radiance data${{{y}}_2}$ , and one for all other observations${{{y}}_1}$ , which are associated with the corresponding error covariances R1 and R2, and observation operators H1 and H2.${{{x}}_{\rm{b}}}$ is the background state vector, and x is the state vector to be analyzed.${{\beta }_{\rm{b}}}$ is the prior estimate of bias correction coefficient.${{{B}}_x}$ and${{{B}}_\beta }$ are the error covariance matrices for the background state and bias parameters, respectively. The radiance data are corrected using the bias correction vector${{b}}$ given by Eq. (2), which depends on predictors in P. Those predictors are usually air-mass dependent (Flobert et al., 1991). Subscript i represents the index of predictors, Np is the number of predictors. pi represents individual column vector i of matrix P. In the version of GSI used in this study, there are five predictors for bias correction, including the global offset, zenith angle predictor, cloud liquid water (CLW) predictor, temperature lapse rate (TLR) predictor, and the square of the TLR predictor. In the version used, the scan-angle component is still updated outside GSI (Derber and Wu, 1998) using the moving average of the weighted differences between the quality controlled radiance observations and the first guess. More recent versions of GSI have incorporated this term into the variational estimation (Zhu et al., 2014).Equations (3) and (4) give an example of how the predictor terms CLW and TLR are calculated within GSI for the Metop AMSU-A sensor:
Here,
${\rm{B}}{{\rm{T}}_{{\rm{sim}}}}$ and${\rm{B}}{{\rm{T}}_{{\rm{obs}}}}$ are the simulated and observed brightness temperature, respectively. The numbers indicate channels 1 and 2; d1 = 0.754 and d2 = −2.265 are weight coefficients. Based on Eqs. (3) and (4), CLW is not channel-dependent, while TLR is. In Eq. (4),$\tau $ is the transmittance and T is the model background temperature. Index k represents the vertical level, while N is the number of vertical levels. T0 is surface temperature. Figure 2 shows an example of TLR calculated using the US standard atmosphere profile for the Metop AMSU-A sensor. The changes of transmittance between levels$\Delta \tau = \tau (k) - $ $ \tau (k + 1)$ are very similar to vertical weighting functions (not shown).$\Delta \tau $ is always positive (Fig. 2a). For TLR, channels 1−6 are positive. Among these, channels 4 and 5 have clearly larger TLR values than other channels. As the temperature gradient$\Delta T = T(k - 1) - T(k + 1)$ turns to be negative at some high levels, the TLR of channels 8−11 become negative. Channel 15 is an exception since it also peaks near the surface.Figure 2. (a) Vertical profiles of transmittance changes for Metop AMSU-A, (b) vertical profile of temperature, and (c) bias corrected term of the temperature lapse rate. All figures are plotted using US standard atmosphere within the CRTM. Only channels that are finally assimilated are given.
The other non-variational bias-correction approach uses the EnKF analysis equation to update the bias-correction coefficients or predictors, and it is performed within each cycle before state estimation (Miyoshi et al., 2010). Setting the derivative of cost function J in Eq. (1) with respect to
$\,{\beta }$ to zero yields the following optimal update equation for the increment of bias correction coefficient vector,$\delta {\beta }$ :where only satellite data in vector
${{y}_2}$ are used. The above formulation can be shown to be mathematically equivalent to the 3DVar solution for a linear and Gaussian system, which also has the same form as the EnKF analysis update equation (with the inclusion of online bias correction) (Kalnay, 2002). Here, Bβ is assumed to be diagonal and the diagonal elements are assumed to be 0.1, which is also the GSI default value. Within the EnKF system, we apply bias correction to observation priors using prior estimated bias correction coefficients. Equation (5) is used to update$\,{\beta }$ after state estimation. After$\delta {\beta }$ is obtained, the bias correction coefficients$\,{\beta }$ are updated according toWe note here that the coefficient estimation in Eq. (5) appears independent of the state vector estimation, while the estimation of the coefficients using the cost function in Eq. (1) involves state vector x. This is because the EnKF algorithm is sequential; the correction coefficients are updated using the updated state vector. In comparison, in Eq. (1), state vector x and the coefficient vector are updated simutaneously. However, because
${{B}_x}$ and${{B}_\beta }$ are uncorrelated, the state and coefficient estimations, even in the variational framework, are not strongly linked. Still, the algorithm differences can lead to differences in results in practice, and it is a goal of this paper to compare their relative performance.In our EnKF experiments, we use either the variational procedure that minimizes the cost function Eq. (1) or the EnKF procedure that updates the correction coefficients using Eq. (5). In Table 1, we label the former procedure “variational” and the latter procedure “EnKF ”. When the variational procedure is used to estimate bias correction coefficients for EnKF experiments (for comparison purposes), GSI 3DVar is run with full observations using the EnKF ensemble mean forecast as its background in each analysis cycle; in other words, the GSI 3DVar cycles piggyback on the EnKF cycles; the GSI 3DVar solver is therefore used purely for the purpose of variationally estimating the bias correction parameters, while its state estimation is discarded. For GSI 3DVar, two outer loops and a maximum of 50 inner iterations are imposed.
As pointed out earlier, in theory, the variational and EnKF parameter updates should be mathematically equivalent, for linear systems with Gaussian errors at least. The key difference lies in that, within a variational framework, the parameter and state estimations are performed simultaneously, while within EnKF they are updated sequentially. In addition, the variational minimization involves linearization of the observation operators.
2.1. Brief review of the GSI-based EnKF system
2.2. Assimilation data
2.3. Introduction to two online bias correction schemes
-
To compare with the results without satellite radiance assimilation, we continue to use the same test period as Z13, starting at 0000 UTC 8 May 2010 and ending at 2100 UTC 16 May 2010, over a total of nine days. This period features a variety of active convective weather, including the 10 May Oklahoma tornado outbreak featuring strong mesoscale forcing, a mesoscale convective system (MCS) on the night of the 11th, a cold-season type Front Range upslope low-pressure precipitation event, and some southeast propagating MCSs across Texas. The DA experiments are run in continuous three-hourly updated cycles. The assimilation domain has a horizontal grid spacing of 40 km, with 207 × 207 × 50 grid points, and covers the whole of North America (Fig. 3). The model top is at 10 hPa. As described in Z13, the use of a relatively short nine-day period was mainly due to computational resource constraints; however, this length appears long enough for the EnKF state estimation and bias correction to stabilize [as discussed in Pan et al. (2014)].
Table 2 lists all experiments presented in this paper. We first investigate initial bias correction coefficient issues and test the two bias correction methods in the EnKF system described in section 2c. AMSU-A radiance data are assimilated since they are available from different satellites and the number of data available within the domain is much more stable than other types of radiance data. Experiments EnKF_AMSUA_VarB0, EnKF_AMSUA_VarB, and EnKF_AMSUA_VarB2 use the variational bias correction with GSI 3DVar. For the variational bias correction,
${y}$ also contains conventional observations during the update of bias correction coefficients$\,{\beta }$ . In experiment EnKF_AMSUA_VarB_SP, we use only radiance observation to update$\,{\beta }$ . Experiments EnKF_AMSUA_EnB use EnKF parameter estimation. Here, the number “0” in the name indicates the use of zero$\,{{\beta }_{\rm{b}}}$ in the first analysis cycle; otherwise, the recycled initial estimate of$\,{{\beta }_{\rm{b}}}$ is used. In the recycled case, the entire nine days of DA cycles are repeated, starting from the$\,{{\beta }_{\rm{b}}}$ estimated at the end of the previous pass of the nine-day cycled DA. The numbers “1” and “2” represent the recycled initial estimate of$\,{{\beta }_{\rm{b}}}$ from the zero and first run, respectively. For simplicity, the number 1 is ignored. This recycled procedure is an attempt to improve the initial guess of$\,{{\beta }_{\rm{b}}}$ when the total length of DA cycles is relatively short. Experiments with “0” are considered as spin-up runs for the correction parameter estimation and are included here mainly to see if a better initial estimate of$\,{{\beta }_{\rm{b}}}$ , which is typically available in long-running DA cycles, can improve the results, given the relatively short nine-day DA period we are using. Note that in continuously cycled operational DA systems, this procedure is not necessary, and is used here to help alleviate the potential problem of the short DA period.Experiment Satellite radiance Bias correction method EnKF_AMSUA_VarB0 AMSU-A Variational procedure, ${{\beta }_{\rm{b}}}$= 0 EnKF_AMSUA_VarB AMSU-A Variational procedure, using recycled ${{\beta }_{\rm{b}}}$ from EnKF_AMSUA_VarB0 EnKF_AMSUA_VarB2 AMSU-A Variational procedure, using recycled ${{\beta }_{\rm{b}}}$ from EnKF_AMSUA_VarB EnKF_AMSUA_VarB_SP AMSU-A Similar to EnKF_AMSUA_VarB, but ${{\beta }_{\rm{b}}}$ is obtained using only radiance observations EnKF_AMSUA_EnB AMSU-A Similar to EnKF_AMSUA_VarB, but using EnKF procedure EnKF_NoSat − − EnKF_ALL AMSU, AIRS,
MHS, HIRS 3/4Similar to EnKF_AMSUA_VarB, but using full radiance Table 2. List of data assimilation experiments. Here,
${{\beta }_{\rm{b}}}$ is the first guess of the bias correction coefficients vector used in the first analysis cycle. “VarB” in the experiment names stands for variational bias correction, while “EnB” stands for EnKF-based bias correction.We then test the impact of the full set of radiance data in the EnKF system. EnKF_ALL assimilates all of the AMSU-A, AMSU-B, AIRS, MHS, HIRS3 and HIRS4 data, and is compared with EnKF_NoSat without satellite data (but with PW data included). The nine-day DA cycles are run twice, with the second time starting from the recycled estimate of
$\,{{\beta }_{\rm{b}}}$ from the first pass, as in the “1” experiments. To see if the DA results produced on the 40-km grid can improve forecasts at the 13-km RAP operational resolution, 13-km forecasts from interpolated analyses of the above two experiments are run twice a day starting at 0000 and 1200 UTC, and the precipitation forecasts are evaluated. -
The 40-km forecasts are verified against sounding and surface observations. As in Z13, root-mean-square error (RMSE) and relative percentage improvement (RPI) are used as the verification metrics. The RPI is defined as
where the subscripts “sat” and “nosat” represent the experiments with and without satellite radiance assimilation. Positive RPI values mean reduced error or improved analyses/forecasts due to radiance assimilation.
Following Pan et al (2014), the statistical significance of RMSE differences and RPIs is determined by using bootstrap resampling (Candille et al., 2007; Buehner and Mahidjiba, 2010; Schwartz and Liu, 2014). The same method with 3000 randomly selected samples are used. For those samples, their mean and a two-tailed 90% confidence interval from 5% to 95% is calculated. For RPIs, if the bounds of a 90% confidence interval between the forecast pair are all higher than zero, it means that the improvement from the satellite assimilation experiment over the corresponding no satellite experiment is statistically significant at the 90% confidence level. In contrast, if the bounds of the 90% confidence level includes zero, it indicates a statistically insignificant difference (Xue et al., 2013; Pan et al., 2014).
For the 13-km forecasts, we focus on precipitation verification using the NCEP stage IV precipitation data as observation. Two verification scores—the Gilbert skill score (GSS) (Gandin and Murphy, 1992) (also known as the equitable threat score) and frequency bias (FBIAS)—are calculated over the contiguous US.