Length Scale Analyses of Background Error Covariances for EnKF and EnSRF Data Assimilation

: Data assimilation (DA) combines incomplete background values obtained via chemical transport model predictions with observational information. Several 3-Dimensional variational (3DVAR) and sequential methods (e.g., ensemble Kalman ﬁlter (EnKF)) are used to deﬁne model errors and build a background error covariance (BEC) and are important factors affecting the prediction performance of DA. The BEC determines the spatial range, where observation concentration is reﬂected in the model when DA is applied to an air pollution transport model. However, studies investigating the characteristics of BEC using air quality models remain lacking. In this study, horizontal length scale (HLS) and vertical length scale (VLS) analyses of a BEC were applied to EnKF and ensemble square root ﬁlter (EnSRF), respectively, and two ensemble-based DA methods were performed; the characteristics were compared with those of a BEC applied to 3DVAR. The results of 6 h PM 2.5 predictions performed for 42 days were evaluated for a control run without DA (CTR), 3DVAR, EnKF, and EnSRF. HLS and VLS respectively exhibited a high correlation with the ground wind speed and with the planetary boundary layer height for diurnal and daily variations; EnKF and EnSRF exhibited superior performances among all the methods. The root mean square errors were 11.9 µ g m − 3 and 11.7 µ g m − 3 for EnKF and EnSRF, respectively, while those for 3DVAR and CTR were 12.6 µ g m − 3 and 18.3 µ g m − 3 , respectively. Thus, we proposed a simple method to ﬁnd a Gaussian function that best described the error correlation of the BEC based on the physical distance.


