Impact of Model Bias Correction on a Hybrid Data Assimilation System

Hybrid data assimilation combines a conventional 3-D or 4-D variational system with background error covariance (BEC) generated from ensemble forecast systems. In order to achieve better BEC, three perturbation schemes, namely, the random combination of multiple physical paramterization schemes (referred to as MP), the MP plus stochastical perturbation on physical process tendencies (MP-SPPT), and the unified perturbation of stochastic physics with bias correction (UPSB, proposed by the authors of this paper in a previous work), were first used in a regional ensemble model, i.e., the Global and Regional Assimilation and Prediction System-Regional Ensemble Prediction System (GRAPES-REPS), and the BECs thus obtained were compared for 7-day ensemble forecasts. The results show that UPSB, which is in fact an MP-SPPT but with the systematic model bias removed, has a better consistency, i.e., the ratio between root-mean-square error (RMSE) and ensemble spread is much closer to 1, especially at low model levels, compared to the other two schemes. Moreover, the BEC derived from UPSB captured more reasonable distributions of forecast errors. Second, performance of a hybrid data assimilation system (the GRAPES-MESO hybrid En-3DVar) was evaluated by using the BECs from the three perturbation schemes for 7-day hybrid data assimilation forecasts, and thus disclosing the effect of the model bias correction (assuming that the random stocastical features are in general offset in the three perturbation schemes) on the hybrid system forecasts. A covariance weight of 0.8 was prescribed, and this value was determined through sensitivity experiments. The forecast results from the hybrid data assimilation system show that UPSB reduced the false correlation between distant points. The quality of analysis fields of the UPSB scheme shows visible improvement, i.e., the analysis fields produced by UPSB have much smaller RMSEs than those of the other two schemes, at all vertical model levels. The quality of the hybrid data assimilation forecast fields was also improved by this scheme. Furthermore, the improvement was much greater in the early stage of the assimilation cycle than in the late stage. Generally, the quality of the hybrid data assimilation of GRAPES-MESO hybrid En-3DVar could be efficiently improved by the model bias correction in the UPSB scheme.


