A Comparison between 3 DVAR and EnKF for Data Assimilation Effects on the Yellow Sea Fog Forecast

The data assimilation method to improve the sea fog forecast over the Yellow Sea is usually three-dimensional variational assimilation (3DVAR), whereas ensemble Kalman filter (EnKF) has not yet been applied to this weather phenomenon. In this paper, two sea fog cases over the Yellow sea, one spread widely and the other spread narrowly along the coastal area, are studied in detail by a series of numerical experiments with 3DVAR and EnKF based on the Grid-point Statistical Interpolation (GSI) system and the Weather Research and Forecasting (WRF) model. The results show that the assimilation effect of EnKF outperforms that of 3DVAR: for the widespread-fog case, the probability of detection and equitable threat scores of the forecasted sea fog area are improved respectively by ~57.9% and ~55.5%; the sea fog formation of the other case completely mis-forecasted by 3DVAR was produced successfully by EnKF. These improvements of EnKF relative to 3DVAR benefit from its flow-dependent background error covariances, resulting in more realistic depiction of sea surface wind for the widespread-fog case and better moisture distribution for the other case in the initial conditions. More importantly, the correlation between temperature and humidity in the background error covariances of EnKF plays a vital role in the response of moisture to the assimilation of temperature, which leads to a great improvement in the initial moisture conditions for sea fog forecast. Extra sensitivity experiments of EnKF indicate that the forecast result is sensitive to ensemble inflation and localization factors, in particular, highly sensitive to the latter.