Introduction
Chemical data assimilation (DA) was proposed for reducing the uncertainty of a chemical transport model (CTM) [1][2][3][4]. Data assimilation combines (mixes) incomplete background values obtained via CTM prediction with observational information, including errors. The obtained results are close to the true values with lower errors compared with the uncertainty of each model and the observation. The DA method is divided into the variational (e.g., three-and four-dimensional variational methods (3DVAR and 4DVAR, respectively)) and sequential methods (e.g., optimal interpolation (OI), ensemble Kalman filter (EnKF)). Recently, efforts were made to increase the predictability of aerosols by applying various DAs to CTM and by utilizing ground-based or satellite observations (OI: [5,6], 3DVAR: [7][8][9], 4DVAR: [10,11], and EnKF: [12,13]). Various DA methods were applied to compensate for the shortcomings of the EnKF, which needs to perturb observation data. For example, the ensemble square root filter (EnSRF) [14,15], local ensemble transform Kalman filter (LETKF) [16,17], and ensemble adjustment Kalman filter (EAKF) [18,19] are representative methods applied to CTMs. All studies on DA suggest that the predictability of aerosols, including PM 2.5 , can be improved if the initial field and model parameters are improved by assimilating the observed data. Methods for defining model errors and building the background error covariance (BEC) are important factors affecting the prediction performance of the DA technology [20].
The easiest approach to building a BEC in a 3DVAR or 4DVAR is using the National Meteorological Centre (NMC) method [21]. The NMC method defines the difference between two different prediction periods that predict the same time targets as climatically determined model errors. In this case, the analysis increment using DA at observation points is distributed in the form of a Gaussian-type circle (isotropic). In contrast, an ensemble-based DA calculates the BEC by assuming that the ensemble average is the true state through short-time ensemble prediction; it defines the degree to which each ensemble member falls from the mean as a model error. Therefore, it is characterized by the time-variant BEC, and it reflects the atmospheric flow at the time of the DA (i.e., it is flow-dependent). Pagowski and Grell [22] compared the DA of the EnSRF to 3DVAR methods to assimilate PM 2.5 ; they pointed out that EnSRF was superior to 3DVAR for a 6 h prediction when examining the advantages. A simple aerosol parameterization method has limitations in predicting PM 2.5 .
Skachko et al. [23] compared the results of applying EnKF and 4DVAR to the CTM simulation of stratospheric ozone; they revealed that EnKF showed comparable performance to 4DVAR and indicated that the predictability of EnKF could be better than that of the 4DVAR for hydrogen chloride (HCl) and nitric acid (HNO 3 ) with a relatively long lifetime. The 4DVAR and EnSRF methods were utilized for the inverse problem estimation of the CO 2 flux considering Chatterjee and Michalak's study [24]. Although observation DA was employed through an ideal experiment, EnSRF provided a more realistic and useful uncertainty assessment; the overall DA performance depended on operational constraints, such as the density and ensemble number of the observational data. Only a few studies have compared each DA method by applying DA to the CTM; there have been very few studies on aerosols.
Recently, Northeast Asia, where aerosols are gaining research interest, has witnessed growth in research on ensemble-based DA. Peng et al. [25,26] suggested that predictability can be enhanced by performing DA on China's ground air pollution observation data using the EnKF method and optimizing the initial conditions (ICs) and emissions of PM 10 , SO 2 , NO 2 , O 3 , and CO, including PM 2.5 , which are the main pollutants. Chu et al. [14] used EnSRF to evaluate the air quality and reduce emissions for the strong emission reduction policy implemented in Beijing in 2015 via modeling; they suggested that more accurate and quantitative contamination concentration and emission estimation were possible. Kong et al. [17] constructed the reanalysis data of air pollutants from 2013 to 2018 with a high resolution (horizontal grid of 15 km × 15 km) using LETKF DA. They argued that these data could be utilized as learning materials for artificial intelligence, health impact assessments of air pollutants, and long-term variability analysis. Recently, Park et al. [13] evaluated the contribution of the improved ICs and boundary conditions (BCs) for the prediction of PM 2.5 by applying EnKF to the community multiscale air quality (CMAQ) model.
Most of the current studies develop and compare each DA method by exclusively focusing on prediction performance and, therefore, studies on the characteristics of BEC-which is a key element of DA in the field-using the air quality model remain lacking. The BEC determines the spatial range in which the observation concentration is reflected in the model when DA is applied to the air pollution transport model. This is a very important factor that determines the initial conditions of transport and diffusion in the prediction time. Thus, it is necessary to interpret the radius of the influence of the DA quantitatively.
In this study, horizontal and vertical length scale (HLS and VLS, respectively) analyses of the BEC were applied to EnKF and EnSRF, respectively, and the two ensemble-based DA methods were performed. Their characteristics were investigated by comparing the results with those of the BEC applied to 3DVAR, i.e., the variational DA. The correlation between the HLS and VLS results and the meteorological variables was also analyzed to explain the flow-dependent BEC used in EnKF and EnSRF. Further, the performance of each DA method was evaluated using the 6 h prediction results. Finally, we suggested a The EnKF method estimates a true state x t k at time k using the "forecast" (or background) state x f k and observation y o k . The optimal state x a k (analysis state) can be obtained based on the condition that the variance of analysis satisfies the minimum by linearly combining the forecast and observation [27]: where the subscript i indicates the ith ensemble member in a total of N ensembles (i = 1, 2, . . . , N); the analysis state (x a i,k ) is obtained for each ensemble member; and N is set to 40, after considering the calculation cost and sampling error. Further, P f k represents the BEC of the EnKF, which includes the error correlations between different variables and those in the space of the model. The BEC includes only the error correlation for space because only PM 2.5 was used as the state variable in this study; the BEC can be modeled using a short-ensemble simulation. Further, the ensemble mean can be obtained easily by defining the degree of deviation from the mean as the error of the model, assuming that the ensemble mean is closer to the true value.
Among components related to the observation, H represents an observational operator, and it is a linear operator that calculates the value of the model corresponding to the observation point. This is necessary because the observation is not exactly located in the grid of the model when calculating the difference between the model and the observations. In this study, H, which linearly interpolates the model values of the four grids around the observation point, is applied. R k represents the observation error covariance matrix, and it includes instrument and representative errors that occur because observations at certain points do not represent the entire model grid. Considering several studies, the no correlations between the observation points was assumed; only diagonal components of the R k matrix were defined and constructed as ε t = ε 2 o + ε 2 r [25,28,29]. Further, ε o represents the instrument error, which is defined as ε o = 1 + 0.0075 * Π 0 ; Π 0 represents the observed value of PM 2.5 . The representative error (ε r = γε o √ ∆x/L) is a function of the horizontal resolution of the model (∆x = 27 km), where γ represents an adjustable scaling factor (here, γ = 0.5) and L denotes the impact radius. Considering the studies by Peng et al. [25] and Ha et al. [8], 3 km is used because most of the observation points are located in the city center. The perturbed observation (y o i ) performs random sampling from the Gaussian probability distribution function based on the defined R k ; it is assimilated with each ensemble member. K k represents the relative ratio of the model and the observation error; it is known as the Kalman gain matrix. A practical implementation was explained in a study conducted by Park et al. [13].
The EnKF builds an ensemble of the model and requires a perturbed observation (y o i ), which results in sampling errors attributed to the limited number of ensemble members. In contrast, the EnSRF does not require random sampling for the observation, and therefore, the sampling error of the observation can be minimized [30]. The difference is that the EnKF analysis ensemble (x a i,k ) is calculated separately using the ensemble means and deviations: Atmosphere 2022, 13, 160 4 of 18 where overbar ( the EnKF analysis ensemble ( , ) is calculated separately using the ensemble means and deviations: where overbar ̅ and prime ( ) represent the ensemble mean and deviation, respectively, and therefore, , = , − . , is updated without the observation data through , i.e., the "reduced Kalman gain." For the subsequent prediction, , = + , are used as initial fields. The remaining observation-related elements, observation error covariance ( ), and linear observation operator ( ), were applied in the same manner given that the random sampling of observation data was not performed.

Variational Method: 3DVAR
Grid-point statistical interpolation (GSI) version 3.6 [31] of the National Centers for Environmental Prediction (NCEP) was utilized to compare the BECs of the EnKF and EnSRF with the BEC of 3DVAR and examine its characteristics. The GSI system provides a 3DVAR process that calculates a best-fit state (or analysis state) by considering the weight of the error information of the model and observation. The state wherein the objective function satisfies the minimum is defined as the analysis state given by Most notations and meanings are the same as those in Equation (1). However, BEC does not change with respect to time in 3DVAR (it is time-invariant); the static BEC ( ) uses the NMC method. In this study, the difference between the 12 and 24 h predictions of PM2.5 was calculated for a total of 42 days to establish the BEC. The DA system developed in [32] was used to apply the GSI to the CMAQ model. Detailed numerical experiments, which include the EnKF and EnSRF, are described in Section 2.3. The GSI provides the generalized background error covariance matrix model (GEN_BE v2.0), a BEC calculation and analysis tool, and it calculates and provides HLS, VLS, and error variance values using the model grid [33].

Covariance Localization and Inflation
EnKF and EnSRF are ensemble-based DAs and they theoretically enable more accurate error estimation by considering the building of ensembles close to infinity. However, the number of ensembles is limited because of the computational cost; this implies that a sampling error always exists, and this error degrades the DA performance with two phenomena: spurious correlation and filter divergence. Localization was applied to the BEC to remove spurious correlations. Localization was used to gradually decrease the error covariances with respect to the distance from the observation point in the BEC using the decorrelation function proposed by Gaspari and Cohn [34], which is applied in several studies: where and represent functions based on distances in the vertical and horizontal directions, respectively. The component product is applied to the BEC ( ) using the Schur product operator (∘). In this study, the influence radius required for localization was set to 100 and 2 km in the horizontal and vertical directions, respectively; these values ) and prime ( ) represent the ensemble mean and deviation, respectively, and therefore, x i,k is updated without the observation data throughK k , i.e., the "reduced Kalman gain". For the subsequent prediction, x a i,k = x a k + x a i,k are used as initial fields. The remaining observation-related elements, observation error covariance (R), and linear observation operator (H), were applied in the same manner given that the random sampling of observation data was not performed.

