Online Nonlinear Bias Correction in Ensemble Kalman Filter to Assimilate GOES‐R All‐Sky Radiances for the Analysis and Prediction of Rapidly Developing Supercells

The present study introduces the online non‐linear bias correction for the assimilation of all‐sky GOES‐16 Advanced Baseline Imager (ABI) channel 9 (6.9 μm) radiances in a rapidly cycled EnKF for convective scale data assimilation (DA). This study is the first to explore the use of the radar reflectivity as the anchoring observation for ABI all sky radiance assimilation. The online and offline nonlinear bias correction methods are compared and evaluated for a case of rapidly developing supercells over Oklahoma and Texas. The analysis and background of the online bias correction perform better than the offline approach during the suppression of spurious clouds and the establishment of non‐precipitating and precipitating regions when the supercell storms are observed to develop. The online approach not only improves the analysis and background over the radar anchored region but also the unanchored non‐precipitating regions compared to the offline approach. Both quantitative and subjective verification of the deterministic forecasts showed consistent superior performance from the online bias correction over the offline approach. Diagnostics reveal that the online bias correction retains useful information in the innovation, which in turn improves subsequent analysis, background and background ensemble spread for both the thermodynamic and dynamic fields. The effect is accumulated during the DA cycling that is responsible for the superior analysis and forecast of the supercells.