Introduction
Sea fog usually refers to the fog that occurs over the ocean or a coastal region [1,2].It causes atmospheric horizontal visibility less than 1 km and even to tens of meters, which has a serious effect on harbor activities and marine transport.Among the seas of China, the Yellow Sea experiences a high frequency of sea fog [3,4].Numerical modeling is already becoming the major approach to both investigating the formation mechanism and developing a forecast method for the Yellow Sea fog [5,6].
Previous studies revealed that the sea fog simulation is extremely sensitive to the errors of initial conditions [3,[7][8][9][10][11][12][13][14][15].Therefore, it is very necessary to provide better initial conditions for sea fog simulation as best as we can via data assimilation.The sea fog simulation result can be improved to some extent by assimilating many types of observations using the three-dimensional variational data assimilation (3DVAR) method based on the Weather Research and Forecasting (WRF) model.These assimilated observations include routine measurements along the coast of the Yellow Sea [15], sea surface wind data aboard QuikSCAT (Quick Scatterometer) [16], satellite radiance [17], Doppler radar radial velocity [18], as well as temperature and moisture profiles derived from satellites [19].Due to a lack of routine observations over sea, a new method was proposed by Wang et al. [5] to assimilate satellite-derived humidity from the observed sea fog over the Yellow Sea, which improves the sea fog nowcasting skill with the increase of equitable threat scores (ETS) ranging from 15% to 20%.
However, due to the employment of static background error covariances, there is a weakness in the 3DVAR method that the above studies feature essentially.The NMC (National Meteorological Center) method [20] is generally used to generate these background error covariances using forecast differences.As pointed out by Bouttier [21], the NMC method is suitable only for estimation of climatological covariances.Because forecast differences are usually calculated over a reasonably long period of time (e.g., half or a month of forecast differences), the variation of the background error in different synoptic situations from one case to the next is neglected.
The limitations of static background error covariances make 3DVAR unable to manifest the flow-dependent feature.The data assimilation (DA) methods, such as Hybrid-3DVAR [22][23][24], 4DVAR [25,26] and EnKF (Ensemble Kalman Filter) [27], can provide flow-dependent background error covariances, and they have been widely employed in simulations of hurricane and typhoon that have a strong flow structure [28][29][30][31][32]. Past researches have shown that the information in the observations can be spread more reasonably (i.e., flow-dependently) by EnKF than by 3DVAR.
The Yellow Sea fog season usually starts in March and ends in August, and its dominant type belongs to advection fog [3].The sea fog generally forms when water vapor condenses within the marine atmospheric boundary layer under an appropriate synoptic system, for instance, an isolated high-pressure system or a transient cyclone/anticyclone couplet [3,33].Atmospheric flow near the sea surface dominated by the synoptic system perhaps plays an important role in reasonably spreading the observed information during the entire assimilation period, especially the information of the numerous routine observations located along the coast.Sea fog has become a severe marine weather and its influences on marine activities might even compete with tropical cyclones [34].However, up till now, application of the advanced flow-dependent DA methods to sea fog numerical modeling is seldom seen.
The National Centers for Environmental Prediction (NCEP) developed an efficient tool-the Grid-point Statistical Interpolation (GSI)/EnKF system-that includes built-in 3DVAR and EnKF [35], which provides a strong research approach for DA study in sea fog numerical modeling.The goal of the present study is to explore and compare the assimilation effects of 3DVAR and EnKF in sea fog forecast based on the GSI/EnKF system with the WRF model.
This paper is organized as follows.Section 2 briefly describes the methods of 3DVAR and EnKF.Section 3 shows how sea fog data assimilation and forecast experiments are conducted, including data, study cases, model configurations, design of DA schemes and numerical experiments.In Section 4, results and analysis are presented, and comparison of DA effects between 3DVAR and EnKF are addressed in detail.Finally, summary and conclusions are given in Section 5.

3DVAR
Let x, x b , x a , and y be the model state at the beginning of the assimilation window, a background or prior estimate of x, analysis of x and observation, respectively.The cost function for 3DVAR can be defined by where H is the observation operator matrix which transforms data from model space to observation space, R is the observation error covariance matrix, and B is the background error covariance matrix.
By using the increment formulation [36], the analysis increment and observation innovation are respectively defined as δx = x a − x b and d = y − H(x b ).Thus, Equation 1 can be rewritten The 3DVAR procedure in the GSI/EnKF system calculates x a by minimizing Equation ( 2) with the iterative solution method [35].
An accurate way to quantify R is to estimate it statistically for each specific model domain and observation type [37].However, for many observation types, R can be much more complicated [38] and is difficult to accurately quantify.Thus, in this study, R is assumed to be static and diagonal, and provided prior to the assimilation.B is calculated by statistics, and the NMC method [20] is used here to generate it as follows where x 12h t and x 24h t are respectively 12 h and 24 h model forecasts at the time of t, and T n is the time period for statistics, which is usually 15-30 days.For the WRF model, the array size of x b is usually ~1.0 × 10 7 for a typical sea fog modeling (e.g., the sizes are respectively ~1.01 × 10 7 and ~0.58 × 10 7 for the two cases in this study), thus directly solving the inversion of B (i.e., B −1 ) requires ∼ O (10 14 ) times of calculations, which is technically impossible.However, B is usually diagonalized by control variable transform [20], which uses length scale coefficients to store the correlations between different grids and regression coefficients to store the correlations between different variables.Note that these coefficients are regionally averaged to decrease the computational cost and B is therefore static and nearly homogeneous and isotropic.

EnKF
The background error covariance that is represented by the matrix P b in EnKF, is calculated using ensemble forecast members.Each ensemble member and Kalman gain matrix are respectively updated by where R and H are the same as in Equation ( 1), x b i is the i-th ensemble member (let m be their total number, so that i = 1, 2, . . ., m), x a i is the i-th updated member, and y i is the i-th perturbed observation for x b i in a way consistent with R. To avoid calculating the large matrix P b , EnKF usually approximately calculates P b H T and HP b H T as According to the least square theory, the analysis error is calculated by In Equation ( 4), the ensemble can obtain increased spread because of the perturbation of observations.However, the introduction of perturbed observations adds an additional source of sampling error related to the estimation of the observation errors.In addition to reducing the accuracy of the analysis error covariance estimate, this extra source of sampling error increases the probability that the analysis error covariance will be underestimated.Whitaker and Hamill [39] proposed the EnSRF (Ensemble Square Root Filter) method, in which perturbation of observations is cancelled and a "reduced" Kalman gain is used instead to update deviations from the ensemble mean.Equation ( 4) can be rewritten as where is used as a coefficient to guarantee that Equation ( 10) is consistent with Equation (8).The EnSRF method is adopted by the GSI/EnKF system, and was chosen to use in this paper.
Limited ensemble size generally brings a sampling error, which leads to systematically underestimated analysis error and eventually results in filter divergence [40].Another consequence of sampling error is unrealistic long-distance correlations in the background error.Thus, ensemble inflation and localization are necessary for EnKF [41].A multiplicative inflation method [42] is applied to inflate posterior ensemble spread: where σ b and σ a are the prior and posterior ensemble standard deviations, respectively; β is the inflation parameter that controls the degree of inflation.If β = 1, the ensemble is inflated completely so the posterior standard deviation becomes the same as the prior one.If β = 0, there is no inflation.The localization can be implemented by inserting a smooth correlation function ρ [43] into the Kalman gain matrix (Equation ( 5)): (12) where the symbol • represents the Schur (element-wise) product of two matrices, and ρ is a function of the distance between the model grid point and the observation location (denoted by z).Let L be the localization scale, and the method proposed by Gaspari and Cohn [44] is used to define ρ: the covariance at the distance of L will be cut to 20.8%, and it will be cut to 0 at the distance of 2L.
EnKF sequentially assimilates each observation, and treats the analysis updated by one observation as the background for the next assimilation.As a result, y, R, and HP b H T become scalers, and P b H T is reduced to a vector.This brings the benefit that not only is matrix inversion avoided but also there is no need to simplify background error covariances (e.g., diagonalizing), which results in intact correlations between physical variables.Since P b H T is updated along with assimilating observations, it means the background error covariances may vary substantially depending on the flow of the day.

Data
The surface synoptic charts were from KMA (Korea Meteorological Administration), demonstrating the weather system and atmospheric flow near the sea surface for the studied sea fog cases.The observed fog patches were retrieved empirically from the infrared, albedo, and visible cloud imageries of MTSAT (the Multifunctional Transport Satellite; [45]) using the method by Wang et al. [5].
The initial and lateral boundary conditions for the WRF simulation later were derived from the NCEP Final Analysis (FNL; 1 • × 1 • , 6-hourly; [46]), and the sea surface temperature (SST) data were extracted from daily NEAR-GOOS (North-East Asian Regional Global Ocean Observing System) dataset [47].Observations for data assimilation include radiosonde and surface measurements, satellite-retrieved sea surface winds, and radiation brightness temperatures from AMSU-A/B (Advanced Microwave Sounding Unit A/B), HIRS-3/4 (High Resolution Infrared Radiation Sounder 3/4) and MHS (Microwave Humidity Sounder) carried by satellites [48].
Besides the radiosonde and surface measurements mentioned above, CCMP (the Cross Calibrated Multi-Platform) global surface wind data [49] and NCEP ship observations in prepBUFR format [50] were used for evaluation.

Sea Fog Cases
Two typical advection sea fog cases over the Yellow Sea, one spread widely over the Yellow Sea and the other spread narrowly along the Shandong Peninsula coast (Figure 1), were chosen for study in detail.They happened on 10 April 2009 and 2 April 2014 (hereafter C-09 and C-14), respectively.Although both the two cases were dominated by a high-pressure system (Figure 1a,c), they have obvious differences.In the case C-09, the widespread fog area was formed as warm moist air moved onto a cold sea area (see the SST in Figure 2a) due to easterly winds in the south of the high (Figure 1a).But in the case C-14, weak southerly winds near the center of the high (Figure 1c) carried warm air from the western Yellow Sea and its neighboring land area northward onto the cool coastal area (see the SST in Figure 2b) along Shandong Peninsula, resulting in sea fog formation.The shade of cloud color shows that the fog in the case C-09 is thicker than that in the case C-14 (cf. Figure 1b,d), implying the case C-14 is not as humid as the case C-09.Perhaps it is due to their obvious difference between the prevailing flows mentioned above (cf.Figure 1a,c).

Model Configuration
The Advanced Research core of the WRF (ARW, version 3.8.1)[51] was employed in this study.According to the work of Lu et al. [52] and Wang et al. [5], a combination of the planetary boundary layer (PBL) and microphysics schemes YSU-LIN was chosen.The details of model configurations, including domains, horizontal and vertical resolutions, and other physical schemes, are listed in Table 1.Both sea fog cases used two domains (D1 and D2) and applied two-way nesting (Figure 2).

Model Configuration
The Advanced Research core of the WRF (ARW, version 3.8.1)[51] was employed in this study.According to the work of Lu et al. [52] and Wang et al. [5], a combination of the planetary boundary layer (PBL) and microphysics schemes YSU-LIN was chosen.The details of model configurations, including domains, horizontal and vertical resolutions, and other physical schemes, are listed in Table 1.Both sea fog cases used two domains (D1 and D2) and applied two-way nesting (Figure 2).Colors show sea surface temperature (SST) within the D2 domains, and the red and black dots respectively locate radiosonde stations and ship measurement sites.The sea area surrounded by thick dashed lines and coastlines in D2 of (a) is used for evaluating forecasted sea surface winds later.

Design of DA Schemes and Experiments
Three DA schemes were designed for comparing the data assimilation effects between 3DVAR and EnKF. Figure 3 shows a schematic illustration of these schemes.DA-1 and DA-3 schemes are the cycling WRF-3DVAR scheme and the cycling WRF-EnKF, respectively; DA-2 is the cycling WRF-3DVAR with multi-backgrounds from the WRF ensemble forecast, and its backgrounds at the beginning of the DA window are the same as DA-3.Among the schemes, B and P b represent the background error covariances for 3DVAR and EnKF, respectively, and obs stands for observations; Each assimilation cycle is connected by the WRF integration (i.e., wrf.exe in Figure 3), and the final assimilation analysis  a or  a ̅̅̅ are the initial conditions for the next WRF forecast.

Design of DA Schemes and Experiments
Three DA schemes were designed for comparing the data assimilation effects between 3DVAR and EnKF. Figure 3 shows a schematic illustration of these schemes.DA-1 and DA-3 schemes are the cycling WRF-3DVAR scheme and the cycling WRF-EnKF, respectively; DA-2 is the cycling WRF-3DVAR with multi-backgrounds from the WRF ensemble forecast, and its backgrounds at the beginning of the DA window are the same as DA-3.Among the schemes, B and P b represent the background error covariances for 3DVAR and EnKF, respectively, and obs stands for observations; Each assimilation cycle is connected by the WRF integration (i.e., wrf.exe in Figure 3), and the final assimilation analysis x a or x a are the initial conditions for the next WRF forecast.The initial (t = t0 in Figure 3) ensembles for DA-2 and DA-3 were created by a random perturbation method [23,24,58].In this study, the built-in method known as RANDOMCV in the WRF-DA system, was employed to generate ensemble initial conditions (ICs) by adding random noise to the analysis in the control variable space.The ensemble ICs were generated 6 hours prior to the starting time of the assimilation window, and then the WRF runs with these ensemble ICs were integrated for 6 hours to produce the initial ensemble.The ensemble size was set at 40 (i.e., m = 40 in Figure 3); the GEN_BE tool (Ver 2.0) [59] was employed to generate the static B with the control variables listed in Table 2.
The major difference between 3DVAR and EnKF is the background error covariances, which cause different assimilation results.However, another factor will also contribute to the assimilation results.It is the ensemble forecast of DA-3 in the DA window, which might lead to a difference of initial conditions even if the update steps are the same.Thus, the purpose of the DA-2 scheme is to isolate the contribution of the ensemble forecast relative to DA-1.
Minimization of the cost function J() in Equation 2for 3DVAR was completed by iterative solution, including the outer loop (OL) and inner loop (IL) iterations with a convergence criterion (e.g., 1.0 × 10 −5 ).For EnKF, the default value in the GSI/EnKF system for the inflation factor (β = 0.9) was taken, and the localization scale was set to 200 km instead of the default value 500 km, according to the results of the hybrid-3DVAR experiments for the Yellow Sea fog by Wang [60].The detailed configurations and differences between 3DVAR and EnKF are listed in Table 2, and the observation error variances (R) for each observation type are given in Table 3.Note that the R provided by the GSI/EnKF system in Table 2 may be sub-optimal, because R is often given by the instrument error The initial (t = t 0 in Figure 3) ensembles for DA-2 and DA-3 were created by a random perturbation method [23,24,58].In this study, the built-in method known as RANDOMCV in the WRF-DA system, was employed to generate ensemble initial conditions (ICs) by adding random noise to the analysis in the control variable space.The ensemble ICs were generated 6 hours prior to the starting time of the assimilation window, and then the WRF runs with these ensemble ICs were integrated for 6 hours to produce the initial ensemble.The ensemble size was set at 40 (i.e., m = 40 in Figure 3); the GEN_BE tool (Ver 2.0) [59] was employed to generate the static B with the control variables listed in Table 2.
The major difference between 3DVAR and EnKF is the background error covariances, which cause different assimilation results.However, another factor will also contribute to the assimilation results.It is the ensemble forecast of DA-3 in the DA window, which might lead to a difference of initial conditions even if the update steps are the same.Thus, the purpose of the DA-2 scheme is to isolate the contribution of the ensemble forecast relative to DA-1.
Minimization of the cost function J(x) in Equation 2for 3DVAR was completed by iterative solution, including the outer loop (OL) and inner loop (IL) iterations with a convergence criterion (e.g., 1.0 × 10 −5 ).For EnKF, the default value in the GSI/EnKF system for the inflation factor (β = 0.9) was taken, and the localization scale was set to 200 km instead of the default value 500 km, according to the results of the hybrid-3DVAR experiments for the Yellow Sea fog by Wang [60].The detailed configurations and differences between 3DVAR and EnKF are listed in Table 2, and the observation error variances (R) for each observation type are given in Table 3.Note that the R provided by the GSI/EnKF system in Table 2 may be sub-optimal, because R is often given by the instrument error variance (and with an ad hoc inflation factor for EnKF) to account for representivity errors for a specific model domain and observation type [38].Six sea fog forecast experiments were conducted, in which Exp-A1-Exp-A3 were designed for the case C-09 and the others were for the case C-14.Table 4 summarizes the DA schemes and assimilated observations for each experiment.The starting times of the prediction were 1200 UTC 9 April 2009 and 00UTC 2 April 2014 for the cases C-09 and C-14, respectively.Cycling 3DVAR/EnKF assimilations were performed in the DA window before the starting points.For the case C-09, the DA window was 12 h with a 6-h cycle interval and the forecast length was 24 h.The case C-14 had a 12-h forecast period with 12-h DA window and 3-h cycle interval.The WRF was initialized every 12 h and run 24 h for a period of 15 days centered at the starting day of the case forecast.The WRF results were then used to generate a new static background error covariance matrix (i.e., B in Figure 3) by the NMC method [20] with the GEN-BE tool [53].The background error uses regression coefficients to correlate wind speed, temperature, and pressure, while relative humidity is independent.Note that the flow-dependent background error covariance matrix for EnKF (i.e., P b in Figure 3), in which any two analysis variables are codependent, was statistically calculated based on the ensemble of x b and updated in every cycle.

Methods for Sea Fog Diagnostics and Evaluation
Fog area is the aspect of most concern of sea fog forecasting.However, it is impossible to determine sea fog area by limited observations from ships and buoys over the Yellow Sea.The method used by Wang et al. [5] was employed here to retrieve the sea fog area (hereafter termed "observed sea fog") using the MTSAT data.
Sea fog area diagnosed from the model results (hereafter termed "forecasted sea fog") was diagnosed by forecasted liquid water content (LWC).The region with LWC at the lowest model level ≥0.016 g/kg or the fog-top height ≤400 m is defined as forecasted sea fog [5,61,62].The value 0.016 g/kg for LWC is equivalent to 1 km horizontal visibility calculated using the formula by Stoelinga and Thomas [63].The fog-top height is the height of the most upper vertical level where the LWC is more than 0.016 g/kg in the marine PBL.
Forecasted sea fog was compared with observed sea fog by the statistical scores for a binary event (yes or no; 1 or 0), including the probability of detection (POD), false alarm ratio (FAR), bias, and equitable threat score (ETS) [5,61,62].These scores are defined by where F, O, and H, respectively, refer to the numbers of forecast points, observed points, and correctly forecast points with fog occurrence (hits; binary value is 1); R = F(O/N) is a random hit penalty, and N is the total number of points in the verification area.POD means the percentage of the observed fog which is correctly forecasted, while FAR refers to the percentage of the forecasted fog not observed.Since the bias is defined here as the ratio of total forecast points to total observed points the best value for the bias is 1.0, while values less than 1.0 indicate under prediction and values greater than 1.0 indicate over prediction.ETS is a positively oriented score; that is, the larger the ETS is, the better a forecast will be.The domain D2 (see Figure 2) was taken as the verification area.Both the forecasted and observed sea fog areas were meshed onto the same grids in the verification area, where point-to-point comparisons were conducted and the above scores calculated statistically.Because the horizontal resolutions are respectively 10 km and 6 km for the cases C-09 and C-14 (see Table 1), the grid sizes were accordingly set to 0.10 • and 0.05 • , respectively.The verification was made for all forecast outputs at 1-h intervals.Both the land and the sea areas covered by high-altitude clouds were excluded from the verification domain; the early hours of the morning and evening were skipped, because the method used by Wang et al. [5] has not the capability to retrieve sea fog patch at these hours due to the low solar zenith angle.

Verification of Sea Fog Area
Forecasted sea fog of the cases C-09 and C-14 are shown in Figures 4 and 5, respectively.Panels in the top row in both Figures 4 and 5 illustrate the observed sea fog, and panels in the other rows give the forecasted sea fog.By visually comparing the forecasted sea fog with the observed sea fog, it is clearly seen that the forecast experiments (i.e., Exp-A3 and Exp-B3) with DA-3 scheme perform best (e.g., cf. Figure 4t,o,j with Figure 4e, respectively; cf. Figure 5r,h, m with Figure 5c, respectively).Particularly, Exp-B1 that uses the DA-1 scheme completely fails to capture the forecasted sea fog (Figure 5f-j), whereas Exp-B3 successfully forecasts the sea fog formation and development in the next several hours.However, the dispersion after 7 h is not captured.Forecasted sea fog of the cases C-09 and C-14 are shown in Figures 4 and 5, respectively.Panels in the top row in both Figures 4 and 5 illustrate the observed sea fog, and panels in the other rows give the forecasted sea fog.By visually comparing the forecasted sea fog with the observed sea fog, it is clearly seen that the forecast experiments (i.e., Exp-A3 and Exp-B3) with DA-3 scheme perform best (e.g., cf. Figure 4t,o,j with Figure 4e, respectively; cf. Figure 5r,h, m with Figure 5c, respectively).Particularly, Exp-B1 that uses the DA-1 scheme completely fails to capture the forecasted sea fog (Figure 5f-j), whereas Exp-B3 successfully forecasts the sea fog formation and development in the next several hours.However, the dispersion after 7 h is not captured.Instead of using eyeball comparison, a quantitative evaluation was conducted using statistical scores defined above.Table 5 outlines the temporally-averaged statistical scores.Compared to Exp-A1, Exp-A3 show obvious improvements, in which POD, FAR, bias, and ETS are all positive and the ETS is improved by 55.5%.Similarly, compared to Exp-B1, all scores of Exp-B3 except for FAR have significant improvements especially with an ETS gain of 8320.0%.However, the experiments Exp-A2 and Exp-B2 have the opposite improvement performance: the former negative, the latter positive.Instead of using eyeball comparison, a quantitative evaluation was conducted using statistical scores defined above.Table 5 outlines the temporally-averaged statistical scores.Compared to Exp-A1, Exp-A3 show obvious improvements, in which POD, FAR, bias, and ETS are all positive and the ETS is improved by 55.5%.Similarly, compared to Exp-B1, all scores of Exp-B3 except for FAR have significant improvements especially with an ETS gain of 8320.0%.However, the experiments Exp-A2 and Exp-B2 have the opposite improvement performance: the former negative, the latter positive.As shown in Figure 3 and Table 4, Exp-A2/B2 and Exp-A1/B1 have different backgrounds at the beginning of the DA windows with the same DA scheme, while Exp-A2/B2 and Exp-A3/B3 have the same background but with different DA schemes.According to the above results of sea fog area verification, it is demonstrated that Exp-A3/B3 performs much better than Exp-A1/B1 and Exp-A2/B2 mainly due to the DA-3 scheme.

Verification by Measurements
Sea fog over the Yellow sea usually forms and develops under a high-pressure weather system with stable marine PBL and suitable moisture conditions [15].The radiosonde, surface, and ship measurements were taken to validate the analysis fields (i.e., ICs).The agreement between the initial analysis and the radiosonde observations (see their sites in Figure 2) were first evaluated, focusing on bias and root-mean-square error (RMSE).Figure 6 shows the average vertical profiles of height, temperature, and moisture (represented by the water vapor mixing ratio) over all radiosonde observations.As shown in Figure 3 and Table 4, Exp-A2/B2 and Exp-A1/B1 have different backgrounds at the beginning of the DA windows with the same DA scheme, while Exp-A2/B2 and Exp-A3/B3 have the same background but with different DA schemes.According to the above results of sea fog area verification, it is demonstrated that Exp-A3/B3 performs much better than Exp-A1/B1 and Exp-A2/B2 mainly due to the DA-3 scheme.

Verification by Measurements
Sea fog over the Yellow sea usually forms and develops under a high-pressure weather system with stable marine PBL and suitable moisture conditions [15].The radiosonde, surface, and ship measurements were taken to validate the analysis fields (i.e., ICs).The agreement between the initial analysis and the radiosonde observations (see their sites in Figure 2) were first evaluated, focusing on bias and root-mean-square error (RMSE).Figure 6 shows the average vertical profiles of height, temperature, and moisture (represented by the water vapor mixing ratio) over all radiosonde observations.7 presents the result for the case C-09, and the result of the case C-14 is quite similar but not shown here.It is clearly seen that Exp-A1 and Exp-A2 have almost the same RMSE and bias, however Exp-A3 greatly decreases its RMSE and bias (Figure 7b,c).Compared with the isobars on the KMA surface chart (Figure 7a), it is also found that the SLP distribution of Exp-A3 is much closer to the observed fact than Exp-A and Exp-B (Figure 7a,d-f).For example, the 1022.5 hPa contour in Exp-A3 is more consistent with that on the chart than either of Exp-A1 and Exp-A2.At the same time, the improved SLP promotes better sea surface wind (SSW).7 presents the result for the case C-09, and the result of the case C-14 is quite similar but not shown here.It is clearly seen that Exp-A1 and Exp-A2 have almost the same RMSE and bias, however Exp-A3 greatly decreases its RMSE and bias (Figure 7b,c).Compared with the isobars on the KMA surface chart (Figure 7a), it is also found that the SLP distribution of Exp-A3 is much closer to the observed fact than Exp-A and Exp-B (Figure 7a,d-f).For example, the 1022.5 hPa contour in Exp-A3 is more consistent with that on the chart than either of Exp-A1 and Exp-A2.At the same time, the improved SLP promotes better sea surface wind (SSW).The forecasted SSW in the region, which is indicated by dashed lines in Figure 2a, is verified by the CCMP wind product.The RMSEs of wind speed (direction) for Exp-A1, Exp-A2, and Exp-A3 are 1.30 m/s (25.74 • ), 1.17 m/s (25.03 • ), and 1.14 m/s (23.70 • ), respectively.The result shows that Exp-A3 yields better sea surface wind assimilation.
Atmosphere 2018, 9, x 14 of 29 The forecasted SSW in the region, which is indicated by dashed lines in Figure 2a, is verified by the CCMP wind product.The RMSEs of wind speed (direction) for Exp-A1, Exp-A2, and Exp-A3 are 1.30 m/s (25.74°), 1.17 m/s (25.03°), and 1.14 m/s (23.70°), respectively.The result shows that Exp-A3 yields better sea surface wind assimilation.As well as the profiles of geopotential height, Figure 6 contains the profiles of temperature and mixing ratio.For the case C-09, the profiles of mixing ratio bias and RMSE for Exp-A3 have best-performing values below 950 hPa (about 600 m; Figure 6c), as well as for temperature (Figure 6b).The average mixing ratios bias and RMSE below 950 hPa of Exp-A1/A3 are 0.19/0.13g/kg and 1.28/1.09g/kg, respectively.Although the profiles of temperature and mixing ratio for Exp-B3 become worse than Exp-B1 and Exp-B2 above 1000 hPa (near 1300 m; Figure 6e,f), there are significant improvements below 1000 hPa.Especially at the 1010 hPa level, the RMSE of the temperature drops from 1.4 K to 0.3 K, and the RSME of the mixing ratio decreases from 0.90 g/kg to 0.55 g/kg.
It is worth mentioning that the layers within the marine PBL where temperature and moisture are improved in Exp-A3 and Exp-B3 have obviously different depths: the former about 600 m, while the latter near 130 m.It should be noted that the fog in the case of C-09 is thicker than that of the case C-14, which was addressed previously at the end of Section 3.2.This consistency between the depths and the fog thicknesses indicates that the marine PBL structures in both Exp-A3 and Exp-B3 are reasonably improved.
Furthermore, the forecasted moisture over sea is verified using ship measurements (see the sites in Figure 2) during the forecast time (Figure 8).For the case C-09 (Figure 8a), the mixing ratio forecasted by Exp-A3 performs better than Exp-A1 and Exp-A2 at S1-4, and performs worse than Exp-A1 and Exp-A2 at S7. Exp-A1 and Exp-A3 have almost the same mixing ratio biases at S5 and S8.The RMSEs of Exp-A1-A3 are 0.59, 0.56, and 0.51 g/kg, respectively.For the case C-14 (Figure 8b), the biases for Exp-B3 are generally comparable or smaller than those for Exp-B1 and Exp-B2.The RMSEs of Exp-B1-B3 are 2.03, 2.04, and 1.74 g/kg, respectively.The verification result indicates that DA-3 produces a better moisture structure over sea.As well as the profiles of geopotential height, Figure 6 contains the profiles of temperature and mixing ratio.For the case C-09, the profiles of mixing ratio bias and RMSE for Exp-A3 have best-performing values below 950 hPa (about 600 m; Figure 6c), as well as for temperature (Figure 6b).The average mixing ratios bias and RMSE below 950 hPa of Exp-A1/A3 are 0.19/0.13g/kg and 1.28/1.09g/kg, respectively.Although the profiles of temperature and mixing ratio for Exp-B3 become worse than Exp-B1 and Exp-B2 above 1000 hPa (near 1300 m; Figure 6e,f), there are significant improvements below 1000 hPa.Especially at the 1010 hPa level, the RMSE of the temperature drops from 1.4 K to 0.3 K, and the RSME of the mixing ratio decreases from 0.90 g/kg to 0.55 g/kg.
It is worth mentioning that the layers within the marine PBL where temperature and moisture are improved in Exp-A3 and Exp-B3 have obviously different depths: the former about 600 m, while the latter near 130 m.It should be noted that the fog in the case of C-09 is thicker than that of the case C-14, which was addressed previously at the end of Section 3.2.This consistency between the depths and the fog thicknesses indicates that the marine PBL structures in both Exp-A3 and Exp-B3 are reasonably improved.Furthermore, the forecasted moisture over sea is verified using ship measurements (see the sites in Figure 2) during the forecast time (Figure 8).For the case C-09 (Figure 8a), the mixing ratio forecasted by Exp-A3 performs better than Exp-A1 and Exp-A2 at S1-4, and performs worse than Exp-A1 and Exp-A2 at S7. Exp-A1 and Exp-A3 have almost the same mixing ratio biases at S5 and S8.The RMSEs of Exp-A1-A3 are 0.59, 0.56, and 0.51 g/kg, respectively.For the case C-14 (Figure 8b), the biases for Exp-B3 are generally comparable or smaller than those for Exp-B1 and Exp-B2.The RMSEs of Exp-B1-B3 are 2.03, 2.04, and 1.74 g/kg, respectively.The verification result indicates that DA-3 produces a better moisture structure over sea.

Feature of Initial Condition Differences
From the results of the evaluation in Section 4.2, it can be shown that the assimilation effect using the DA-3 scheme is better than those using the other DA schemes.The two sea fog cases, one spread widely over the Yellow Sea and the other spread narrowly along the Shandong Peninsula coast, can both obtain obvious forecast improvement by using EnKF other than 3DVAR.Next, we discuss the issue of how 3DVAR and EnKF work for the improvement and how they differ in the assimilation process for the two sea fog cases.
Differences between the initial conditions of different experiments are investigated.Figure 9 shows the initial condition difference near the sea surface of Exp-A2 minus Exp-A1, and Exp-A3

Feature of Initial Condition Differences
From the results of the evaluation in Section 4.2, it can be shown that the assimilation effect using the DA-3 scheme is better than those using the other DA schemes.The two sea fog cases, one spread widely over the Yellow Sea and the other spread narrowly along the Shandong Peninsula coast, can both obtain obvious forecast improvement by using EnKF other than 3DVAR.Next, we discuss the issue of how 3DVAR and EnKF work for the improvement and how they differ in the assimilation process for the two sea fog cases.
Differences between the initial conditions of different experiments are investigated.Figure 9 shows the initial condition difference near the sea surface of Exp-A2 minus Exp-A1, and Exp-A3 minus Exp-A1 for the case C-09.It is similar in Figure 10 but here for the case C-14.It is clear there are rather small differences when the DA-2 scheme is used instead of the DA-1 scheme (see the upper panels in both Figures 9 and 10).However, when the DA scheme is changed from DA-1 to DA-3, the differences are significant (see the bottom panels in both Figures 9 and 10).In the areas that overlap or neighbor where sea fogs occur (sea fog areas are denoted by the closed thick solid lines in Figures 9 and 10), the differences become obviously larger.For instance, in Figure 10d-f, the differences of temperature, mixing ratio, wind speed, and SLP are about 1.0 K, 0.6 g/kg, 3.0 m/s, and 4 hPa, respectively.However, it is different in Figure 9, the differences in the area over the sea fog area are not significant.Instead, the area marked with a dashed frame, which locates northeast of the fog area, has large differences of mixing ratio and wind (Figure 9e,f).
After having an insight into the evolution of forecasted sea fog in the sea fog cases (Figures 4 and 5), it is found that the forecast improvement is embodied in the development stage for the case C-09 (e.g., Figure 4d,i,n,s), but for the case C-14 the forecast improvement is already spotted in the initial stage (e.g., Figure 5a,f,k,p).These different forecast improvements might be associated with the areas that have the significant differences mentioned above.Next, we will explore this aspect.For convenience of reference, the above two areas stamped by dashed frames are briefly called significant zones.
Combining the above information, two questions for further analysis are proposed here: Are the forecast improvements closely related to the significant zones?If the answer is yes to the above question, then how are the significant zones formed in the assimilation process?

Reason for the Forecast Improvements
For the case C-14, Exp-B1 completely fails, resulting in no sea fog in the whole forecast period (Figure 5).However, Exp-B3 successfully captures the sea fog process from the beginning.The positive mixing ratio difference at 1000 hPa in the significant zones (Figure 10b,e), which shows it is about 0.6 g/kg wetter in Exp-B3 that Exp-B1, indicates that the failure of Exp-B1 may be due to lack of enough moisture.On the contrary, the success of Exp-B3 results from enough moisture (Figure 10e) In the areas that overlap or neighbor where sea fogs occur (sea fog areas are denoted by the closed thick solid lines in Figures 9 and 10), the differences become obviously larger.For instance, in Figure 10d-f, the differences of temperature, mixing ratio, wind speed, and SLP are about 1.0 K, 0.6 g/kg, 3.0 m/s, and 4 hPa, respectively.However, it is different in Figure 9, the differences in the area over the sea fog area are not significant.Instead, the area marked with a dashed frame, which locates northeast of the fog area, has large differences of mixing ratio and wind (Figure 9e,f).
After having an insight into the evolution of forecasted sea fog in the sea fog cases (Figures 4  and 5), it is found that the forecast improvement is embodied in the development stage for the case C-09 (e.g., Figure 4d,i,n,s), but for the case C-14 the forecast improvement is already spotted in the initial stage (e.g., Figure 5a,f,k,p).These different forecast improvements might be associated with the areas that have the significant differences mentioned above.Next, we will explore this aspect.For convenience of reference, the above two areas stamped by dashed frames are briefly called significant zones.
Combining the above information, two questions for further analysis are proposed here: Are the forecast improvements closely related to the significant zones?If the answer is yes to the above question, then how are the significant zones formed in the assimilation process?

Reason for the Forecast Improvements
For the case C-14, Exp-B1 completely fails, resulting in no sea fog in the whole forecast period (Figure 5).However, Exp-B3 successfully captures the sea fog process from the beginning.The positive mixing ratio difference at 1000 hPa in the significant zones (Figure 10b,e), which shows it is about 0.6 g/kg wetter in Exp-B3 that Exp-B1, indicates that the failure of Exp-B1 may be due to lack of enough moisture.On the contrary, the success of Exp-B3 results from enough moisture (Figure 10e) additionally with a strengthening high-pressure (Figures 10f and 1c), which means the atmospheric stratification is more stable and suitable for sea fog occurrence.
Figure 11 further gives the observation innovation (i.e., y − H(x b )) and analysis increment (i.e., x a − x b ) at the lowest model level for Exp-B1 and at the last DA cycle (t = t 2 in Figure 3).Prior to the assimilation of this cycle, the observation innovations and analysis increments of temperature and mixing ratio differ very little between Exp-B1 and Exp-B3 (cf.Figures 11 and 11d; cf. Figure 11b,e).It means that Exp-B1 and Exp-B3 have almost the same backgrounds (i.e., x b ) of temperature and mixing ratio.However, after assimilating new observations, Exp-B1 and Exp-B3 obtain distinctly different analysis increments (cf. Figure 11c,f).In the significant zone, both Exp-B1 and Exp-B3 get a negative temperature increment, but the amplitude of the former is much smaller than that of the latter; Exp-B1 has a gain of mixing ratio about 0.1 g/kg over the land part of the significant zone, whereas the gain of Exp-B3 increases up to at least 0.2 g/kg almost over the whole significant zone.
Atmosphere 2018, 9, x 18 of 29 additionally with a strengthening high-pressure (Figures 10f and 1c), which means the atmospheric stratification is more stable and suitable for sea fog occurrence.Figure 11 further gives the observation innovation (i.e. ,  − ( b )) and analysis increment (i.e.,  a −  b ) at the lowest model level for Exp-B1 and at the last DA cycle (t = t2 in Figure 3).Prior to the assimilation of this cycle, the observation innovations and analysis increments of temperature and mixing ratio differ very little between Exp-B1 and Exp-B3 (cf. Figure 11 and Figure 11d; cf. Figure 11b,e).It means that Exp-B1 and Exp-B3 have almost the same backgrounds (i.e.,  b ) of temperature and mixing ratio.However, after assimilating new observations, Exp-B1 and Exp-B3 obtain distinctly different analysis increments (cf. Figure 11c,f).In the significant zone, both Exp-B1 and Exp-B3 get a negative temperature increment, but the amplitude of the former is much smaller than that of the latter; Exp-B1 has a gain of mixing ratio about 0.1 g/kg over the land part of the significant zone, whereas the gain of Exp-B3 increases up to at least 0.2 g/kg almost over the whole significant zone.In summary, the marine PBL below 1000 hPa of Exp-B3 is much wetter than that of Exp-B1 in the initial condition over the significant zone (Figure 9), contributing to the successful forecast of Exp-B3.However, Exp-A3 perhaps owes its forecast improvement in the developing stage (Figure 4) to the moisture advection from its significant zone, which has both higher mixing ratio and larger easterly wind (Figure 9e,f).
To further understand the relationship between the significant zone and the forecast improvement in Exp-A3, evolution of the averages of temperature, mixing ratio, and u-and v-components of wind over the significant zone during the DA window and forecast period are In summary, the marine PBL below 1000 hPa of Exp-B3 is much wetter than that of Exp-B1 in the initial condition over the significant zone (Figure 9), contributing to the successful forecast of Exp-B3.However, Exp-A3 perhaps owes its forecast improvement in the developing stage (Figure 4) to the moisture advection from its significant zone, which has both higher mixing ratio and larger easterly wind (Figure 9e,f).
To further understand the relationship between the significant zone and the forecast improvement in Exp-A3, evolution of the averages of temperature, mixing ratio, and u-and v-components of wind over the significant zone during the DA window and forecast period are demonstrated by the time series in Figure 12.At the forecast beginning (t = 0), compared to Exp-A1, Exp-A3 has lower temperature (Figure 12a), higher mixing ratio (Figure 12b), and stronger southeasterly wind (smaller uand v-components; Figure 12c,d), which means a moisture advection coming from the significant zone to the fog area in Exp-A3.This situation is maintained during the whole forecast period, resulting in continual moisture advection from the significant zone to the fog area (Figure 9e,f) and promoting the forecast improvement.
Atmosphere 2018, 9, x 19 of 29 demonstrated by the time series in Figure 12.At the forecast beginning (t = 0), compared to Exp-A1, Exp-A3 has lower temperature (Figure 12a), higher mixing ratio (Figure 12b), and stronger southeasterly wind (smaller u-and v-components; Figure 12c,d), which means a moisture advection coming from the significant zone to the fog area in Exp-A3.This situation is maintained during the whole forecast period, resulting in continual moisture advection from the significant zone to the fog area (Figure 9e,f) and promoting the forecast improvement.

Impact of the Background Error Covariances
As seen in the above reasonable analysis, the sea fog forecast improvements are clearly related to the significant zones in the two fog cases.Next, we try to explain the formation of the significant zones, and compare the impact of background error on the formation by 3DVAR and EnKF.
In the assimilation process, analysis increment represents the gain of assimilating observations jointly controlled by the background error and observation error covariances.Since observation error is usually confirmed in advance, background error covariances become the key factor to determine the analysis increment.If we focus on inspecting the analysis increments in Figures 11  and 12, it is found that there are two noteworthy facts.One fact is that increments of the 3DVAR are generally larger than the EnKF in the case C-09 (see the bars in Figure12) except for humidity.The other one is that in the case C-14 there is a see-sawing relationship of the increments between temperature and humidity in the EnKF but not in the 3DVAR.Namely that negative temperature increment corresponds to positive humidity increment, and vice versa (Figure 11f).According to the

Impact of the Background Error Covariances
As seen in the above reasonable analysis, the sea fog forecast improvements are clearly related to the significant zones in the two fog cases.Next, we try to explain the formation of the significant zones, and compare the impact of background error on the formation by 3DVAR and EnKF.
In the assimilation process, analysis increment represents the gain of assimilating observations jointly controlled by the background error and observation error covariances.Since observation error is usually confirmed in advance, background error covariances become the key factor to determine the analysis increment.If we focus on inspecting the analysis increments in Figures 11 and 12, it is found that there are two noteworthy facts.One fact is that increments of the 3DVAR are generally larger than the EnKF in the case C-09 (see the bars in Figure 12) except for humidity.The other one is that in the case C-14 there is a see-sawing relationship of the increments between temperature and humidity in the EnKF but not in the 3DVAR.Namely that negative temperature increment corresponds to positive humidity increment, and vice versa (Figure 11f).According to the theoretical introduction in Section 2, these two facts involve the variance and covariance of the background error.
Formula ( 3) is taken to calculate the background error variances for 3DVAR.Removing H T from the formula (6) yields the formula to estimate the background error variances for EnKF: Forecast differences of a half-month-long term are used in formula ( 3), but instead formula (18) employs member differences at a certain time.Figure 13 shows the background error variances of temperature, relative humidity, and wind speed at 1000 hPa for the case C-09. Figure 14 is as the same as Figure 13 but for the case C-14, and the variances are calculated in the last cycle.Note that the variances for 3DVAR do not vary within the whole assimilation window.Overall, the variances of 3DVAR are larger than EnKF over sea (cf. the upper row and the bottom row in Figures 13 and 14).Since observation error can be considered fixed, larger variances make the assimilation result approach closer to observation than to background.This may explain the first fact mentioned above: why the variable increments of 3DVAR are generally larger than those of EnKF in the assimilation window (see the bars in Figure 12).
Atmosphere 2018, 9, x 20 of 29 theoretical introduction in Section 2, these two facts involve the variance and covariance of the background error.Formula ( 3) is taken to calculate the background error variances for 3DVAR.Removing  T from the formula (6) yields the formula to estimate the background error variances for EnKF: Forecast differences of a half-month-long term are used in formula ( 3), but instead formula (18) employs member differences at a certain time.Figure 13 shows the background error variances of temperature, relative humidity, and wind speed at 1000 hPa for the case C-09. Figure 14 is as the same as Figure 13 but for the case C-14, and the variances are calculated in the last cycle.Note that the variances for 3DVAR do not vary within the whole assimilation window.Overall, the variances of 3DVAR are larger than EnKF over sea (cf. the upper row and the bottom row in Figures 13 and 14).Since observation error can be considered fixed, larger variances make the assimilation result approach closer to observation than to background.This may explain the first fact mentioned above: why the variable increments of 3DVAR are generally larger than those of EnKF in the assimilation window (see the bars in Figure 12).To clearly reveal the impacts of the background error variances and covariances from 3VAR and EnKF, six assimilation experiments with single synthetic observation were conducted in the last cycle.Because the u-component wind increments in the significant zone of the case C-09 are larger in Exp-A3 than Exp-A1, the single observation is set to u-component wind.In the significant zone of the case C-14, although the increments of both temperature and mixing ratio are much more obvious in Exp-B3 than Exp-B1, temperature is chosen as the single observation due to the fact that mixing ratio is not an observational variable.Their description is given in Table 6.Two points (points A and C) are deliberately arranged located in the significant zones for the two cases, respectively.In addition, another point (point B) is placed in the Taiwan Strait to further demonstrate the flow-dependent feature of EnKF.The values of the observations are given based on the real relation between the observation and the background.Figure 15 shows that the different responses of 3DVAR and EnKF in assimilating a single synthetic u-component wind observation of a value of 2 m/s above the background flow (arrows in Figure 15).The increment of EnKF is smaller than that of 3DVAR due to the difference of background error variances.The former is about 0.8 m/s while the latter is about 1.2 m/s at the point A (the top right point in Figure 15).The corresponding background error variances of 3DVAR are larger than that of EnKF, which can be seen by comparing the background error variances in the significant zones of Figure 13e,f.Because the background flow is easterly, a smaller u-component To clearly reveal the impacts of the background error variances and covariances from 3VAR and EnKF, six assimilation experiments with single synthetic observation were conducted in the last cycle.Because the u-component wind increments in the significant zone of the case C-09 are larger in Exp-A3 than Exp-A1, the single observation is set to u-component wind.In the significant zone of the case C-14, although the increments of both temperature and mixing ratio are much more obvious in Exp-B3 than Exp-B1, temperature is chosen as the single observation due to the fact that mixing ratio is not an observational variable.Their description is given in Table 6.Two points (points A and C) are deliberately arranged located in the significant zones for the two cases, respectively.In addition, another point (point B) is placed in the Taiwan Strait to further demonstrate the flow-dependent feature of EnKF.The values of the observations are given based on the real relation between the observation and the background.Figure 15 shows that the different responses of 3DVAR and EnKF in assimilating a single synthetic u-component wind observation of a value of 2 m/s above the background flow (arrows in Figure 15).The increment of EnKF is smaller than that of 3DVAR due to the difference of background error variances.The former is about 0.8 m/s while the latter is about 1.2 m/s at the point A (the top right point in Figure 15).The corresponding background error variances of 3DVAR are larger than that of EnKF, which can be seen by comparing the background error variances in the significant zones of Figure 13e,f.Because the background flow is easterly, a smaller u-component wind increment of EnKF makes a stronger easterly wind than that of 3DVAR, which may interpret the wind increments in the significant zone in Figure 9f.Almost the same result is shown at the point B. Beyond that, the flow-dependent characteristics of EnKF are clearly found in Figure 15.Even though the observation measures the values of the u-component at one point, the data assimilation spreads out this information.However, the spreads of 3DVAR and EnKF differ very much.For 3DVAR, the spread patterns at points A and B are almost the same though the background flow is markedly different (Figure 15a).However, for EnKF, the spread patterns are closely correlated with the flow resulting from its flow-dependent background error covariances (Figure 15b).
wind increment of EnKF makes a stronger easterly wind than that of 3DVAR, which may interpret the wind increments in the significant zone in Figure 9f.Almost the same result is shown at the point B. Beyond that, the flow-dependent characteristics of EnKF are clearly found in Figure 15.Even though the observation measures the values of the u-component at one point, the data assimilation spreads out this information.However, the spreads of 3DVAR and EnKF differ very much.For 3DVAR, the spread patterns at points A and B are almost the same though the background flow is markedly different (Figure 15a).However, for EnKF, the spread patterns are closely correlated with the flow resulting from its flow-dependent background error covariances (Figure 15b).Moisture increments in the assimilation process especially benefit from the multivariate correlations of EnKF, which was well certified by the result of Exp-B3.Here, this benefit is simply illustrated in Figure 16.Given a single synthetic temperature with a value of 2 K below the background at point C, 3DVAR only responds to the temperature itself (Figure 16a), however the response of EnKF includes not only the temperature but also the mixing ratio (Figure 16b).The amplitude pattern of the temperature increments is almost the same as that of the background error variances of temperature (cf. the color shadings in Figure 16b,c), indicating that the temperature assimilation strongly depends on the background error variances of temperature.However, the assimilation of mixing ratio is chiefly controlled by the background error covariances between relative humidity and temperature (cf. the color shading in Figure 16d and the contours in f Figure 16b).Moisture increments in the assimilation process especially benefit from the multivariate correlations of EnKF, which was well certified by the result of Exp-B3.Here, this benefit is simply illustrated in Figure 16.Given a single synthetic temperature with a value of 2 K below the background at point C, 3DVAR only responds to the temperature itself (Figure 16a), however the response of EnKF includes not only the temperature but also the mixing ratio (Figure 16b).The amplitude pattern of the temperature increments is almost the same as that of the background error variances of temperature (cf. the color shadings in Figure 16b,c), indicating that the temperature assimilation strongly depends on the background error variances of temperature.However, the assimilation of mixing ratio is chiefly controlled by the background error covariances between relative humidity and temperature (cf. the color shading in Figure 16d and the contours in f Figure 16b).

Sensitivity of Localization Scale and Inflation Factor
To explore the impacts of localization scale and inflation factor on the performance of EnKF for the two sea fog cases, a series of EnKF assimilation and forecast experiments (hereafter sensitivity experiments) were conducted.These experiments are the same as Exp-A3/Exp-B3, except for different localization scales and inflation factors that are taken from the combination of five scale values (50 km, 100 km, 200 km, 500 km, and 800 km) and five factor values (0.6 to 1.0 with 0.1 interval).Thus, for each sea fog case, there are 25 experiments to conduct, among which one experiment is Exp-A3 or Exp-B3 which is called the control experiment and needs not to be done.
Figures 17 and 18 show the differences of the initial conditions between the sensitivity experiments and Exp-A1 for the case C-09.For a given localization scale (L), e.g., 50 km, the differences of mixing ratio are almost the same with different inflation factors (see the panels in the extreme left column in Figure 17).However, for an inflation factor (β), the differences vary much with different localization scales.The result of differences for SLP and 10-m wind is quite similar (Figure 18).It indicates that the localization scale likely dominates the impact on the performance of EnKF.The same conclusion can also be drawn for the case C-14 (figures not shown).

Sensitivity of Localization Scale and Inflation Factor
To explore the impacts of localization scale and inflation factor on the performance of EnKF for the two sea fog cases, a series of EnKF assimilation and forecast experiments (hereafter sensitivity experiments) were conducted.These experiments are the same as Exp-A3/Exp-B3, except for different localization scales and inflation factors that are taken from the combination of five scale values (50 km, 100 km, 200 km, 500 km, and 800 km) and five factor values (0.6 to 1.0 with 0.1 interval).Thus, for each sea fog case, there are 25 experiments to conduct, among which one experiment is Exp-A3 or Exp-B3 which is called the control experiment and needs not to be done.
Figures 17 and 18 show the differences of the initial conditions between the sensitivity experiments and Exp-A1 for the case C-09.For a given localization scale (L), e.g., 50 km, the differences of mixing ratio are almost the same with different inflation factors (see the panels in the extreme left column in Figure 17).However, for an inflation factor (β), the differences vary much with different localization scales.The result of differences for SLP and 10-m wind is quite similar (Figure 18).It indicates that the localization scale likely dominates the impact on the performance of EnKF.The same conclusion can also be drawn for the case C-14 (figures not shown).In order to directly show the impacts of localization scale and inflation on the sea fog forecast, statistical results of the ETS score were calculated for all sensitivity experiments.Tables 7 and 8 list the ETS score for each sensitivity experiment and its improvement relative to the control experiment (Exp-A3 or Exp-B3).Table 7 shows that the ETS score with the smaller localization scale is better than that with the larger scale.However, the result in Table 8 is reversed.The ETS score with the smaller localization scale is worse than that with the larger one.For example, the ETS scores of the sensitivity experiments with the 50 km scale compared to those with the 200 km scale, are improved by at least 26% for the case C-09 (Table 7) while are decreased by at least 20% for the case C-14.Note that 200 km is apparently the dividing scale.The best localization scale is 50 km for the case C-09, while 500 km is the best for the other case.
Unlike the obvious impact of the localization scale, the inflation factor affects the ETS score slightly, but also a best value exists: 0.7 for the case C-09 and 1.0 for the case C-14.The results of ETS scores, the initial differences in Figures 17 and 18 as well, suggest that localization scale and inflation are important for the EnKF performance in sea fog modeling and they should be tuned for particular cases.In order to directly show the impacts of localization scale and inflation on the sea fog forecast, statistical results of the ETS score were calculated for all sensitivity experiments.Tables 7 and 8 list the ETS score for each sensitivity experiment and its improvement relative to the control experiment (Exp-A3 or Exp-B3).Table 7 shows that the ETS score with the smaller localization scale is better than that with the larger scale.However, the result in Table 8 is reversed.The ETS score with the smaller localization scale is worse than that with the larger one.For example, the ETS scores of the sensitivity experiments with the 50 km scale compared to those with the 200 km scale, are improved by at least 26% for the case C-09 (Table 7) while are decreased by at least 20% for the case C-14.Note that 200 km is apparently the dividing scale.The best localization scale is 50 km for the case C-09, while 500 km is the best for the other case.Unlike the obvious impact of the localization scale, the inflation factor affects the ETS score slightly, but also a best value exists: 0.7 for the case C-09 and 1.0 for the case C-14.The results of ETS scores, the initial differences in Figures 17 and 18 as well, suggest that localization scale and inflation are important for the EnKF performance in sea fog modeling and they should be tuned for particular cases.

Conclusions
Numerical forecasting of sea fog is undoubtedly challenging.It can be tricky due to many imperfect aspects closely related to sea fog simulation, such as the microphysics scheme for vapor condensation with aerosol effect, the turbulence scheme for air-sea interaction in the bottom of marine PBL, the subsidence inversion caused by a high, warm-moisture advection controlled by a cyclone/anticyclone couplet, and the initial conditions.Particularly, the last reason is vital for sea fog modeling.Some work was done previously for it based on the 3DVAR method.
In view of the theoretical advantages of EnKF and its effective applications on some weather phenomena except sea fog, two sea fog cases over the Yellow sea were studied in detail.There are distinguishable differences between the two sea fogs, one spread widely and the other spread narrowly along the coastal area.A series of data assimilation and forecast experiments were conducted aiming at a comparison between 3DVAR and EnKF for the data assimilation effects on the Yellow Sea fog forecast.The results of the experiments, including forecasted sea fog area, vertical profiles of temperature, and moisture in the initial conditions, were evaluated by retrieved sea fog patch and measurements.The main conclusions of this research work are summarized below: (a) The assimilation effect of EnKF obviously excels that of 3DVAR, performing not only forecasted sea fog but also the distribution of temperature, moisture, and wind in the initial conditions.For the widespread-fog case, the assimilation by EnKF significantly improves the forecasted sea fog area, raising POD and ETS by up to about 57.9% and 55.5%, respectively.Especially for the case that spreads along the coast, the assimilation by EnKF successfully produces the sea fog formation which is completely mis-forecasted by 3DVAR.(b) The analysis increments strongly depend on the background error covariances.The flow-dependent background error covariances of EnKF out-compete that of 3DVAR, as evidenced by more realistic depiction of sea surface wind for the widespread-fog case and better existence of moisture for the other case in the initial conditions.
(c) Compared with 3DVAR, the multivariate correlations (e.g., correlation between temperature and humidity) in the background error covariances of EnKF play a key role in adjusting/generating moisture through assimilation of temperature.This helps greatly to improve the moisture conditions for sea fog forecast.(d) A series of sensitivity experiments of EnKF shows that the forecast result is sensitive to ensemble inflation and localization factors, in particular, highly sensitive to the latter.These factors need to be tuned case by case or an adaptive localization method should be considered.
Further studies are still needed, though encouraging results have been achieved in this study.For instance, perturbation of sea surface temperature might be considered in the EnKF ensemble; the reason needs to be sought as to why in the initial condition of the widespread-fog case, no sea fog is formed responding to the assimilation by EnKF as that in the other case.It is not enough to clarify in detail the advantages of EnKF on sea fog assimilation by only using two examples of sea fog cases, hence more cases are required in the next study.In addition, the forecast experiments with EnKF assimilation are deterministic, in which the mean of the ensemble of assimilation analyses is taken as the initial condition.As a matter of fact, it is easy and convenient to carry out ensemble forecast of sea fog initialized by the members of EnKF, and therefore the sea fog ensemble forecast based on EnKF will be desirable.

Figure 1 .
Figure 1.Korea Meteorological Administration (KMA) surface synoptic charts (a,c) and Multifunctional Transport Satellite (MTSAT) visible satellite images (b,d).Upper panels (a,b) are the case C-09 at t at 0000 UTC 10 April 2009, and lower panels (c,d) are for the case C-14 at 0300 UTC 2 April 2014.

Figure 1 .
Figure 1.Korea Meteorological Administration (KMA) surface synoptic charts (a,c) and Multifunctional Transport Satellite (MTSAT) visible satellite images (b,d).Upper panels (a,b) are the case C-09 at t at 0000 UTC 10 April 2009, and lower panels (c,d) are for the case C-14 at 0300 UTC 2 April 2014.

Figure 2 .
Figure 2. Weather Research and Forecasting (WRF) domains for the cases: (a) C-09 and (b) C-14.Colors show sea surface temperature (SST) within the D2 domains, and the red and black dots respectively locate radiosonde stations and ship measurement sites.The sea area surrounded by thick dashed lines and coastlines in D2 of (a) is used for evaluating forecasted sea surface winds later.

Figure 2 .
Figure 2. Weather Research and Forecasting (WRF) domains for the cases: (a) C-09 and (b) C-14.Colors show sea surface temperature (SST) within the D2 domains, and the red and black dots respectively locate radiosonde stations and ship measurement sites.The sea area surrounded by thick dashed lines and coastlines in D2 of (a) is used for evaluating forecasted sea surface winds later.

Figure 4 .
Figure 4. Comparison between the forecasted and observed sea fog areas for the case C-09.Panels in the top row (a-e) are the observed fog patches, and panels in the other rows from up to down (f-t) are the forecasted fog patches for Exp-A1, Exp-A2, and Exp-A3, respectively.

Figure 4 .
Figure 4. Comparison between the forecasted and observed sea fog areas for the case C-09.Panels in the top row (a-e) are the observed fog patches, and panels in the other rows from up to down (f-t) are the forecasted fog patches for Exp-A1, Exp-A2, and Exp-A3, respectively.

Figure 6 .
Figure 6.Vertical profiles of the root-mean-square errors (RMSEs) (solid lines) and biases (dashed lines) between the initial analysis and the radiosonde observations of geopotential height (a,d), temperature (b,e), and mixing ratio (c,f).The upper and lower panels are for the cases C-09 and C-14, respectively.

Figure 6 .
Figure 6.Vertical profiles of the root-mean-square errors (RMSEs) (solid lines) and biases (dashed lines) between the initial analysis and the radiosonde observations of geopotential height (a,d), temperature (b,e), and mixing ratio (c,f).The upper and lower panels are for the cases C-09 and C-14, respectively.

Figure 7 .
Figure 7. Evaluation result of sea level pressure (SLP) in the initial analysis for the Case C-09.Top panels respectively show (a) KMA surface synoptic chart, (b) RMSE, and (c) bias of SLP by using surface measurements.Bottom panels illustrate SLP of (d) Exp-A1, (e) Exp-A2, and (f) Exp-A3, respectively.

Figure 7 .
Figure 7. Evaluation result of sea level pressure (SLP) in the initial analysis for the Case C-09.Top panels respectively show (a) KMA surface synoptic chart, (b) RMSE, and (c) bias of SLP by using surface measurements.Bottom panels illustrate SLP of (d) Exp-A1, (e) Exp-A2, and (f) Exp-A3, respectively.

Figure 8 .
Figure 8. Biases between the forecasted and observed mixing ratio for the cases (a) C-09 and (b) C-14.Numbers in parentheses are the mixing ratios of ship measurements.The improvements (%) of Exp-A3 relative to Exp-A1 are in parentheses below the S# (identified mark of ship) along horizontal coordinates, as well as Exp-B3 relative to Exp-B1.

Figure 8 .
Figure 8. Biases between the forecasted and observed mixing ratio for the cases (a) C-09 and (b) C-14.Numbers in parentheses are the mixing ratios of ship measurements.The improvements (%) of Exp-A3 relative to Exp-A1 are in parentheses below the S# (identified mark of ship) along horizontal coordinates, as well as Exp-B3 relative to Exp-B1.

Atmosphere 2018, 9 , x 16 of 29 Figure 9 .
Figure 9.Initial differences of temperature (a,d) and mixing ratio (b,e) at 1000 hPa, SLP and 10-m wind (c,f) between Exp-A2 and Exp-A1 (a-c), Exp-A3 and Exp-A1 (d-f), respectively.The dashed frame indicates the significant zone, and the closed thick solid line denotes the forecasted sea fog at 1200 UTC 10 April 2009.

Figure 9 .
Figure 9.Initial differences of temperature (a,d) and mixing ratio (b,e) at 1000 hPa, SLP and 10-m wind (c,f) between Exp-A2 and Exp-A1 (a-c), Exp-A3 and Exp-A1 (d-f), respectively.The dashed frame indicates the significant zone, and the closed thick solid line denotes the forecasted sea fog at 1200 UTC 10 April 2009.

Figure 11 .
Figure 11.Comparison of observation innovation ( −   ) and analysis increment (  −   ) at the lowest model level between Exp-B1 (a-c) and Exp-B3 (d-f) for the last DA cycle.Observation innovations of temperature and mixing ratio before and after assimilation are illustrated in the left and middle panels, respectively; and analysis increments of temperature (colors) and mixing ratio (contours, unit: g/kg) are shown in the right panels.The dots indicate the observation locations, and the dashed frame indicates the significant zone.

Figure 11 .
Figure 11.Comparison of observation innovation (y − Hx b ) and analysis increment (x a − x b ) at the lowest model level between Exp-B1 (a-c) and Exp-B3 (d-f) for the last DA cycle.Observation innovations of temperature and mixing ratio before and after assimilation are illustrated in the left and middle panels, respectively; and analysis increments of temperature (colors) and mixing ratio (contours, unit: g/kg) are shown in the right panels.The dots indicate the observation locations, and the dashed frame indicates the significant zone.

Figure 12 .
Figure 12.Time series of the average (a) temperature, (b) mixing ratio, (c) u, and (d) v wind components over the area marked by the dashed frame in Figure 8 for Exp-A1 (red) and Exp-A3 (blue).The bars show the analysis increments within data assimilation window.

Figure 12 .
Figure 12.Time series of the average (a) temperature, (b) mixing ratio, (c) u, and (d) v wind components over the area marked by the dashed frame in Figure 8 for Exp-A1 (red) and Exp-A3 (blue).The bars show the analysis increments within data assimilation window.

Figure 13 .
Figure 13.Comparison between the variance distributions at 1000 hPa for temperature (left), relative humidity (middle) and wind speed (right) for 3DVAR (upper; Exp-A1) and EnKF (lower; Exp-A3) of the case C-09.The dashed frame indicates the significant zone.

Figure 13 .
Figure 13.Comparison between the variance distributions at 1000 hPa for temperature (a,d), relative humidity (b,e) and wind speed (c,f) for 3DVAR (upper; Exp-A1) and EnKF (lower; Exp-A3) of the case C-09.The dashed frame indicates the significant zone.

Figure 14 .
Figure 14.(a-f) As in Figure 13, but for Exp-B1 and Exp-B3 of the case C-14.

Figure 14 .
Figure 14.(a-f) As in Figure 13, but for Exp-B1 and Exp-B3 of the case C-14.

Figure 15 .
Figure 15.Comparison between (a) 3DVAR and (b) EnKF for u wind component (colors) increments at 1000 hPa in single synthetic observation tests from Exp-AS1 to Exp-AS4.Vectors show background flow.

Figure 15 .
Figure 15.Comparison between (a) 3DVAR and (b) EnKF for u wind component (colors) increments at 1000 hPa in single synthetic observation tests from Exp-AS1 to Exp-AS4.Vectors show background flow.

Figure 16 .
Figure 16.Results of (a) Exp-BS1 and (b-d) Exp-BS2.Colors in (a,b) show temperature increments (K) and contours in (b) represent mixing ratio increments (g/kg); Covariances between temperature at the dot and (c) temperature or (d) RH at each model grid are illustrated.

Figure 16 .
Figure 16.Results of (a) Exp-BS1 and (b-d) Exp-BS2.Colors in (a,b) show temperature increments (K) and contours in (b) represent mixing ratio increments (g/kg); Covariances between temperature at the dot and (c) temperature or (d) RH at each model grid are illustrated.

Atmosphere 2018, 9 , x 24 of 29 Figure 17 .
Figure 17.(a-y) Initial differences of mixing ratio at 1000 hPa between the sensitivity experiments and Exp-A1 for the case C-09.The dashed frame indicates the significant zone.

Figure 17 .
Figure 17.(a-y) Initial differences of mixing ratio at 1000 hPa between the sensitivity experiments and Exp-A1 for the case C-09.The dashed frame indicates the significant zone.

Figure 18 .
Figure 18.(a-y) As in Figure 17, but for SLP and 10-m wind.

Figure 18 .
Figure 18.(a-y) As in Figure 17, but for SLP and 10-m wind.

Table 4 .
Design of experiments.

Table 5 .
Statistical results of the experiments.The improvements (%) in Exp-A2 and Exp-A3 relative to Exp-A1 are in parentheses and set in boldface, as well as Exp-B2 and Exp-B3 relative to Exp-B1.

Table 5 .
Statistical results of the experiments.The improvements (%) in Exp-A2 and Exp-A3 relative to Exp-A1 are in parentheses and set in boldface, as well as Exp-B2 and Exp-B3 relative to Exp-B1.

Table 6 .
Description of single observation tests.

Table 6 .
Description of single observation tests.

Table 7 .
Statistical results of the equitable threat score (ETS) for the sensitivity experiments for the case C-09.Their improvements (%) relative to Exp-A3 are in parentheses and set in boldface.

Table 7 .
Statistical results of the equitable threat score (ETS) for the sensitivity experiments for the case C-09.Their improvements (%) relative to Exp-A3 are in parentheses and set in boldface.

Table 8 .
As in Table7, but for the sensitivity experiments for the case C-14 relative to Exp-B3.