Variational Method: 3DVAR
Grid-point statistical interpolation (GSI) version 3.6 [31] of the National Centers for Environmental Prediction (NCEP) was utilized to compare the BECs of the EnKF and EnSRF with the BEC of 3DVAR and examine its characteristics. The GSI system provides a 3DVAR process that calculates a best-fit state (or analysis state) by considering the weight of the error information of the model and observation. The state wherein the objective function J satisfies the minimum is defined as the analysis state given by Most notations and meanings are the same as those in Equation (1). However, BEC does not change with respect to time in 3DVAR (it is time-invariant); the static BEC (B) uses the NMC method. In this study, the difference between the 12 and 24 h predictions of PM 2.5 was calculated for a total of 42 days to establish the BEC. The DA system developed in [32] was used to apply the GSI to the CMAQ model. Detailed numerical experiments, which include the EnKF and EnSRF, are described in Section 2.3. The GSI provides the generalized background error covariance matrix model (GEN_BE v2.0), a BEC calculation and analysis tool, and it calculates and provides HLS, VLS, and error variance values using the model grid [33].

Covariance Localization and Inflation
EnKF and EnSRF are ensemble-based DAs and they theoretically enable more accurate error estimation by considering the building of ensembles close to infinity. However, the number of ensembles is limited because of the computational cost; this implies that a sampling error always exists, and this error degrades the DA performance with two phenomena: spurious correlation and filter divergence. Localization was applied to the BEC to remove spurious correlations. Localization was used to gradually decrease the error covariances with respect to the distance from the observation point in the BEC using the decorrelation function proposed by Gaspari and Cohn [34], which is applied in several studies: where ρ v and ρ h represent functions based on distances in the vertical and horizontal directions, respectively. The component product is applied to the BEC (P f k ) using the Schur product operator (•). In this study, the influence radius required for localization was set to 100 and 2 km in the horizontal and vertical directions, respectively; these values efficiently reflect the flow-dependent features while removing the spurious correlation through sensitivity experiments.
Sampling errors attributed to the limited number of ensembles underestimate the model errors. The phenomenon that appears because of an underestimation of these errors is called filter divergence, wherein the model error is low if the predictive ensemble converges together. The Kalman filter indicates that the model is more reliable and does not Atmosphere 2022, 13, 160 5 of 18 reflect the observation information in a model; therefore, the DA effect does not appear. For the CTM, which simulates a limited area, the IC effects quickly disappear. The concentration of pollutants is largely dependent on emissions and removal processes, and it is affected by the BC. Therefore, the spread (i.e., the error size of the model) of the initial ensemble member decreases rapidly over time.
A stiff system, such as a chemical reaction, is so stable that a small perturbation can progress rapidly to a quasi-steady state and disappear [35]. In this study, a Gaussian random perturbation is applied by considering a 50% uncertainty in the IC and input emissions; however, filter divergence still occurs. To this end, various measures have been proposed for solving this problem. The relaxation-to-prior (or previous) spread (RTPS) method proposed by Whitaker and Hamil [36] was applied using where σ N MC and σ f represent the background error value of the NMC method applied to 3DVAR and the spread of the ensemble forecast, respectively. The ensemble member is expressed as X f i because it may be inflated and the spread may increase. The conventional RTPS method uses the spread at the time of ensemble prediction (prior spread) as the inflation criterion at the end of the prediction; however, we were able to set the background error calculated from the NMC method of 3DVAR as the inflation criterion. The error variance was set equally in the variational DA and BEC of the ensemble-based sequential DA. Accordingly, the correction intensity of the model could be applied equally at the observation point. However, the effects around the observation point can be reflected differently from 3DVAR because the covariance can utilize each DA feature. β pri represents the inflation strength, which was set to 1.0 in this study. The inflation was applied only before the DA after the ensemble prediction.