The combined assimilation of radar and geostationary infrared observations were first performed by several OSSE studies (Cintineo et al., 2016;Jones et al., 2013Jones et al., , 2014. Their results confirmed that the radar reflectivity observations combined with the water vapor or cloud sensitive BT observation results in a better analysis and forecast for humidity, cloud and precipitation hydrometeors than assimilating radar reflectivity alone. Most of the above OSSEs focus on extra-tropical cyclone systems. Most recently, several studies were conducted to assimilate real ABI radiance with radar observations for rapidly developing thunderstorm and supercell events (e.g., Johnson et al., 2021;Zhang et al., 2019Zhang et al., , 2018 at convection permitting resolution. These studies found that for convective scale data assimilation (DA) systems, the frequent assimilation of all-sky BT along with radar observations results in better representation of storms in the analysis and improved forecast lead time than radar only experiments.
Satellite radiance (Radiance and BT are interchangeably used in the paper) observations may contain systematic biases. Past studies suggested identification and removal of such biases during DA can be critical or beneficial for the analysis and subsequent forecasts (Aravéquia et al., 2011;Dee & Uppala, 2009;Kazumori, 2014;Zhu et al., 2019). Initial efforts were steered toward identification and removal of this bias only in clear-sky regions (Kelly & Flobert, 1988). The air-mass and scan angle dependent nature of the clear-sky biases motivated studies to account for the spatial variation in the bias as a function of various bias predictors (e.g., Eyre, 1992). With the emergence of cloudy radiance observation assimilation, in addition to examining the use of various predictors, studies have also begun to explore the use of various linear or nonlinear bias correction functions for bias correction (Otkin et al., 2018;Otkin & Pothast, 2019).
Satellite radiance bias correction can also be broadly classified into offline and online methods. The general idea behind an offline bias correction method is to fit the bias correction function with the accumulated innovation (observation minus background) statistics (e.g., Harris & Kelly, 2001). This approach has been shown to be effective in improving the data quality and reducing the air-mass dependent bias. In the presence of model bias, however, an offline approach may have difficulty distinguishing between the model bias and observation bias (Auligne et al., 2007;Eyre, 2016). An online radiance bias correction method was proposed and implemented first in variational DA (Dee, 2004;Derber & Wu, 1998;Zhu et al., 2014) and then in ensemble Kalman filter (EnKF; Fertig et al., 2009;Miyoshi et al., 2010). In such an online bias correction approach, the bias correction coefficients and the model states are updated simultaneously during the DA cycling. The simultaneous assimilation of bias-free anchoring observations and bias-bearing satellite radiances in the online approach allows the bias of the radiance observation to be separately estimated from the bias of the model (Auligne et al., 2007;Eyre, 2016). In addition, compared to the offline approach, the online approach has been proven to be more adaptive to the changes in the model and observation characteristics (Auligne et al., 2007).
Although early studies have suggested the importance of correcting the bias of the radiance during DA, limited efforts were made to explore the radiance bias correction for convective scale all-sky infrared radiance DA. Most studies assimilating real all-sky infrared radiance for convective scales either perform no bias correction (e.g., Zhang et al., 2018Zhang et al., , 2019 or use a simple offline approach where domain-wide averaged innovation was used to approximate the bias (e.g., Johnson et al., 2021). Otkin and Pothast (2019) implemented an offline approach where a nonlinear function was used to fit innovations. This study showed the positive impact of nonlinear bias correction on upper-level cloud systems associated with widespread precipitation.
Although early studies using a simple model (e.g., Eyre, 2016) or in a clear air radiance assimilation context (Auligne et al., 2007) both portrayed the superior behaviors of online bias correction approach relative to an offline approach in the presence of anchoring observations, there has not yet been an attempt to use online bias correction methods for all-sky radiance assimilation in convective scale DA. In the present study, we implement online nonlinear bias correction for all-sky ABI radiance DA for a case of rapidly developing supercells. Such an online bias correction approach is implemented with radar reflectivity as the anchoring observations and with a rapidly updated EnKF. Radar reflectivity is explored as the anchoring observation given its matching high temporal resolution as the ABI BTs. In addition, different from early online bias correction studies, the nonlinear polynomial form of bias correction function was adopted (Otkin et al., 2018). We term this approach as "online nonlinear bias correction." As a first study, we assimilate the water vapor and cloud sensitive ABI channel 9 10.1029/2021MS002711 3 of 25 radiances sensitive primarily to clouds and water vapor in the middle and upper troposphere. The online nonlinear bias correction is then comprehensively compared with the corresponding offline approach. Furthermore, the optimal choice of predictor for nonlinear bias correction is examined for the rapidly developing supercell case. The proposed online nonlinear bias correction method for the all-sky radiance DA for convective scale prediction has not been published in earlier studies. To demonstrate this new approach in depth, we first use a single high impact weather event and its associated detailed diagnostics to facilitate the understanding of and to reveal the impact of the approach. Such a single case study approach has been used in early studies when a new approach is introduced (e.g., Honda et al., 2018;Minamide & Zhang, 2018;Wu et al., 2019;Zhang et al., 2018). The limit of a single case study and the need to perform experiments with more cases for future studies are discussed in Section 5.
The paper is organized as follows: Section 2 discusses the theory behind offline and online bias correction methods used in this study followed by a brief description of the rapidly developing supercell event and experiment design in Section 3. Section 4 discusses the effect of using different bias predictor variables on the convective scale DA system. In addition, the performance of offline and online bias correction methods for the analysis and forecast of the rapidly developing supercell are evaluated and diagnosed. Conclusions and scope for future work are presented in Section 5.

Offline and Online Bias Correction Methods
In this study, a nonlinear online bias correction method is proposed and examined for GOES-R ABI all-sky radiance DA. In this section, we first review a theory derived from a simplified framework to understand online and offline bias corrections (Eyre, 2016). We then discuss practical considerations examining the theory in a realistic convective scale DA system. Finally, we describe the nonlinear form of offline and online bias correction methods adopted for all-sky radiances in this study.

Theoretical Differences Between Offline and Online Bias Correction Methods and Their Practical Application
The bias correction methods in this study aim to correct the observation biases rather than the biases that exist in the model. Distinguishing the observation and model biases requires unbiased reference observations to be assimilated along-side the bias-bearing observation. These unbiased observations help differentiate the model and observation component of the biases thereby preventing the drifting of the analysis toward the model climatology. These unbiased reference observations are termed as "anchor observations" (Cucurull et al., 2014). The simplified framework presented in Eyre (2016) establishes a theory for understanding the differences between offline and online bias correction methods. The framework considers two observations "y 1 " and "y 2 " which observe the model state "x". "y 1 " is an anchor observation and "y 2 " is the bias corrected observation. Anchor observations are expected to be unbiased and therefore can help distinguish the model and observation component of the biases. The analysis at any DA cycle is written as a weighted sum of model background and the observations. One major assumption regarding the bias in the simplified framework is that they do not vary with time. The difference of the analysis and background biases between the offline and online approach derived by Eyre (2016) are summarized in Table 1. It can be seen that, in the absence of any anchoring observation (when w 1 = 0 in Table 1), both offline and online bias correction behaves the same way where the bias in background and analysis is equal to model bias. For both the offline and online bias correction approaches, the background and analysis biases are reduced when anchoring observations are assimilated. For a given weight to the anchoring and bias corrected observation, the online bias correction method results in less biases in the analysis and background compared to the offline method.
This simplified framework reveals the need of the co-existence of an unbiased anchor observation along with the bias corrected observation and the potential advantage of an online bias correction approach relative to an offline approach. In the past studies, radiosonde observations and GPS Radio Occultation observations have been used as anchor observations (Cucurull et al., 2014). The utility of such observations as anchoring observations in convective scale DA with frequent cycling (e.g., 10 min in this study), however, is limited, given their relatively low time resolution. Although, the radar reflectivity observations can have biases as a result of antenna side lobes, atmospheric absorption and attenuation (Vivekanandan et al., 2003), past studies suggest the quality 4 of 25 controlled MRMS reflectivity observations can serve as the anchoring observations for the ABI radiance assimilation. First, the US NEXRAD WSR-88D system requirement (Zrnic, 2012) specifies a precision in reflectivity as less than 1 dBz, The NEXRAD observations are then processed through several quality control procedures (Tang et al., 2020) to form 3D MRMS reflectivity mosaic, which are assimilated in this study. During the quality control, errors are further reduced. Past studies have used these MRMS reflectivity observation (Biswas & Chandrasekar, 2018;Gourley et al., 2003;Ice et al., 2005) or precipitation (Carr et al., 2015;Upadhyaya et al., 2020) as ground truth in order to quantify biases in the satellite reflectivity or precipitation observations. We therefore explore the use of radar reflectivity as anchor observations for the assimilation of GOES-R ABI radiances. There are a few aspects that need to be considered. The simplified framework assumes that both the anchor and bias corrected observations are directly observing the same model state. In the realistic convective scale DA system, radar reflectivity and ABI radiances are remote sensing observations that indirectly measure atmospheric or model states. Different observation operators are used to relate what radar reflectivity and ABI "see" with the model states. Additionally they are sensitive to different variables and contribute to the analysis increment in different regions of the cloud. However, radar reflectivity and ABI cloud sensitive radiances are physically connected. For example, the increments created by radar reflectivity on precipitating hydrometeors and updrafts at low-and mid-levels are physically linked to the formation of high-level clouds which are seen by ABI radiances. Therefore, we expect radar reflectivity observations can still serve as an anchor over cloudy regions where there are positive reflectivity. The time scale for the formation of a rapidly evolving cumulus cloud system is approximately 10 min or less (Bluestein, 1992). With a 10 min DA cycling period used in this study, the lag between radar reflectivity and radiance observations is included and the two observations are dynamically linked over the precipitating regions. We also expect zero reflectivity can serve as an anchor over spurious precipitating regions because the increments created by radar reflectivity on spurious precipitation also affects the suppression of high-level spurious clouds which are seen by ABI radiances. Radar reflectivity observations cannot act as anchor observation in the following regions. In the regions where both the model and the observation are clear air, the humidity increments created by ABI radiances are not constrained by the radar reflectivity observations as radar reflectivity is sensitive only to precipitating hydrometeors. In the regions where both the model and observations have non-precipitating clouds, radar reflectivity cannot serve as an anchor because the innovation for radar reflectivity observation is zero.

Non-linear Offline Bias Correction
The non-linear bias correction method introduced in Otkin et al. (2018) is an offline method wherein the statistics of observation departure from the background (dY) is used to estimate the coefficients of a polynomial. In this study, the coefficients of the polynomial are calculated using the first-guess radiance departure values collected from the current DA cycle. Following Otkin et al. (2018), the observation departure from the model first-guess is expressed as a cubic polynomial function of a predictor as, Offline bias correction Online bias correction Statistics of firstguess departures (O-B) are used to estimate the bias.
Bias is estimated at each analysis step along with the estimated state The bias correction to "y 2 " is given by, 2 = 2 − The bias correction to "y 2 " is given by, 2 = 2 − Assuming the anchor observations are unbiased that is, b 1 = 0, then the ratio background/analysis to model bias is given by (1− 2 )+ 1 = (1− 1 ) (1− 2 )+ 1 Note. The bias in model, background and analysis are represented as "b m ", "b b " and "b a " respectively. Similarly, "b 1 " and "b 2 " refers to the bias in anchor and bias corrected observation. "w b ", "w 1 " and "w 2 " refers to the weights given to the background, anchor observation and bias corrected observation respectively. The relaxation rate of the model is represented by " ′ ." Table 1 Summary of Major Differences Between Offline and Online Bias Correction Based on Eyre (2016) In the above equation, i denotes each observation. "m" represents the total number of observations, "dY" represents the observation departure vector, "z" is the bias predictor, "c" is a constant which is set to the average value of the bias predictor and "b" represents the coefficients of the polynomial. The above equation can be represented in the matrix form as, = Here, "A" represents the matrix containing "z i -c" terms for each of the observations, "b" represents coefficients of the polynomial and is the vector containing "dy i " elements. The solution to the above system is given by Otkin et al. (2018), Here, "α" is a constant regularization parameter which is set to a small value (10 −9 ).
The above method represents fitting a cubic polynomial function of bias predictors to the model first-guess departure. For all the experiments performed in the paper, the order of polynomial is set to three following Otkin and Pothast (2019).

Non-Linear Online Bias Correction
In an online bias correction approach, the estimation of bias correction coefficients is performed along with the state update step in the DA system. This was achieved in the variational DA system by augmenting the control vector with the bias coefficients (Derber & Wu, 1998). The EnKF version of the online bias correction was implemented by Miyoshi et al. (2010) in Japan Meteorological Agency operational global analysis and prediction system. An equivalent online bias correction approach for clear-sky radiances was previously implemented in Gridpoint Statistical Interpolation based serial Ensemble Square Root Filter (EnSRF; Hamill et al., 2011) where the bias ( , ) is expressed as a linear function of predictors p and bias correction coefficients β. Different from Hamill et al., 2011, this study assimilates all sky radiances. Therefore the linear function is replaced with the cubic polynomial form similar to the non-linear offline method in Equation 4 The state and bias coefficient updates are coupled through the inner loop iterations. The following steps are followed in the actual implementation: 1. With the latest available estimate of the bias coefficients, the bias is removed from the ensemble mean innovations prior to assimilation (see innovation term in Equation 4). 2. Observation prior and perturbations are updated using the serial EnSRF Equations 4 and 5 (Liu et al., 2018;Whitaker & Hamill, 2002),̄=̄+ where K and ̃ represent the Kalman gain and reduced Kalman gain. H and represent the observation operator and observation respectively. ̄ , ̄ , ′ and ′ represent the ensemble mean analysis, ensemble mean background, analysis perturbation and background perturbation respectively.
3. Once all the observations are assimilated, using all the updated observation prior values, the bias coefficient update is performed using Equation 6: where denotes the covariance of the bias coefficients and it is assumed diagonal, p denotes the matrix of bias predictors at each observation location, and represents the increment to the observation prior, which is obtained from step 2 4. The process from step 1 to 3 is repeated for specified "n" inner loop iterations to obtain the final and therefore b 10.1029/2021MS002711 6 of 25 5. The ensemble mean state (̄ ) and perturbations ( ′ ) are update in the last inner loop iteration using the following serial EnSRF Equations 7 and 8.
It should be noted that in the EnSRF formulation, the bias correction is applied to the ensemble mean innovation as shown in Equations 4 and 8. It is also possible to apply the bias correction approach on a member by member basis using the stochastic EnKF, that is, EnKF with perturbed observations (Houtekamer & Mitchell, 1998).

Overview of Case Study, Model and DA Configuration
In order to test the online non-linear bias correction and compare it with an equivalent offline approach, a rapidly developing supercell event of 18 May 2017 was selected. The case is characterized by a convective dry line along with the high-level trough which lead to the formation of multiple supercells across Texas and Oklahoma. Figure 1 shows the observed channel 9 ABI BTs and composite radar reflectivity observations from 1750 to 1900 UTC. Both ABI and radar reflectivity observations show the initiation and rapid intensification of the storms. The ABI observations sense the formation of towering cumulus clouds well before the strong reflectivity signatures (≧35 dBz, Mecikalski & Bedka, 2006) are seen by the radar observations. The two storms initiated around 1800 UTC (labeled as A and B in Figure 1) developed into long-lived supercells. At around 1900 UTC, elevated convection further north initiated another thunderstorm. This study focuses on the analysis and shortterm prediction of the two long-lived supercells. This event was also used by Johnson et al. (2021) in the context of studying the impact of channel 9 and 10 ABI radiance assimilation with adaptive observation error and additive inflation. Studies have shown (e.g., Burghardt et al., 2014;Cintineo & Stensrud, 2013;Kain et al., 2013;Zhang, Minamide, & Clothiaux, 2016) the challenges for an accurate prediction of the convection initiation and early development of the storm. Such challenges are partially attributed to the limited availability of observations associated with this stage of the storm. The present work therefore focuses on the initiation and rapid intensification stage of the supercell.
The analysis of the rapidly developing supercells is facilitated by a convective scale DA system whose configuration includes 40-member ensemble whose initial and boundary conditions are from the NCEP Global Ensemble Forecast System (GEFS Zhou et al., 2017) with hourly assimilation of conventional observations from 0000 to 1700 UTC. The details of the GSI-EnKF system extended for convective scales can be found in Johnson et al. (2015) and Wang and Wang (2017). As discussed in Section 4.1, the first step is to determine the bias predictor variables which will be used to perform comparison of the online and offline bias correction methods. The conventional DA is followed by 10-min assimilation of radar and ABI from 1710 to 1900 UTC to achieve the first step. Experiments using the CONUS domain and the limited domain covering the supercells are conducted for this step. The former was employed because it includes a variety of weather features and more ABI statistics can be obtained. As discussed in Section 4.1, the results of the optimal predictors remain the same for both domains. Further, the experiments are run up to 1900 UTC in order to ensure that the statistical result stabilizes during the DA cycling period. For the evaluation of non-linear offline and online bias correction methods, assimilation experiments are performed over a limited domain covering Oklahoma from 1750 to 1900 UTC in order to focus specifically on the impact of assimilation on the supercells. A cycling starting time of 1750 is selected because such a time is the first that the radar reflectivity observations are able to anchor the satellite observations (Section 2.1). Advanced Research WRF model version 3.9.1 (Skamarock et al., 2008) is used in this study with physics configuration consisting of Thompson et al. (2008) scheme for microphysics parametrization, Mellor-Yamada-Nakanishi-Niino (Nakanishi & Niino, 2004 for planetary boundary layer, RUC land surface model (Benjamin et al., 2004), and rapid radiative transfer model (RRTMG; Mlawer et al., 1997) for shortwave and long wave radiation parameterizations. In addition to evaluating the analyses of the rapid initiation and strengthening of the storms, the skill of deterministic forecasts initialized from the ensemble mean analyses every 10 min from 1820 UTC to 1900 UTC is evaluated. The length of the DA cycling and forecast used here is consistent with the earlier published supercell DA and prediction studies (e.g., Wang & Wang, 2017Dowell et al., 2011;Putnam et al., 2019). The horizontal covariance localization adopts the Gaspari and Cohn (1999) function with a cut-off radius of 300 km for conventional observations. The ABI Community Radiative Transfer Model (CRTM; Han, 2006) is employed as a forward operator to create model first-guess radiances. In order to simulate the cloudy radiances, apart from temperature and humidity, vertical profiles of snow (Qs), rain (Qr), ice (Qi), hail (Qh) and graupel (Qg) mixing ratios are given as input to CRTM. The scattering properties which are functions of the hydrometeor species, wavelength and effective radii are read from the CRTM look-up tables (LUTs). The LUTs for thermal infrared wavelengths are constructed based on Mie scattering assuming spherical shape for liquid cloud hydrometeors and a modified gamma size distribution (Chen et al., 2008;Simmer, 1994). The ice phase hydrometeors are assumed to be non-spherical hexagonal column with a gamma size distribution (Yang et al., 2005). The primary source of bias in the simulation of cloudy radiance is the assumption of Mie scattering. Apart from this, more bias is introduced in forward operator with the assumed shape and size distribution of ice phase hydrometeors (Yi et al., 2016). In the present study, these biases in the forward operator along with the instrument bias are together considered as observation bias.
Simulated radar reflectivity in WRF are calculated using the Rayleigh scattering approximation (Han et al., 2013;Smith, 1984). The particle shape and size distribution used in the calculation of simulated reflectivity are selected based on the underlying microphysics scheme. For Thompson microphysics scheme used in this study, the graupel and rain hydrometeors are assumed to be spherical and snow is assumed to be fractal-like aggregates. The particle size distribution for all the species is assumed to be exponential. Several studies like Lang et al. (2011), Tao et al. (2016 and Stanford et al. (2017) have highlighted the contribution from approximations used in the microphysics scheme to the bias in simulated reflectivity. The bias arising out of microphysics scheme are the projection of model bias on to the reflectivity space. The main assumption used in the forward operator calculation is the Rayleigh scattering approximation, which is considered a good assumption for S band radars (Doviak, 2006) and therefore the radar forward operator is expected to introduce negligible bias to the simulated reflectivity.

Bias Predictor Selection
Past studies have found that the effectiveness of bias correction for all-sky radiance assimilation can be sensitive to the choice of bias predictors (e.g., Otkin & Pothast, 2019;Otkin et al., 2018). Therefore, a set of experiments are conducted to first determine the best predictors for the nonlinear bias correction. Given the offline approach serves as a baseline, three predictors, namely Observed BT, Simulated BT and Symmetric BT are examined in the offline bias correction mode. The predictors are then used in the online approach. Symmetric BT is defined as an average of observed and simulated BT (Geer & Bauer, 2011;Otkin & Pothast, 2019). Details of the experiments can be found in Table 2.
The non-linear offline bias correction is applied over observed and spuriously simulated cloud region. These regions were identified using the ABI channel 9 BT threshold of 239 K, which was applied to both observed and simulated BT. This threshold was determined using the methodology described in Harnisch et al. (2016) and Johnson et al. (2021)  regions were bias corrected with a constant offset value, which was calculated using the diagnostics from a no bias correction DA cycled experiment, while other regions use the non-linear bias correction to account for the non-linear impacts of clouds.
The large negative bias in ABI observation at 1710 UTC ( Figure 2) is due to the presence of a large spurious cloud region over Oklahoma. This spurious cloud which appears as cold BT region (as shown in Figures 3e, 3i and 3m) is characterized by deep cumulus surrounded by anvil like cloud structure. When observed BT is used as predictor, analysis fits the observation less compared to the background (Figure 2 solid blue curves). Diagnostics show that over parts of spurious cloud region bias is overestimated (not shown). Therefore the sign of innovation changes from positive to negative, leading to spurious humidity and hydrometeor increments. As a result, the domain average RMSI (Root Mean Square Innovation, which is calculated as where N refers to the number of samples, BT obs and BT sim refers to the observed and simulated BTs respectively) of analysis remains higher than the first guess for most of the early DA cycles. Whereas by using simulated or symmetric BT predictor, such spurious increments do not occur. From the seventh DA cycle (1810 UTC) onwards, both bias and RMSI stabilizes in all the three bias predictor experiments. From the stabilized cycles, it can be seen that using simulated BT as a predictor for the non-linear offline bias correction results in the least bias and RMSI compared to other predictors (red line in Figure 2). Before the storm signature appears in the observations, the simulated BT predictor experiment reduces the spurious clouds the most effectively (1730 UTC and 1800 UTC panels in Figure 3). Once the storm signature appears in the observations, the simulated BT predictor experiment better captures the cold BT region seen in reality and therefore better analyzes the storm compared to other predictors (e.g., 1830 UTC and 1900 UTC panels in Figure 3). Additionally, at or near the mature stage of the storm, the analysis from simulated BT predictor experiment is able to create more high level cloud over the anvil region compared to other predictors (Figures 3h, 3l and 3p). Based on these experiments, the simulated BT predictor is chosen for further study. It is noted that the superior performance of the simulated BT predictor is not dependent on the size of the domain as discussed in Section 3. As discussed earlier, the offline approach serves as a baseline to be compared with the online approach. Therefore, the optimal bias predictor from the offline approach is also used for online bias correction experiments. The range of predictors experimented is based on early studies of Otkin et al. (2018) and Otkin and Pothast (2019). Attempts to use other variables as bias predictor for all-sky radiance bias correction should be planned in a future study.

Online and Offline Bias Correction Experiment Results
Experiments were conducted over a limited domain covering Oklahoma in order to evaluate the difference between the online and offline methods of bias correction specifically on the case of rapidly developing supercells. The DA cycling starts at 1750 UTC and ends at 1900 UTC. 1750 UTC is selected as the start of the DA cycling for the following reasons. Before 1750 UTC, over the DA domain, the model contains non-precipitating spurious clouds and the radar reflectivity observations show no precipitation. As discussed in Section 2.1, radar reflectivity observations are unable to anchor the ABI radiance. From 1750 UTC onwards, the model began to show a narrow spurious precipitation region in region of no observed precipitation. Then non-zero positive radar reflectivity observations start to appear as the storm starts to initiate. As discussed in Section 2.1, radar reflectivity is expected to serve as the anchor beginning at 1750 UTC. Moreover, in order for radar reflectivity to provide constraint on calculation of bias correction coefficients, it is necessary to consider grid points where both radar reflectivity and ABI BT are contributing to the analysis increments. This is conceptually similar to the previous work on bias correction by Auligne et al. (2007) which uses radiance observations in the vicinity of radiosonde observations to calculate the bias correction coefficients. Therefore, the bias correction coefficients for both non-linear offline and online bias correction are calculated and applied over ABI radiances in the vicinity of radar anchored regions.
In addition to the experiments with bias correction, a cycled DA experiment was performed from 1750 to 1900 UTC assimilating radar reflectivity and ABI ch9 radiances without any bias correction. Comparison of the analysis BT from the no bias correction experiment with the online and offline bias correction experiments shows that bias correction more effectively suppresses spurious clouds and spinups the real storm ( Figure 6). This result is consistent with the early studies such as Otkin and Pothast (2019). Forecasts corresponding to the experiments with bias correction are also more skillful than those without bias correction (Figures 7 and 8). Given the objective of the study is the comparison of two bias correction approaches rather than the comparison with and without the bias correction, detailed observation and physical space diagnostics as shown in Section 4 are performed only for the offline and online bias correction experiments.

Impact on the Analysis and First Guess
Sawtooth plots of domain averaged innovation bias and RMSI can quantitatively compare the online and offline bias correction approaches during DA cycling. Figure 4 demonstrates such results from 1750 UTC to 1900 UTC. Since radar and ABI observes different regions of cloud, bias and RMSI with respect to both radar reflectivity and ABI radiances are presented. To separately diagnose the radar reflectivity anchoring effects on different regions of the clouds, the innovation bias and RMSI over radar anchored region (observed and spurious precipitation regions), non-precipitating cloud region, observed cloud region and spurious cloud region are presented.
The 8 cycles shown in Figure 4 can be broadly divided into the suppression of spurious cloud stage (cycles 1-4, 1750-1820 UTC) and establishment of the storms (cycle 3-8, 1810-1900 UTC). At the start of DA cycling, due to the presence of spurious cloud cover, the difference between the observations and the model shows a positive bias in radar reflectivity observation space and a negative bias in ABI cloudy radiance observation space (Figures 4a and 4c). As the spurious clouds are suppressed in the follow on DA cycles (cycles 2-4), the biases are reduced. In the subsequent DA cycling, when the storm is developing rapidly, the model spins up the storm slower compared to observations, as a result the innovation bias in ABI radiance becomes positive and similarly the innovation bias in radar reflectivity observation becomes negative due to less cloud cover (Figures 4a and 4c). Starting from cycle 4, the DA tries to add the observed storms in the analysis but due to the rapid intensification, there is a lag between the observed and analyzed storms. This lag between the analysis and observation is evident in the reflectivity sawtooth plot as a negative bias. Further, the apparent bias of reflectivity even at the end of the DA cycling period is contributed to by the continued initiation of new storms that the DA then needs to add, despite the two primary storms of interest being well analyzed by this time. The initial absence of these storms in the analysis as they first appear in observations contributes to a continued negative "bias" in the DA diagnostics.
Starting at 1810 UTC, the experiment using the online bias correction approach results in smaller bias and RMSI in the first guess verified against the anchored radar observations than that using the offline approach (Figures 4a  and 4b). It should be noted that the bias plotted in Figure 4a is the bias of first guess minus the observed reflectivity. This result suggests that the online approach directly improves the radar anchored regions. Similar improvements from online bias correction against offline method was observed when the relative bias in reflectivity was used (not shown). Diagnostics in ABI identified spurious cloud regions are presented in Figures 4c and 4d. Major improvement from online approach in the suppression of spurious clouds occur from the second to the fourth DA cycles (from 1800 UTC to 1820 UTC). In the first few cycles, more high level spurious solid hydrometeors along with a narrow spurious precipitating region are present in the background. The no-precipitation radar reflectivity observations over the narrow spurious precipitation region can provide an anchoring. As a result, the online bias correction suppresses the spurious clouds more effectively than the offline approach.
The bias and RMSI in the observed cloudy region is separately evaluated in Figures 4e and 4f. This result primarily evaluates the region associated with the observed, rapidly developing supercells. As shown later (Figure 6), the model delays capturing the supercells in general. As more positive radar observations become available to serve as an anchor, the innovation bias starts to stabilize from 1820 UTC (cycle 4). In contrast, the bias in offline approach stabilizes only from 1840 UTC. Moreover, the stabilized bias and RMSI in the online bias correction is smaller compared to the offline approach. This result suggests that as the storm intensifies the observed cloud region is better analyzed using the online approach than the offline approach.
To understand the effect of the online and offline bias correction over the un-anchored region, the statistics over the non-precipitating cloud region is evaluated in Figures 4g and 4h. Although radar only anchors the precipitation region, the first guess of the online bias correction fits the non-precipitating cloud observations more closely than the offline approach. For example, from the second DA cycle (1800 UTC) up to the fourth or fifth cycles (1820-1830 UTC), the spurious non-precipitating clouds are suppressed more effectively in the online bias correction, consistent with a smaller RMSI and a slightly reduced bias compared to offline approach. From 1840 UTC (cycle 6), although the offline approach fails to spin up the observed non-precipitating cloud seen as a positive innovation bias and a large RMSI, the online bias correction is able to capture the observed non-precipitating clouds faster seen as a near-zero bias and reduced RMSI in Figures 4g and 4h. Probabilistic distribution of first guess deviation from the ABI radiance observations as a function of observed channel 9 BT is analyzed to further understand the effect of online and offline bias correction approaches following Otkin et al. (2018). Data accumulated from the 1750 to 1900 UTC DA cycles during the rapid development of the supercells were used to calculate the probability distribution shown in Figure 5. As there are more clear  sky than cloudy pixels, there is high percentage of probability of higher BTs. Furthermore, due to the presence of spurious cloud pixels, there is a spike in the probability at clear-sky BTs (Figures 5a and 5b). As expected, for the offline experiments, the systematic deviation in cloudy radiance are larger compared to that of clear-sky radiances (Otkin & Pothast, 2019). The online approach exhibits less systematic deviation over the cloudy region (observed BT≤239 K) than the offline approach (Figures 5a and 5b). This result suggests the effectiveness of using the radar reflectivity anchoring observations to analyze the observed cloudy regions. As discussed in Figure 4, although radar only anchors the precipitation regions, the non-precipitation cloud regions can be indirectly affected. Probabilistic distribution of the first guess deviation for the non-precipitation ABI radiance observations is shown in Figures 5c and 5d. The region with observed BT > 239 K represents the spurious non-precipitating cloud pixels whereas the region below 239 K denotes the ABI observed non-precipitating clouds. For both regions, the online bias correction exhibits less systematic deviation compared to the offline approach (dashed curves in Figures 5c and 5d). Similar to the ABI plots, the probability distribution of firstguess deviation from radar reflectivity observation is shown in Figures 5e and 5f. The probability plot of radar reflectivity shows that in an online bias correction experiment more percentage of observations have bias close to zero compared to the offline method. In an offline bias correction experiment, more percentage of observations are negatively biased. The systematic deviation (shown as dotted line) is seen to be larger for strong reflectivity observations. The online approach results in less systematic deviation compared to the offline method. In summary, results in Figure 5 are consistent with those in Figure 4.
The analysis BT plots from 1750 to 1900 UTC reveals visual improvements of the online bias correction approach in comparison to the offline bias correction approach. Figure 1 is used as a qualitative verifying reference. Consistent with Figures 4c and 4d, the analyzed ABI radiances by the online approach over the spurious cloud regions show major improvement compared to the offline approach from 1800 to 1820 UTC. In addition, the rapid development of the supercell storms is better captured by the analysis of the online bias correction approach (Figures 6l-6p vs. Figures 6t-6x vs. Figures 1d-1h). Specifically, in comparison to the offline approach, the online bias correction reduces the delay in capturing the deepening phases of the storms. At 1900 UTC when the development of supercells reaches the mature stage and the sawtooth plots reach their stability (Figures 1 and 4), the offline approach is not able to develop higher level clouds which are responsible for the cold BT observations. The visual evaluation in Figure 6 is consistent with the statistical results in Figures 4 and 5. In summary, the results in Figures 4-6 suggest that the online bias approach improves the analysis and first guess of clouds during the suppression of spurious cloud and the development of the supercell storms. The online bias correction approach not only directly improves the radar anchored region but also indirectly improves the spurious and observed non-precipitating cloud regions.

Impact on the Forecast
Apart from the analysis differences, performance of deterministic forecasts from offline and online bias correction experiments are evaluated. For quantitative comparison of the forecasts, fraction skill score is calculated with 30 km neighborhood radius and 35 dBz composite reflectivity threshold (Mecikalski & Bedka, 2006) over the limited domain. It has been shown by Johnson et al., 2021 that useful skill can be realised starting from the smallest radius of 15 km. For this study, FSS values were calculated with both the 15 (not shown) and 30 km neighborhood radii. The relative qualitative performance of the offline and online bias correction experiments were found to be the same for both radii. One hour forecasts initialized every 10 min from the ensemble mean EnKF analyses valid at 1820 UTC to 1900 UTC are verified against the gridded Multi Radar Multi Sensor composite reflectivity observations (Zhang, Howard, et al., 2016).
The qualitative and quantitative verification shows that the forecast from the two bias correction experiments show more skillful forecasts compared to the experiments without bias correction (Figures 7 and 8). Online bias correction experiment is consistently better than the offline bias correction at all initialization times (Figures 7 and 8). For example, for the forecasts initialized from 1820 UTC (Figures 7k and 7p and red lines in Figure 8), compared to the offline bias correction, the online bias correction is able to show the first hints of the northern storm although the simulated storm is short lived compared to reality (Figure 7a). The southern storm in the online approach is seen strengthening and much longer lived compared to the offline approach (Figures 7k and 7p). Consistently, the FSS of the offline approach has a low value around and below 0.2 for less than 10-20 min and drops to zero afterward (red lines in Figure 8). In comparison, the FSS of the online approach starts with a skill of ∼0.5, increases to a value of ∼0.6 before dropping to a value of 0.2 around 1900 UTC. For 1830 UTC initialized forecasts (Figures 7l and 7q), the online approach shows strengthening swaths of reflectivity corresponding to both the northern and southern storms in the first 30 min before the simulated storms start to dissipate. On the other hand, offline approach only shows the rapid strengthening of the southern storm for about 30 min with the northern simulated storm short-lived (Figures 7l and 7q). This subjective evaluation is consistent with the FSS (black lines in Figure 8). Although the skills from both experiments decay rapidly in the first ∼30 min, the online approach still results in superior skill. The forecast from the online approach is able to provide skillful forecasts (FSS > 0.5) in the first 20 min compared to the offline approach which can provide FSS > 0.5 for only 10 min. For the 1840 UTC and 1850 UTC initialized forecasts, both storms are seen strengthening very rapidly and sustained longer than the earlier initialization times. The online approach is able to sustain stronger reflectivity swaths over longer forecast lead times compared to the offline approach especially for the northern storm. The resultant FSS (green and brown lines in Figure 8) shows consistently higher skill from the online approach than the offline approach. For the 1850 UTC forecasts (brown lines in Figure 8), the skill from the offline experiment decreases rapidly in the first 20 min (till 1910 UTC) after which the FSS drops below 0.5. The online approach is able to sustain the useful forecast until about 1930 UTC. For the 1900 UTC initialized forecast, the online approach is seen to develop both the northern and southern storms with the track, longevity and morphology more similar to the reality than the offline approach. In the offline approach the north storm weakens after about 30 min and redevelops which is seen as a small break in the corresponding reflectivity swath. Due to enhanced convection further north of two long lived supercells, several new severe thunderstorms initiate at around 1900 UTC. Although the rapid development of the thunderstorm immediately north of the northern supercell is seen in both the online and offline bias correction experiments, the online approach captures the storm near the Oklahoma Panhandle better than the offline approach. The FSS of the online approach starts with a skill of about 0.8 against 0.7 from the offline bias correction experiment (purple curves of Figure 8). For about 50 min into the forecast, the online approach is able to maintain its skill above 0.5, whereas the offline approach could not maintain its skillful forecast for more than 40 min. Verification of forecast using other reflectivity thresholds and neighborhood radii (not shown) revealed similar improvement in forecast skills gained by the online bias correction approach. Verification was also performed for the probabilistic forecasts derived from the 40 ensemble members and the results (not shown) are consistent and support those of the deterministic forecast shown in Figures 7 and 8.

Physical Diagnostics
In this section, physical diagnostics are performed during the DA cycling to examine the reasons behind the improved analysis by the online approach during the development stage of the observed supercells. The difference in the bias correction approach has very little effect in the first two cycles. The background BT departure valid at 1800 UTC from both offline and online bias correction experiments show similarly large negative innovations due to the absence of clouds in the ensemble mean background over the observed storm location (not shown). Beginning at the 1800 UTC analysis cycle, the offline bias correction approach removes a large part of the innovation as bias as shown by a smaller bias corrected innovation (Figure 9a). However, the online approach retains most of the negative innovation over the observed storm region (Figure 9c). The differences of the bias corrected innovations are more apparent during later cycles as the supercell continues to develop. At 1810 UTC, the online approach again estimates a smaller bias and therefore retains more negative innovations over the Figure 9. Bias corrected radiance innovation in K from offline (a), (b) and online (c), (d) bias correction experiment at 1800 and 1810 UTC also shown is the plane AB which is used to plot cross-section along the storm.
northern and southern storm locations compared to the offline approach (Figures 9b and 9d). The difference in the estimated bias from the offline and the online approach is consistent with the theoretical expectation presented in Section 2.1. The estimated bias in the offline approach includes a larger contribution from the model bias, which in this case is the underestimation of storm coverage by the model.
The different impacts of the online and offline bias correction are further understood by the resulting analysis increments of humidity, hydrometeor mixing ratios and updraft velocity. At 1800 UTC, the model background has no hydrometeors in the region of observed storms. A cross-section (line AB in Figure 9 for 1800 UTC) of background ensemble spread of RH shows similar background ensemble spreads in offline and online bias correction experiments (Figures 10a and 10b). Given the larger negative bias corrected innovation in the online approach (Figure 9c), the online bias correction experiment produces larger analysis increments to RH at the observed northern and southern storm locations (which is denoted by red and blue circles in Figure 10). The resultant differences of the RH analyses influence the hydrometeor field in the subsequent background forecast at 1810 UTC. As shown in Figure 11, the background hydrometeor ensemble spread at 1810 UTC in the online approach is larger than the offline approach especially at the observed southern storm location and the upper portion of the northern storm. Such a larger background hydrometeor spread together with the greater bias corrected innovation results in a larger positive hydrometeor increment at both lower and higher model levels for the southern storm and at upper levels for the northern storm by the online approach. At 1820 UTC, the hydrometeor spread and increment continue to grow. The online approach shows again larger background hydrometeor spread and greater positive increments at the observed northern and southern storm locations. The greater addition of the positive hydrometeor increments lead to rapid formation of the supercells with a coherent storm structure in the online approach, consistent with Figure 6. At 1830 UTC, the online approach is able to spin up both northern and southern storms with a deep towering cloud structure with additionally greater increments of the hydrometeors. In summary, Figures 10 and 11 show that the online approach is able to more appropriately increment the moisture and hydrometeor fields during the DA cycling. As a result, the analysis of the online approach is able to capture the developing storms more rapidly, leading to the reduction of conditional bias over observed cold BT regions (Figure 5b). Such a difference is rooted back to the fundamental differences of the two approaches in estimating observation biases.
In addition to the thermodynamic field, the online and offline bias correction experiments analyze dynamical fields differently. In view of this, the background updraft velocity at 500 hPa from offline and online bias correction experiment at 1810, 1820 and 1830 UTC is shown in Figure 12. At 1810 UTC, there is no significant difference in background updraft velocity between the offline and online bias correction experiment. As the storm intensifies, at 1820 UTC, the online bias correction experiment shows stronger ensemble mean background updraft (represented by red color in the plot) at northern and southern storm location compared to the offline bias correction experiment. Similar differences in updraft velocity can be observed at 1830 UTC. This stronger updraft along with the more proper hydrometeor and moisture field analysis by the online approach contribute to the superior forecast of the two supercells shown in Figure 8. In summary, the online bias correction retains useful information in the innovation, which in turn produces different analyses and background ensemble spread for both thermodynamic and dynamical fields. The effect is accumulated during the DA cycling that is responsible for the superior analysis and forecast of the supercells. online bias correction experiment at 1800 UTC with ensemble mean analysis increments shown as contour lines from −30% to 30% in 3% intervals. The positive and negative increments are shown in solid and dotted lines respectively. The x-axis denotes the distance in km measured along the plane AB. The locations of northern and southern storm are shown in red and blue colored circles respectively. Figure 11. Cross-section (along AB as in Figure 9) plot of background ensemble spread in total (sum of solid and liquid phase) hydrometeor (in g/kg) with contour lines (positive increments in solid line and negative increments in dotted line) representing the ensemble mean analysis increments to hydrometeor at 1810, 1820 and 1830UTC.

Conclusions
This study implemented an online nonlinear bias correction for all-sky ABI radiance DA with a rapidly updated EnKF for a case of rapidly developing supercells. The online nonlinear bias correction was comprehensively compared with the corresponding offline approach through examining the analysis and prediction of a rapidly developing supercell event of 18 May 2017 over the U.S. Great Plains.
The use of the radar reflectivity as the anchoring observation is developed and explored. Although the radar reflectivity and ABI radiances are sensitive to different variables and contribute to the analysis increment in different regions of the cloud, radar reflectivity and ABI cloud sensitive radiances can be physically connected. In this study, the observed and spurious precipitating cloud regions are identified as the regions utilizing the anchoring information from radar. The bias correction coefficients for both the nonlinear offline and online bias correction experiments were calculated and applied using ABI radiances in the vicinity of radar anchored regions. A set of experiments were first conducted to determine the best predictor for the bias. Three predictors, namely observed BT, simulated BT and symmetric BT were examined for the offline bias correction. The simulated BT predictor was found to produce the least innovation bias and RMSI and was therefore chosen for further study. Objective and subjective evaluation were performed during the rapid DA cycling over the radar anchored region, ABI spurious cloud, ABI observed cloud and non-precipitating cloud regions. The online bias correction performs better than the offline approach during the suppression of the spurious clouds and the development of non-precipitating and precipitating cloud regions when the supercell storms develops in reality. Apart from the direct improvement over the radar anchored region, the online bias correction was able to provide a firstguess with less bias and RMSI than the offline approach over observed and spurious non-precipitating cloud regions. In addition to evaluating the analyses and the background forecasts during the DA cycling, quantitative and subjective verification of the deterministic forecasts showed consistent superior performance from the online bias correction over the offline approach. To further the understanding of the effect of various bias correction methods, several physical diagnostics were performed during the development of the supercells. Consistent with the online bias correction theory, the estimated bias in the offline approach included a larger contribution from the model bias, which in this case is the underestimation of storm by the model. Whereas, the online bias correction retains useful information in the innovation, which in turn produces different subsequent analysis, background and background ensemble spread for both the thermodynamic and dynamic fields. The effect is accumulated during the DA cycling that is responsible for the superior analyses and forecast of the supercells.
This study is the first to explore the use of radar reflectivity as the anchoring for online nonlinear bias correction to assimilate GOES-R ABI all-sky infrared radiances for the analysis and prediction of rapidly developing supercells. As a proof of concept, the present study utilizes only channel 9 water-vapor sensitive radiances from ABI. Given this is the first time such a study is presented and the need to introduce the details of the approach and the diagnostics, the present study focuses on a single high impact weather event for the demonstration. Further studies should evaluate the effect of the online and offline bias correction assimilating additional water-vapor and surface sensitive channels and extend to multiple cases. Furthermore, work is underway to implement the nonlinear online and offline bias correction approach in EnVar.

Data Availability Statement
All data required for running the experiments and produced by the experiments are archived locally and available upon request. The GSI EnKF source code can be found from https://dtcenter.org/community-code/ gridpoint-statistical-interpolation-gsi.