Introduction
Data assimilation is a method typically used to seek the optimal combination of background (usually in the form of numerical model results) and newly received observations. It offers updated and accurate initial condi-tions for a numerical model to make better forecasts. A new trend in the development of data assimilation schemes in recent years is the so-called hybrid approach, which attempts to achieve more accurate and flow-dependent ensemble estimates of the background error covariance (BEC), and which can then be combined with statistical estimates of static BEC (Buehner, 2005;Wang, 2011). Hybrid data assimilation combines a conventional 3-D or 4-D variational system with BEC generated from ensemble forecast systems. Hamill and Snyder (2000) used a linear method to introduce ensemble background error statistics to a 3DVar (three-dimensional variational) algorithm to build a hybrid scheme, and then they preliminary tested it in a simple model. Lorenc (2003) introduced an ensemble-estimated BEC matrix to the 3DVar system via the Extending Control Variable (ECV) method and proved its advantages. It was also confirmed that the linear combination method is equivalent to the ECV method (Wang et al., 2008a). Wang et al. (2008a, b) developed an ETKF-3DVar hybrid data assimilation approach (ETKF; ensemble transform Kalman filter) by putting the ETKF into the 3DVar system, and then applied it to the Weather Research and Forecasting (WRF) data assimilation system. In addition, many scientists have conducted extensive research and experiments on hybrid data assimilation schemes based on different models (Liu et al., 2009;Zhang et al., 2009;Buehner et al., 2010a, b;Wang et al., 2013;Ma et al., 2014;Chen et al., 2015;Zhang et al., 2015). It can be concluded that hybrid data assimilation schemes perform very well in global and regional numerical weather prediction systems, especially for model levels in the troposphere and in areas with relatively poorer observational coverage (Wang, 2011;Liu and Xiao, 2013;Xia et al., 2018).
As mentioned above, the hybrid method applies flowdependent, ensemble-based estimates of BEC to the 3DVar/4DVar system, which allows the improvement of the ensemble perturbation schemes to benefit the data assimilation skill. Therefore, a choice of perturbation scheme is crucial for a better ensemble BEC, which may have positive effects on ensemble estimates of BEC. Currently, there are many ensemble perturbation schemes, such as the initial perturbation methods of singular vector schemes (Buizza and Palmer, 1995;Molteni et al., 1996), ETKF (Bishop et al., 2001;Wang and Bishop, 2003;Wei et al., 2006;Ma et al., 2008), random Monte Carlo perturbations (Hollingsworth, 1980;Mullen and Baumhefner, 1994), and the breeding method Kalnay, 1993, 1997), as well as the model perturbation methods of Stochastically Perturbed Parameterization Tendencies (SPPT), which can also be interpreted as stochastical perturbation on physical process tendencies (SPPT), and Stochastic Kinetic Energy Backscatter (Shutts, 2005). At present, SPPT is one of the most popular approaches, and has been successfully applied in global and regional hybrid assimilation system experiments. The basis of the SPPT scheme is to randomly perturb the tendencies from the physical parameterization schemes, which aims to represent uncertainties in the effects of the sub-grid physics processes that the atmospheric physics parameterization schemes are designed to describe. Berner et al. (2009), Charron et al. (2009), Tennant et al. (2010, and Yuan et al. (2016) showed that the ensemble perturbation schemes can increase the ensemble spread and reduce the impact of systematic errors of the model to some extent. Toth et al. (2003) and Wang et al. (2018) found that systematic errors can induce an inaccurate ensemble distribution, which renders the ensemble BEC unable to accurately represent the forecast error. This may negatively impact the quality of hybrid data assimilation. The main approach to dealing with systematic errors is statistical post-processing, via methods such as Ensemble Model Output Statistics (EMOS; Gneiting et al., 2005), Analog Bias Correction (Ren and Chou, 2005), Bayesian Model Averaging (BMA; Raftery et al., 2005;Hamill et al., 2011), Kalman filter predictor bias correction (Monache et al., 2006), the trade-off in bias correction method (Cui et al., 2006;Du, 2007), statistical downscaling (Wilby and Wigley, 1997;Wang and Zhi, 2015), and so on. These methods are effective at minimizing the impact of systematic errors on the ensemble predictions, but cannot deal with random errors. Considering these model shortcomings, Chen et al. (2019) introduced a method to remove the systematic bias of each ensemble member in the model integration process. With this in mind, might it be possible to consider the impact of the model systematic errors together with stochastic errors in SPPT scheme to optimize the performance of the ensemble forecast system and improve the quality of BEC estimation in the hybrid assimilation cycle? This idea is worthy of investigation.
In this study, we tested this hybrid implementation idea by using unified perturbation of stochastic physics with bias correction (UPSB, or Exp3 in Table 1; also see Xia et al., 2019) with the GRAPES-REPS (Global and Regional Assimilation and Prediction Enhanced System-Regional Ensemble Prediction System) and GRAPES-MESO hybrid En-3DVar-TD-HLS (Ensemble three-dimensional hybrid data assimilation for GRAPES with Topographic Dependent Horizontal Localization Scale scheme), and compared the results with the traditional multi-physical process scheme and multi-physical process with SPPT-only methods. We adopted these three ensemble perturbation schemes to obtain background statistics coupled with the static background covariance to investigate the impact of model bias correction on our hybrid data assimilation system. In Section 2, the GRAPES-REPS system, the GRAPES-MESO hybrid En-3DVar-TD-HLS system, and the design of our experi-

Introduction to the GRAPES hybrid data assimilation system
The GRAPES-MESO hybrid En-3DVar-TD-HLS system includes the GRAPES-MESO 3DVar data assimilation system and the GRAPES-REPS regional ensemble forecast system. The GRAPES-REPS system has 15 ensemble members, which includes 14 ensemble perturbation forecasts and 1 control forecast. The horizontal resolution is 0.15° × 0.15°, with 502 × 332 grids on 49 vertical model levels. The domain of the system is 15°-64.35°N, 70°-145.15°E, which covers conterminous China. Forecasts are made available every 6 hours. The valid forecast time is 72 h, with an integration time of 60 s. The background and boundary conditions are provided by the Numerical Weather Prediction Center of the China Meteorological Administration. The horizontal resolution of the data assimilation system is the same as the ensemble forecast system. The hybrid data assimilation combines the ensemble BEC with the 3DVar static BEC via ECV. The cost function of the hybrid data assimilation system is as follows: where x is the analysis increment, is the weight of the climatological and statistical BEC, is the ensemble estimated BEC, and . In our study, is 0.2 and is 0.8 (based on previous research on the GRAPES-MESO hybrid En-3DVar-TD-HLS system). H is the observational operator, R is the observational error covariance matrix, and T represents the matrix transpose.
The ensemble BEC matrix of GRAPES-MESO hybrid En-3DVar-TD-HLS is composed of the 12-h forecasts of GRAPES-REPS with ensemble perturbations. The method is as follows: Here, K=14, which indicates the total number of ensemble members, is the ensemble mean of those members, x k is the forecast of K members, and denotes the ensemble BEC matrix with flow-dependence.
In Eq. (2), the estimated BEC matrix is not a full rank matrix, which brings about large sampling errors of estimated covariance. Thus, it is necessary to use a localized correlation matrix C to resolve this problem, and C and must have the same dimension size. Then, the localization of can be done by multiplying them together: • where the symbol " " is the Schur product. In this paper, the horizontal localization scale of the static BEC matrix is 500 km (default value of the GRAPES-MESO 3DVar data assimilation system), while the ensemble BEC matrix adopts the Topography Dependent Horizontal Localization (TDHL) scale scheme (the horizontal localization only changes with the horizontal direction) [refer to Xia et al. (2019)]. Figure 1 is the horizontal distribution of the TDHL scale scheme over different terrain heights, from which we can see that the horizontal localization scale of the BEC is proportional to the terrain height. In other words, the higher the terrain, the larger the horizontal localization scale. The maximum horizontal localization scale of the system is 1500 km. The vertical localization scale of both the static and ensemble BEC is 5 km (default value of the system). Besides, the BEC information on the humidity is completely static in the GRAPES-MESO En-3DVar hybrid assimilation system, due to the high uncertainty of water vapor in the GRAPES-REPS system.

Method
The GRAPES-REPS regional ensemble prediction system uses the downscaling initial perturbation method and model perturbation method of the MP scheme. The latter includes two planetary boundary layer (PBL) parameterization schemes and four cumulus convective (CC) MP + stochastical perturbation on physical process tendencies (SPPT) As above As above Exp3 MP + unified perturbation of stochastic physics with bias correction (UPSB) As above As above 402

Journal of Meteorological Research
Volume 34 parameterization schemes. Then, it randomly combines two of them ( Table 2). The PBL parameterization schemes we used were the MRF (Medium Range Forecast) scheme (Hong and Pan, 1996) and YSU (Yonsei University) scheme (Hong et al., 2006). The CC parameterization schemes were the shallow convection Kain-Fritsch-Eta scheme (Kain and Fritsch, 1993;Kain, 2004), Betts-Miller-Janjic scheme (Betts, 1986), simplified Arakawa-Schubert scheme (Pan and Wu, 1995), and original Kain-Fritsch scheme (Kain and Fritsch, 1990). The integration formula in the GRAPES-REPS regional ensemble prediction model is where is the integration result of the jth ensemble member from t 0 to t, where ; t 0 is the start time, and t is the integration duration from t 0 ; j = 0 represents the control forecast; A is the integration tendency term for the dynamical process; and P describes the physical process.
The SPPT scheme was incorporated into the GRAPES-REPS system based on Exp1 (Table 2), and its integration formula is: where represents the random perturbation of the jth ensemble member. The other variables carry the same meaning as in Eq. (5). For the jth ensemble member, the expression of can be obtained by: where is the random perturbation with no boundary; is the mean value of the random perturbations; is the spectral coefficient, which changes with time t; λ and ϕ are longitude and latitude; is the spherical harmonics, where l and m respectively stand for the horizontal total wave number and the latitudinal wave number; and L denotes the horizontal truncation of the random perturbations. The evolution of is based on a Markov chain, which is shown in Eq. (8): ∆t where is the specified time interval (in this study, 60 s corresponds to the integration length of the GRAPES-REPS system), τ is the timescale of the scalar product for the random field, meets a Gaussian distribution with a variance of 1 and a mean value of 0, and σ is the standard deviation of .
Based on the random field defined by Eqs. (7) and (8), a stretching function is introduced to generate a that can set the range of variations (given the upper and lower boundaries) and can change the PDF distribution: (9) In this paper, is the same as in Yuan et al. (2016) : where , in which and are the upper and lower boundaries of . A recent study showed that the performance and verification of the GRAPES-REPS ensemble is highly sensitive to model bias, and Chen et al. (2019) and Xia et al. (2019) introduced a new method to remove systematic bias of each ensemble member in the model integration process. Based on SPPT and this new method, a UPSB scheme was designed to correct the linear tendency bias (based on the control forecast from GRAPES-REPS) of potential temperature in the procedure of tendency integration for each ensemble member. The integration formula is: where is the tendency of systematic bias at each grid in three dimensions for each integration length, and can be obtained as follows: where Δ is the valid forecast time, is the bias of the atmospheric state at time t, and is the bias at t + Δ. We define the difference between the forecast and the "truth" as the model systematic bias, wherein the analysis fields from the T639 global forecast system are used as the "truth" in this paper. Here, δt is the integration length (60 s in this paper).
Through statistical analysis of the systematic error of various variables of the control forecast in the GRAPES-REPS system, it is found that the wind and pressure fields have nonsignificant and nonlinear systematic errors. The magnitude of their model error is also small. The following is the integration formula of potential tem-perature with linear tendency bias correction: where represents the linear tendency bias of the control forecast.
In this paper, the linear tendency bias of the potential temperature has the same value for each integration length during the entire valid forecast time. For a detailed description of the method, such as how to calculate the bias tendency and how it impacts the forecasts, readers are referred to Chen et al. (2019) and Xia et al. (2019).

Experimental design
Seven-day ensemble perturbation experiments of the three schemes were proposed to explore their impact on the hybrid data assimilation cycle with the dynamicaldownscaled background and boundary conditions from the T639 global ensemble prediction system. The ensemble experiments were from 1200 UTC 5 July to 1200 UTC 11 July 2015, with 72-h forecasts for each ensemble forecast experiment, and the 12-h forecasts were used to calculate the spread, ensemble RMSE (rootmean-square error), ensemble mean forecast error, correlation coefficient between spread and ensemble mean forecast error, and the ensemble BEC. The design of the ensemble experiments can be found in Table 1. In addition, sensitivity experiments with varied ensemble covariance weight were conducted via the 7-day hybrid data assimilation experiments (Table 3) to obtain the best weight to be used in BEC coupling. Then, 7-day hybrid data assimilation experiments were conducted based on the ensemble BEC generated by the 7-day ensemble perturbation experiments. The hybrid data assimilation experiments were conducted from 0000 UTC 6 July to 0000 UTC 12 July 2015, with 72-h forecasts for each experiment. Figure 2 illustrates the design of the experiments. For example, if the experiment time of the hybrid data assimilation was 0000 UTC 12 July 2015, we needed to start the ensemble experiments at 1200 UTC 11 July 2015 to obtain 12-h forecasts, which were used to form the ensemble BEC.
The hybrid data assimilation experiments were only applied to conventional observations (Fig. 3). The variables U, V, pressure, and relative humidity were contained in the surface observations at land stations (SYN-OP), radiosondes from land and ships (TEMP), aircraft reports (AIREP), and surface observations on ships

404
Journal of Meteorological Research Volume 34 (SHIPS) and were assimilated by the GRAPES-MESO hybrid En-3DVar system; only U and V were assimilated by surface observations at land stations (SATOB). The analysis fields from the T639 global forecast system (resolution: 0.2815 × 0.2815) were interpolated to the regional study area with a resolution of 0.15 × 0.15. This could be used to evaluate the bias of the control forecast (72 h) of GRAPES-REPS, the result of the ensemble mean RMSE, and the results produced by the hybrid data assimilation system. The model bias was estimated from the difference between the control forecasts and the T639 global analysis field. For example, if we want to estimate the model bias at around 1200 UTC 11 July 2015, 10day forecasts from 1200 UTC 28 June 2015 to 1200 UTC 8 July 2015 were conducted to estimate the model bias, and then we calculate the average model bias of these 10 days.

Ensemble spread and RMSE
The ensemble spread and ensemble mean RMSE are important indicators to evaluate the performance of the ensemble prediction system. The best match between them is when their ratio (consistency) becomes 1. Figure 4 shows the vertical profile of the consistency between ensemble spread and mean RMSE over 12-h forecasts. It can be seen that Exp1 had worse consistency than Exp2 and Exp3 for U at all levels. Also, it is clear that Exp3 performed best in the match of ensemble spread and mean RMSE from model levels 1 to 20. For other levels, there was only a small difference between Exp2 and Exp3, and Exp3 was slightly better than Exp2. The consistency of V had similar vertical structures as U. The advantages of Exp3 for pressure (PI) were visible at model levels 1-22. On the whole, Exp3 had the best consistency between ensemble spread and mean RMSE among the three perturbation ensemble experiments.

Correlation between spread and forecast error
The ensemble estimates of the BEC matrix were coupled with the static BEC matrix in the hybrid data assimilation cycle. The coupled BEC matrix can be used to map the forecast error with respect to the relevant atmospheric state, meaning that it should be flow-dependence. Thus, it is necessary to examine the correlation between the ensemble spread and the forecast error. Figure 5 shows vertical profiles of the correlation coefficient between the ensemble spread and the ensemble mean forecast error for 12-h ensemble forecasts (the mean of 7-day experiments) with the three ensemble

.2 Correlation between a single point and grid points of 975 hPa
In terms of current computing costs, the ensemble BEC can only be estimated with a limited set of members. This may result in a false correlation between longdistance grid points, thereby overestimating the ensemble BEC and weakening the impact of observations on the assimilation system. The correlation coefficients of ensemble spread between a single point (6th layer, 34°N, 123°E) and the grid points of the 6th layer for U and V are shown in Fig. 6. The single point selected was the center of a vortex. The correlation coefficients of both U and V were relatively large near the vortex, basically in the range of 0.8-1. Also, its distribution was closely connected to that of the vortex. This indicates that the combination of the ensemble BEC with the static Vertical profiles of the ratio of ensemble spread to mean RMSE (consistency) over 12-h forecasts with the three ensemble perturbation schemes, where the green line is for Exp1, blue for Exp2, red for Exp3, and gray is the "perfect" line, i.e., the optimal ratio for (a) U, (b) V, and (c) pressure (dimensionless quantity, PI).

Journal of Meteorological Research
Volume 34 BEC can make the hybrid assimilation system have an obvious flow-dependent characteristic in the study area. It also means that the BEC generated by the three ensemble perturbation schemes can reasonably represent the analysis error with respect to the weather changes.
Although the distributions of the single point correlation coefficients of the three experiments had similar patterns near the vortex, the coefficients in the red box, which is far away from the vortex, had some differences. Exp3 had relatively larger coefficients in the red box, and this may have been due to the limited ensemble size. Compared with Exp1, Exp2 and Exp3 had smaller coefficients in this area. All the above results show that UPSB can moderately reduce the false correlation. Table 4 lists the results of the ensemble covariance weight (7-day average) sensitivity experiments. It is evident that, when the ensemble covariance weight is 0.8, the U RMSE of the analysis, 6-, and 12-h forecast fields is smaller than the other ensemble covariance weights. Besides, the V RMSE is smaller in the analysis and forecast fields when the ensemble covariance weight is 0.8. Therefore, we defined the optimal ensemble covariance weight as 0.8 in our hybrid data assimilation system and used it in later experiments. Figure 7 shows vertical profiles of the analysis field RMSE (mean of 7-day experiments) from the hybrid data assimilation with the three ensemble perturbation schemes. The wind and temperature RMSEs after the hybrid assimilation of Exp3 were much smaller than in Exp2 and Exp1 at all model levels. This suggests that Epx3 improved the quality of the analysis fields pro-duced by the hybrid assimilation system. The maximum difference for U between Exp3 and Exp1 (Exp2) was 0.18 (0.2) m s −1 , respectively. The maximum difference for V, meanwhile, was 0.15 (0.17) m s −1 , respectively, and for temperature, it was 0.08 (0.1) K. The improvement in the performance of the UPSB scheme in hybrid data assimilation was much better at lower model levels than at high levels.

RMSE analysis
The vertical distribution of the RMSE of the 6-h (mean of 7-day experiments) forecast fields is shown in Fig. 8, revealing that the RMSE of U from Exp3 was smaller than that from Exp1 and Exp2 below 200 hPa, but similar to Exp1 and Exp2 above 200 hPa. Also, the UPSB scheme featured less improvement for the 6-h wind forecasts when compared with the improvement for the analysis field. Above 500 hPa, Exp3 had a smaller RMSE of temperature, while below this level the performances of Exp1, Exp2, and Exp3 were almost the same.
The vertical profiles of the RMSE of the 12-h forecast fields (mean of 7-day experiments) were provided in Fig.  9. Below 400 hPa, Exp3 had a smaller RMSE of U, while above that the performances of Exp1, Exp2, and Exp3 were almost the same. The three experiments had small differences across the vertical model levels for temperature field.
In order to check the improvement of Exp3, the ratio of its RMSE to Exp1 and Exp2 is given in Fig. 10 (mean of 7-day experiments). Also, the percentage is calculated as: where u stands for U-wind.
The results show that Exp3 performed beneficially for both the wind and the temperature fields in the hybrid data assimilation cycle. Also, improvement of the analysis APRIL 2020 fields was visible, at 2.8% and 4.4% compared with Exp1 and Exp2, respectively, for U-wind, and 3.5% and 4.6% for V-wind. The improvement ratio of Exp3 for temperature reached 4.4% and 5.3% compared to Exp1 and Exp2, respectively. However, this benefit reduced with the length of the valid forecast time.
Overall, compared with Exp2 and Exp1, Exp3 featured smaller RMSE, especially for the analysis fields and in the early forecast in the hybrid data assimilation cycle. Table 4. Total RMSE of U and V in the analysis, 6-, 12-, and 24-h forecast fields produced by the ensemble covariance weight sensitivity experiments. The ensemble covariance weights were 0.0, 0.2, 0.5, 0.8, and 1.0, and the bold number is the least RMSE of these five experiments

Conclusions and discussion
Three model perturbation methods, which are unified perturbation of stochastic-physics with model bias removed (UPSB), random combination of multiple physical parameterization schemes (referred to as MP), the MP plus stochastical perturbation on physical process tenden-cies (MP-SPPT), were used in GRAPES-REPS and GRAPES-MESO hybrid En-3DVar system to investigate the impact of ensemble model systematic bias on ensemble BEC and hybrid data assimilation system forecast results. Seven-day ensemble forecasts (1200 UTC 5 July to 1200 UTC 11 July 2015) were conducted by using the above three perturbation schemes to analyze the   Fig. 7, but for the 6-h forecast field RMSE from hybrid data assimilation.  V (m s −1 ) Temperature (K) Fig. 9. As in Fig. 7, but for the 12-h forecast field RMSE produced by hybrid data assimilation.
characteristics of the ensemble BEC. The ensemble BEC was coupled with static BEC via ECV for the GRAPES-MESO hybrid En-3DVar system, and then, 7-day hybrid data assimilation forecasts (0000 UTC 6 July to 0000 UTC 12 July 2015) were conducted. The results can be summarized as follows.
(1) The correlation between the BEC constructed by the UPSB scheme and the ensemble forecast error of the 12-h prediction is better than the other two schemes at the middle and lower model levels for the wind field. The correlation coefficients of the UPSB scheme for PI were higher than those of the other two schemes at model level 1-35.
(2) In the single point correlation at the 6th layer of the model, all the ensemble BEC matrices generated by the three ensemble perturbation schemes were flow dependent and could accurately represent the weather states near the vortex area. Nonetheless, the UPSB scheme and SPPT scheme can reduce the possibility of false correlations caused by small sample sizes.
(3) UPSB also demonstrated improvement in the hybrid data assimilation cycle for both the analysis and forecast fields, based on RMSE evaluations. However, this benefit was reduced by the length of the valid forecast time.
In general, the UPSB scheme can improve the quality of ensemble perturbations, thereby enabling the construction of a more accurate ensemble BEC matrix that better represents the characteristics of the forecast error, and thus significantly improving the quality of the product produced by the hybrid data assimilation system. Cui, B., Z. Toth, Y. J. Zhu, et al., 2006: The trade-off in bias correction between using the latest analysis/modeling system with a short, vs. an older system with a long archive. Proceedings of the First THORPEX International Science Symposium,

Journal of Meteorological Research
Volume 34