Numerical Models and Data
WRF version 3.8.1 [37] was used for the hourly weather input data; the initial and boundary conditions of WRF included the NCEP Final (FNL) operational global analysis data [38], which have a spatial resolution of 0.25 • at 6 h intervals. The CMAQ version 5.2.1 [39] model was used to calculate the generation/destruction and advection/diffusion of PM 2.5 . The target regions for which DA studies were conducted with the WRF-CMAQ modeling system are part of Northeast Asia, which include China, Japan, and Korea ( Figure 1). The specific physical options and parameterization methods of WRF and CAMQ were the same as those described by Park et al. [13]. The grid composition and time of the simulation model are listed in Table 1. The KORUS v5.0 inventory [40] was utilized for the anthropogenic emissions of CMAQ, and MEGAN v2.1 [41,42] was used for the biogenic emissions; further, FINN v1.5 [43] was used for the fire emissions. Ground-based observations of Korea and China were used as observation data for the DA and model evaluation. The observational data for China were obtained from the China Urban Air Quality Real-Time Data Release Platform (https://106.37.208.233:20035, last accessed 3 March 2020) provided by the Chinese Ministry of Ecology and Environment. Korea's data were obtained from the National Ambient Air Quality Monitoring Information System (NAMIS) of Korea (https://www.airkorea.or.kr, last accessed 3 March 2020). The numerical simulation was conducted from 1 May to 11 June 2016; this period was the Korea-United States Air Quality (KORUS-AQ) campaign period, which was suitable for conducting numerical simulation research as it included various weather and pollution conditions, such as air stagnation, yellow dust cases, and long-distance pollutant transportation [44]. The spin-up period for the meteorological initial fields and emissions was set for five days from 00 UTC on 26 April 2016 to 30 April 2016. The analysis was performed for a total of 42 days (6 weeks), which corresponded to the KORUS-AQ period. The experiments in which DA was not performed were defined as the control run (CTR); experiments in which DA was applied every 6 h were referred to as 3DVAR, EnKF, or EnSRF according to the applied DA technique. To evaluate the short-time predictions, we focused solely on Korea; thus, the observation data and the simulation results were collected in Korea, while the DA was carried out using all the available data of the study domain in Figure 1.
conducting numerical simulation research as it included various weather and pollution conditions, such as air stagnation, yellow dust cases, and long-distance pollutant transportation [44]. The spin-up period for the meteorological initial fields and emissions was set for five days from 00 UTC on 26 April 2016 to 30 April 2016. The analysis was performed for a total of 42 days (6 weeks), which corresponded to the KORUS-AQ period. The experiments in which DA was not performed were defined as the control run (CTR); experiments in which DA was applied every 6 h were referred to as 3DVAR, EnKF, or EnSRF according to the applied DA technique. To evaluate the short-time predictions, we focused solely on Korea; thus, the observation data and the simulation results were collected in Korea, while the DA was carried out using all the available data of the study domain in Figure 1.

Length Scale Analysis
The BEC was estimated through the 6 h predicted ensemble member because the DA was performed at 6 h intervals. The of Equation (1), which was the BEC of each DA time, was used to analyze the HLS and VLS of the EnKF and EnSRF.
The matrix describes the error covariance of the model at the observation point and all grid positions of the model. This indicates the error correlation after normalization to the maximum variance value for each observation point. The length scale is obtained through function fitting using Equation (9) as a process that best describes the average value for each distance interval (i.e., minimum residual) via plotting the error correlation according to the distance for all observation points by setting the x-axis to a distance:

Length Scale Analysis
The BEC was estimated through the 6 h predicted ensemble member because the DA was performed at 6 h intervals. The P f k H T of Equation (1), which was the BEC of each DA time, was used to analyze the HLS and VLS of the EnKF and EnSRF.
The P f k H T matrix describes the error covariance of the model at the observation point and all grid positions of the model. This indicates the error correlation after normalization to the maximum variance value for each observation point. The length scale is obtained through function fitting using Equation (9) as a process that best describes the average value for each distance interval (i.e., minimum residual) via plotting the error correlation according to the distance for all observation points by setting the x-axis to a distance:

Differences in Analytical Fields
The difference in the analysis-an initial field improved through DA-was investigated for each DA method. Typically, analysis increments, i.e., analysis minus background, are used to investigate the impact of the assimilated observations. However, differentiating two analysis fields from different DA methods is also suitable to characterize not only the different impacts of observations on the improved initial conditions with respect to DA methods but also flow-dependent corrections using the ensemble-based DA. Figure 2 shows only the surface layer of the model by averaging differences in the analysis fields for different DA methods at each DA time point (168 times in total) for the entire period. The difference in the analysis field between the EnKF and EnSRF was very small, with an average of 0.01 µg m −3 (a maximum of 3.00 µg m −3 and a minimum of −2.28 µg m −3 ) in the surface layer. Therefore, only the difference between 3DVAR for EnKF and EnSRF is shown in Figure 2a,b, respectively. The difference in the concentration of the analysis field at the observation point was close to zero because the model errors at the point where the observation data were located were considered equally. The 3DVAR modified the initial field of the model to an isotropic shape based on the HLS set for the grid around the observation data; however, considering the EnKF and EnSRF, which are ensemble-based DAs, the observation data improved the model flow dependently. Thus, the area where the difference in the analysis field occurred appeared around the observation point. The average difference in the analytical field was 4-10 µg m −3 in the sea around the Korean Peninsula. This was reiterated by the absence of an observation point. The difference in the concentration of the analysis field had a positive value in most areas, which implied that EnKF or EnSRF predicted a higher value of PM 2.5 compared with that of 3DVAR based on the transportation and diffusion for the prediction time after the DA (analysis time). The difference was greater in China, which had a higher concentration of PM 2.5 .

Length Scale Analysis and Relationships with Meteorological Variables
The HLS and VLS analyses of the BEC were performed to understand the difference between the analysis fields of the two ensemble-based methods and those of the 3DVAR discussed above. Figure 3 illustrates how much the error correlation of the model decreased with an increase in the distance of the horizontal grid distance away from the observation point and the number of vertical layers along the x-axis using the BEC of the EnKF. Only the results at the time when the estimated HLS and VLS recorded the minimum and maximum values are presented because the BEC was calculated every 6 h. Moreover, the results for all the observation points at the corresponding distance were used, and therefore, they were displayed using boxes and ticks to refer to the minimum, first quartile, median, third quartile, and maximum values. The solid blue line represents a function that best describes the average value; it includes both the HLS and VLS. For example, for the HLS, was 3.5 (for VLS, was 6.4). The effect of the ground observation value obtained via DA indicated that the error correlation of the model was reduced to in the surrounding 3.5 grid (6.4 upper layers). Further, it can be observed that the HLS had a very wide range of data at the corresponding grid distance because it included the results of all regions; the error correlation could be generated in a specific direction in the two-dimensional plane, which depended on the wind direction. This feature was unlike the features of the BEC of 3DVAR. Although the HLS range for the entire period was 1.8-3.5, that of the VLS was 5.4-6.4, and it indicated only a one-layer difference. However, considering the physical distance, HLS had a 48.6-94.5 km change in scale; the heights of the fifth, sixth, and seventh floors from the ground were 630, 980, and 1490 m, respectively. Therefore, the difference between 5.4 and 6.4 in VLS is an important difference in determining whether the effect of ground observation data affected only the lower layer of the whole planetary boundary layer. The EnSRF results of the same analysis are illustrated in Figure 4; the date and time when the maximum and minimum HLSs and VLSs were shown during the numerical experiment period were consistent. Although the difference between the VLS of EnSRF and that of EnKF was very small, the HLS of EnSRF showed a difference that was as large as the 0.4 grid distance compared with that of the EnKF. Considering the average difference for the entire DA time, the HLS of EnSRF was as large as the 0.25 grid distance compared with that of the EnKF. However, the VLS of the EnSRF was as high as 0.04 layers, which indicated that the influence of the

Length Scale Analysis and Relationships with Meteorological Variables
The HLS and VLS analyses of the BEC were performed to understand the difference between the analysis fields of the two ensemble-based methods and those of the 3DVAR discussed above. Figure 3 illustrates how much the error correlation of the model decreased with an increase in the distance of the horizontal grid distance away from the observation point and the number of vertical layers along the x-axis using the BEC of the EnKF. Only the results at the time when the estimated HLS and VLS recorded the minimum and maximum values are presented because the BEC was calculated every 6 h. Moreover, the results for all the observation points at the corresponding distance were used, and therefore, they were displayed using boxes and ticks to refer to the minimum, first quartile, median, third quartile, and maximum values. The solid blue line represents a function that best describes the average value; it includes both the HLS and VLS. For example, for the HLS, L h was 3.5 (for VLS, L v was 6.4). The effect of the ground observation value obtained via DA indicated that the error correlation of the model was reduced to e −1 in the surrounding 3.5 grid (6.4 upper layers). Further, it can be observed that the HLS had a very wide range of data at the corresponding grid distance because it included the results of all regions; the error correlation could be generated in a specific direction in the two-dimensional plane, which depended on the wind direction. This feature was unlike the features of the BEC of 3DVAR. Although the HLS range for the entire period was 1.8-3.5, that of the VLS was 5.4-6.4, and it indicated only a one-layer difference. However, considering the physical distance, HLS had a 48.6-94.5 km change in scale; the heights of the fifth, sixth, and seventh floors from the ground were 630, 980, and 1490 m, respectively. Therefore, the difference between 5.4 and 6.4 in VLS is an important difference in determining whether the effect of ground observation data affected only the lower layer of the whole planetary boundary layer. The EnSRF results of the same analysis are illustrated in Figure 4; the date and time when the maximum and minimum HLSs and VLSs were shown during the numerical experiment period were consistent. Although the difference between the VLS of EnSRF and that of EnKF was very small, the HLS of EnSRF showed a difference that was as large as the 0.4 grid distance compared with that of the EnKF. Considering the average difference for the entire DA time, the HLS of EnSRF was as large as the 0.25 grid distance compared with that of the EnKF. However, the VLS of the EnSRF was as high as 0.04 layers, which indicated that the influence of the observational data appeared in the horizontal direction because of the difference between the two DA methods.   Figures 3 and 4, the date and time when the maximum/minimum of the HLS or VLS appeared were different. The HLS was analyzed in association with the wind speed (WS) and the VLS was analyzed in association with the planetary boundary layer height (PBLH) because the causes of this difference may be attributed to meteorological factors. Changes in the daily average HLS, VLS, WS, and PBLH are shown in Figure 5. The dates 23 and 25 May and 11 June were excluded because all ground observation data of China were missing more than twice out of the four DAs daily. The average of all points of the daily average WS tended to significantly decrease until 9 May, followed by an increase and a gradual decrease (Figure 4a). The tendency of the HLS to change was closely related to these fluctuations in the WS; the HLS of EnSRF was consistently larger than that of EnKF. Figure 4b shows changes in the daily average VLS and PBLH at all points. From May to June, the altitude of the boundary layer tended to gradually increase as summer approached and the solar radiation became stronger. From 7 to 10 May, precipitation appeared as the low-pressure system passed from China to Korea; therefore, the altitude of the boundary layer decreased and subsequently increased. On 27 and 28 May, the altitude of the boundary layer decreased because of the precipitation in central and southern China. Owing to these daily changes, the VLS increased or decreased, i.e., the VLS continued to increase in early June, which was less than the difference between the HLS of EnSRF and that of EnKF. However, the difference widened in early June when the VLS increased compared with that in May.  As shown in Figures 3 and 4, the date and time when the maximum/minimum of the HLS or VLS appeared were different. The HLS was analyzed in association with the wind speed (WS) and the VLS was analyzed in association with the planetary boundary layer height (PBLH) because the causes of this difference may be attributed to meteorological factors. Changes in the daily average HLS, VLS, WS, and PBLH are shown in Figure 5. The dates 23 and 25 May and 11 June were excluded because all ground observation data of China were missing more than twice out of the four DAs daily. The average of all points of the daily average WS tended to significantly decrease until May 9, followed by an increase and a gradual decrease (Figure 4a). The tendency of the HLS to change was closely related to these fluctuations in the WS; the HLS of EnSRF was consistently larger than that of EnKF. Figure 4b shows changes in the daily average VLS and PBLH at all points. From May to June, the altitude of the boundary layer tended to gradually increase as summer approached and the solar radiation became stronger. From 7 to 10 May, precipitation appeared as the low-pressure system passed from China to Korea; therefore, the altitude of the boundary layer decreased and subsequently increased. On 27 and 28 May, the altitude of the boundary layer decreased because of the precipitation in central and southern China. Owing to these daily changes, the VLS increased or decreased, i.e., the VLS continued to increase in early June, which was less than the difference between the HLS of EnSRF and that of EnKF. However, the difference widened in early June when the VLS increased compared with that in May. The average diurnal variation in the HLS and VLS is presented in Figure 6 to examine the horizontal/vertical correlation between day and night model errors. The HLS was analyzed using the daily change in WS; the VLS was analyzed using the daily change in PBLH. The box plot displays the data distribution for 42 days after taking the average of all points, showing the minimum, first quartile, median, third quartile, and maximum values from the bottom. The symbol represents an average of 42 days. The results are displayed at the bottom of the graph to compare the length scale values of the 3DVAR derived from the NMC. Further, 00 and 06 UTC correspond to daytime in China (08 and 14 in local time) and Korea (09 and 15 in local time), and 12 and 18 UTC corresponds to nighttime. Both WS and PBLH showed the highest values at midday and appeared low at night when the atmosphere was stabilized. This daily change tendency of the weather variables was consistent with that of the HLS and VLS. As indicated earlier, the length scale of the EnSRF was larger than that of the EnKF; the difference was noticeable in the HLS. Further, the length scale difference between the two DA methods increased at night compared with that during daytime. The length scale of 3DVAR was very small compared with those of the EnKF and EnSRF. The HLS of 3DVAR, which can be observed from the BEC diagnostic results provided by the GSI package, was a 1.5 grid distance, approximately half of the overall mean values of EnKF and EnSRF (2.79 and 3.05, respectively). The VLS had 3.6 layers in the 3DVAR. However, the EnKF and EnSRF averaged 5.94 and 5.99 layers, respectively, which indicates that the effect of the ground observation data on the upper layer could be limited in the 3DVAR. Further, this suggests that the difference between the analytical field of the variational method (3DVAR) and the ensemble-based method (EnKF and EnSRF) examined in Figure 2 was attributed to the difference between HLS and VLS. The average diurnal variation in the HLS and VLS is presented in Figure 6 to examine the horizontal/vertical correlation between day and night model errors. The HLS was analyzed using the daily change in WS; the VLS was analyzed using the daily change in PBLH. The box plot displays the data distribution for 42 days after taking the average of all points, showing the minimum, first quartile, median, third quartile, and maximum values from the bottom. The symbol represents an average of 42 days. The results are displayed at the bottom of the graph to compare the length scale values of the 3DVAR derived from the NMC. Further, 00 and 06 UTC correspond to daytime in China (08 and 14 in local time) and Korea (09 and 15 in local time), and 12 and 18 UTC corresponds to nighttime. Both WS and PBLH showed the highest values at midday and appeared low at night when the atmosphere was stabilized. This daily change tendency of the weather variables was consistent with that of the HLS and VLS. As indicated earlier, the length scale of the EnSRF was larger than that of the EnKF; the difference was noticeable in the HLS. Further, the length scale difference between the two DA methods increased at night compared with that during daytime. The length scale of 3DVAR was very small compared with those of the EnKF and EnSRF. The HLS of 3DVAR, which can be observed from the BEC diagnostic results provided by the GSI package, was a 1.5 grid distance, approximately half of the overall mean values of EnKF and EnSRF (2.79 and 3.05, respectively). The VLS had 3.6 layers in the 3DVAR. However, the EnKF and EnSRF averaged 5.94 and 5.99 layers, respectively, which indicates that the effect of the ground observation data on the upper layer could be limited in the 3DVAR. Further, this suggests that the difference between the analytical field of the variational method (3DVAR) and the ensemble-based method (EnKF and EnSRF) examined in Figure 2 was attributed to the difference between HLS and VLS.

As shown in
(a) HLS and WS (b) VLS and PBLH Figure 6. Diurnal variations in the horizontal length scales for 3DVAR, EnKF, and EnSRF with wind speed at each time of the day (a). Same as (a) but includes those of vertical length scales with the PBL height (b). The fixed-length scales estimated from the NMC method and used in 3DVAR are also shown.

Short-Term Predictions and Vertical Length Scale Adjustment for the 3DVAR
The results were compared with the observation results to examine how the characteristics of the analyzed length scale affected the DA performance. The short-term PM2.5 prediction performance was evaluated based on the results of prediction times (from +01 Figure 6. Diurnal variations in the horizontal length scales for 3DVAR, EnKF, and EnSRF with wind speed at each time of the day (a). Same as (a) but includes those of vertical length scales with the PBL height (b). The fixed-length scales estimated from the NMC method and used in 3DVAR are also shown.

Short-Term Predictions and Vertical Length Scale Adjustment for the 3DVAR
The results were compared with the observation results to examine how the characteristics of the analyzed length scale affected the DA performance. The short-term PM 2.5 prediction performance was evaluated based on the results of prediction times (from +01 h to +06 h) because the DA was performed four times a day. Further, the tendency of daily concentration change in PM 2.5 may be different because there was a difference in the time zone between China and Korea alongside the weather conditions based on sunrise and sunset. Thus, the prediction performance presented only the model results corresponding to the observation points in Korea. Figure 7 shows the diurnal variation in the PM 2.5 , averaged for 169 stations and 6 weeks in Korea. Along with the observed values, the results of the 3DVAR, EnKF, and EnSRF in the experiments where DA was not performed and those where DA was performed at 6 h intervals are displayed simultaneously. The PM 2.5 tended to decrease during the daytime (00 and 06 UTC) because the boundary layer altitude increased after sunrise and the wind speed became stronger. Furthermore, a weak wind speed and low boundary layer altitude appeared because the atmosphere became stable after sunset; PM 2.5 tended to increase again at night (12 and 18 UTC). These characteristics could be confirmed in the observations, and they were more clearly reflected in the model. The average concentration was underestimated when DA was not performed (CTR experiment) because not all emission sources were considered, the emission calculation methods were inaccurate, and there was uncertainty in the generation/destruction mechanism of the aerosol. Through DA, these errors were improved in all techniques. Although the predictability of the 3DVAR in the daytime was better than that of the CTR, its predictability was less than that of the EnKF and EnSRF. The improved initial fields at 00 UTC and 06 UTC through the 3DVAR DA tended to approach the CTR because the effects of the improved ICs disappeared rapidly. In contrast, the initial fields of 00 and 06 UTC for the EnKF and EnSRF exhibited very little bias at approximately 06 and 12 UTC; these were +06 h predictions, and the effect persisted. The predictability of the EnSRF was slightly better than that of the EnKF; however, the difference was insignificant. Hence, the two ensemble-based DA methods exhibited similar performances.
Atmosphere 2022, 13, x FOR PEER REVIEW 13 of 19 and the effect persisted. The predictability of the EnSRF was slightly better than that of the EnKF; however, the difference was insignificant. Hence, the two ensemble-based DA methods exhibited similar performances. We confirmed that the cause of the rapid disappearance of the influence of the analysis fields of 3DVAR was the low VLS value. Therefore, the DA was additionally performed by increasing the VLS of 3DVAR to approximately six layers, which was the average value in the ensemble DA method. Figure 7 shows the results of performing 3DVAR with 1.5 and 2.0 times the VLS calculated by the NMC method. The most notable change We confirmed that the cause of the rapid disappearance of the influence of the analysis fields of 3DVAR was the low VLS value. Therefore, the DA was additionally performed by increasing the VLS of 3DVAR to approximately six layers, which was the average value in the ensemble DA method. Figure 7 shows the results of performing 3DVAR with 1.5 and 2.0 times the VLS calculated by the NMC method. The most notable change shown after adjusting by 1.5 and 2.0 times was the improvement in the short-term prediction performance during the daytime. Because of the increased VLS, the influence of the ground observation data on the analytical sites of 00 and 06 UTC affected up to the altitude of the boundary layer, which spread to the surrounding areas through the process of transportation and diffusion at the prediction time, thereby enabling the prediction of PM 2.5 to be close to the observation.
The average PM 2.5 vertical distribution at all observation points was analyzed to investigate changes in the upper boundary layer in the analytical field ( Figure 8). Along with the results of the CTR, the +06 h prediction (dotted line) results and +00 h analysis (solid line) results for each DA method are displayed simultaneously to identify the changes in the distribution of the vertical PM 2.5 before and after the DA. Considering the results of the nonadjusted VLS of the 3DVAR (Figure 8a), the ground analysis value was close to the average observation concentration and shows the same value as that of the EnKF and EnSRF. However, the analysis profile showed a lower concentration distribution compared with those of the two ensemble-based DAs from approximately 200 to 1000 m. The analysis profile (solid green line) of 3DVAR gradually approached the ensemble DA experiment because the VLS was increased (Figure 8b,c). Accordingly, the distribution of the predicted concentration of +06 h was increased, and it approached the observed concentration on the ground, which resulted in improved predictability. The increase in the concentration from the dotted to the solid line was the average analysis increment at the time of the DA. The analysis increments of EnKF and EnSRF were lower because of the accumulated DA effect, whereas 3DVAR using the basic values of VLS showed poor short-term predictability, which indicated that significant increments were required. the accumulated DA effect, whereas 3DVAR using the basic values of VLS showed poor short-term predictability, which indicated that significant increments were required. The statistical analysis was performed at all times and points excluding 00, 06, 12, and 18 UTC when the DA was performed. This means that the observation data at the prediction time from +01 h to +06 h were used for statistical evaluations with independent observations over time. The results are summarized in Table 2. The definitions of each statistical metric are presented in the table below. When DA was applied instead of CTR, the figures improved for all metrics. The ensemble DA performed better than 3DVAR, and the EnSRF showed slightly better results than that of EnKF. An error of −30.7% was The statistical analysis was performed at all times and points excluding 00, 06, 12, and 18 UTC when the DA was performed. This means that the observation data at the prediction time from +01 h to +06 h were used for statistical evaluations with independent observations over time. The results are summarized in Table 2. The definitions of each statistical metric are presented in the table below. When DA was applied instead of CTR, the figures improved for all metrics. The ensemble DA performed better than 3DVAR, and the EnSRF showed slightly better results than that of EnKF. An error of −30.7% was observed before the DA was applied considering the mean bias normalized by the mean concentration of the observed PM 2.5 (normalized mean bias (NMB)) in Korea. However, the errors were systematically reduced to −13.0%, −4.4%, and −2.5% because 3DVAR, EnKF, and EnSRF were applied, respectively. There was also a systematic improvement in all the statistical metrics because the VLS of 3DVAR was applied 1.5 and 2.0 times. This implies that the VLS may be underestimated and needs to be when the BEC is estimated using the NMC method and applied to the 3DVAR. However, even the experiment with double VLS did not outperform the two ensemble-based DA methods.

Discussion
A length-scale analysis was performed to investigate the BEC characteristics of the EnKF and EnSRF, and the DA methods of the ensemble-based KF. The length scale diagnosed in the NMC method for calculating the BEC of the 3DVAR-a variational DA-was compared, and the correlation with the meteorological variables was examined. The difference in the analysis fields was negligible at ground observation points because the inflation of the Kalman filter was applied based on the 3DVAR's model error variance. However, it was noticeable in the surrounding areas, especially over the sea.
The results of the length scale analysis in the Gaussian function fitting method indicated that the HLS showed an average of 2.788 and 3.050 grids for the EnKF and EnSRF, respectively. This implied that the EnSRF had a wider impact on the observation data. Considering the VLS, the EnKF averaged 5.943 and EnSRF 5.989 layers, which suggests that was no significant difference between the two DA methods; the ground observation data were significantly corrected to approximately 600 m, i.e., the sixth layer. For the two ensemble-based DA methods, the HLS showed a high correlation with the ground WS, and the VLS showed a high correlation with the PBLH for diurnal and daily variations. This result implies that the DA could be performed by applying "errors of today", reflecting the weather conditions at the time of the DA when estimating the BEC through an ensemble prediction.
Our results confirmed a clear behavior of flow-dependent BEC in the CTM, whereas previous studies examined the predictability using DA techniques. The 6 h PM 2.5 prediction results performed for 42 days were evaluated simultaneously for the CTR without DA, 3DVAR, EnKF, and EnSRF. Although all DA experiments significantly reduced the negative bias compared with CTR, EnKF and EnSRF showed superior performances. Further, the VLS of 3DVAR was lower than that of the two ensemble DA methods with 3.6 layers; therefore, the sensitivity to short-term predictability evaluation was investigated further by increasing the VLS. The uncontrolled experiments showed an NMB of −13.0%, whereas the predictability increased to −10.9% and −9.8% when the VLS was applied 1.5 and 2.0 times. This was a low bias compared with the CTR of −30.7%; however, it fell short of the EnKF and EnSRF of −4.4% and −2.5%, respectively. This implied that the horizontal flow-dependent characteristic differences were also important. In terms of the VLS for 3DVAR, we suggest that the VLS should be two times larger than that estimated by the NMC method based on our sensitivity simulations.
The length scale analysis method presented in this study is a simple method for finding a Gaussian function that best describes the error correlation of the BEC based on the physical distance. However, the GSI's gen_be program estimates the length scale by directly solving the Gaussian function [33]. Although there were differences in the analysis methods in this study, it is more appropriate to apply the best-fitting method because various BECs were calculated for each scenario using the ensemble-based DA. Thus, the results of this study were not limited to predictability comparisons based on DA methods.
The results were significant in providing crucial basic information because they presented a correlation between meteorological variables and the scale analysis of the BEC, which is the core of the DA technology. However, it is necessary to use the same estimation method for direct comparison through strict length scale analysis. We plan to evaluate the predictability based on the DA method through a long-term prediction longer than the 6 h prediction (e.g., 24 h or 48 h) considered in this study. Other approaches for dimension-reduced DA other than the EnKF can be used for the CTM applications to reduce the sampling error, such as the reduced rank square root filter (RRSQRT) [45], singular evolutive Kalman filter (SEEK) [46], Karhunen-Loeve-based Kalman filter (KLKF) [47], and EnKF based on generalized polynomial chaos (gPC) [48]. In particular, to capture non-Gaussian features in a nonlinear model, dynamically orthogonal (DO) field equations [49], a mixture ensemble filter (MEnF) [50], and polynomial chaos expansion (PCE) [51] can be promising advanced DA techniques using the Gaussian mixture model (GMM) for the dimension-reduced methods. Furthermore, we plan to improve the PM 2.5 predictability by applying a hybrid DA that combines the advantages of the ensemble-based DA (i.e., time-evolving forecast errors) and those of the 3DVAR by calculating the analysis fields using a variational method [28].

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.