Impact of Multivariate Background Error Covariance on the WRF-3DVAR Assimilation for the Yellow Sea Fog Modeling

Numerical modeling of sea fog is highly sensitive to initial conditions, especially to moisture in the marine atmospheric boundary layer (MABL). Data assimilation plays a vital role in the improvement of initial MABL moisture for sea fog modeling over the Yellow Sea. In this study, the weather research and forecasting (WRF) model and its three-dimensional variational (3DVAR) data assimilation module are employed for sea fog simulations. Two kinds of background error (BE) covariances with diﬀerent control variables (CV) used in WRF-3DVAR, that is, CV5 and multivariate BE (CV6), are compared in detail by explorative case studies and a series of application experiments. Statistical veriﬁcation metrics including probability of detection (POD) and equitable threat scores (ETS) of forecasted sea fog area are computed and compared for simulations with the implementations of CV5 and CV6 in the WRF-3DVAR system. The following is found: (1) there exists a dominant negative correlation between temperature and moisture in CV6 near the sea surface, which makes it possible to improve the initial moisture condition in the MABL by assimilation of observed temperature; (2) in general, the performance of the WRF-3DVAR assimilation with CV6 is distinctly better, and the results of 10 additional sea fog cases clearly suggest that CV6 is more suitable than CV5 for sea fog modeling. Compared to those with CV5, the average POD and ETS of forecasted sea fog area using 3DVAR with CV6 can be improved by 27.6% and 21.0%, respectively.


Introduction
Sea fog occurs within the marine atmospheric boundary layer (MABL). It is a type of fog that forms over an ocean surface and spreads over coastal areas. It can seriously affect many maritime activities (e.g., transport, navigation, fishery, and rescue) due to its low horizontal visibility [1,2]. e Yellow Sea has the largest number of foggy days among the four seas of China [3][4][5], with a long fog season starting from March and ending in September. e Yellow Sea is located to the north of the Kuroshio Current and under influences of warm-moist air mass that is transported northward as the East Asia summer monsoon gradually intensifies. Due to the cold sea surface temperature (SST), an atmospheric inversion layer can easily develop in the MABL over the Yellow Sea, leading to the formation of advection sea fog in the inversion layer [3]. erefore, the sensitivity of sea fog modeling to the initial condition is actually largely related to temperature and moisture in the MABL [6][7][8][9][10][11][12].
After the textbook Sea Fog published in 1985 by Wang [4], Koračin and Dorman [2] published a book in 2017, in which they summarize the challenges and advancements in observations, modeling, and forecasting of marine fog over the world. Modeling and forecasting have become an important research content of sea fog. In order to mitigate socioeconomic impacts of sea fog over the Yellow Sea and meet increasing practical demands, it is urgent to improve the forecast skill of sea fog. Undoubtedly, a high-quality initial condition is very beneficial to the numerical forecasting for sea fog. A vital way to improve the initial condition in numerical modeling is data assimilation [13][14][15][16][17][18][19]. For the assimilation of given observations, different data assimilation methods yield different assimilation effects. In general, better effect can be obtained when using more advanced data assimilation method. For instance, in their modeling study of sea fog using the weather research and forecasting (WRF) model, Gao et al. [20] demonstrated that ensemble Kalman filter (EnKF) performs better than threedimensional variational (3DVAR) assimilation, which is attributed to the flow-dependent background error (BE) covariance used in the EnKF. However, the EnKF requires much more computing resources compared to the 3DVAR, which is an issue that has to be considered for real-time forecast of sea fog. e WRF-3DVAR [18] has been widely used for sea fog modeling over the Yellow Sea due to its low computational cost [21][22][23][24].
3DVAR minimizes the analysis error variance by minimizing the cost function, which is defined by the weighted sum of analysis increment and observation innovation [25], respectively, with BE and observational error (OE) covariance. A better estimate of BE can definitely lead to a better assimilation analysis.
us, how to conduct reasonable statistics of BE is an important issue, which has spawned a lot of related research works [26][27][28][29][30][31]. Currently, the National Meteorology Center (NMC; now called the National Centers for Environmental Prediction, NCEP) method [32] is the most commonly used to generate BE for the 3DVAR. It first conducts transformation of model variables to a set of control variables and then estimates the relationships between different variables by linear regression. e NMC method is implemented in the WRF-3DVAR [18]. ere are two types of BE with different control variables (CV): one defines the control variables in the physical space (CV3), and the other defines the control variables in eigenvector space (CV5, CV6, and CV7). CV3 uses a vertical recursive filter to model the vertical covariance, while the others use an empirical orthogonal function (EOF) to represent the vertical covariance. CV3 is a global BE and can be used for any regional domain, while CV5, CV6, and CV7 are domaindependent and therefore must be generated based on forecast or ensemble data over the same domain. Due to the poor performance of CV3, it is not suitable for sea fog modeling. In CV5, the moisture control variable is not related to any other variables. In CV6, the moisture control variable is the unbalanced portion of the pseudo-relative humidity, and six additional correlation coefficients in the definition of the balanced part of analysis control variables are introduced. Namely, CV6 is a multivariate background error covariance, especially for the moisture. CV7 is similar to CV5, but it uses a different set of control variables and there exist no correlations between the moisture and other control variables.
e issue on the strengths and shortcomings of CV5 and CV6 has been discussed. Based on the work to use multiple linear regressions to estimate multivariate analysis of humidity in a limited-area model, Berre [33] pointed that the relationships between forecast errors of humidity and those of mass and wind fields cannot be negligible. Dhanya and Chandrasekar [34,35] found that the forecasts of rainfall caused by monsoon depressions and tropical cyclones can be improved by applying data assimilation with CV6, which is better than CV5. However, among existing research work on sea fog over the Yellow Sea involving numerical modeling [24,[36][37][38][39][40], quite a few studies, especially those works with focus on the mechanism of sea fog evolution, have little or no mention of data assimilation. Even if data assimilation is employed in the numerical simulations for sea fog, CV5 is generally used [22,23,[41][42][43]. e reason might be that CV5 is the default option except CV3 in the WRF-3DVAR, which results in the fact that little attention has been paid to CV6 in the application of data assimilation for sea fog modeling.
As stated previously, a successful sea fog modeling definitely requires the assimilation of observations over the sea [22]. Unfortunately, in spite of a large amount of satellite-retrieved data available over the sea, the accuracy of the retrieved water vapor profiles is unsatisfactory. Divakarla et al. [44] performed an evaluation of satellite-retrieved data using global radiosonde observations and found that the root mean square difference is larger than 15% for relative humidity at 1000 hPa. Such big observational errors lead to a malfunction of the retrieved humidity in data assimilation. e approach to correct errors in the observed humidity field and make it more realistic by increments of other control variables is probably an efficient way to modify the MBL moisture structure. Hence, the purpose of the present study is to investigate the impact of multivariate background error covariance on the WRF-3DVAR assimilation for fog modeling over the Yellow Sea. e assimilation effects with CV5 and CV6 are examined. e remaining content is arranged as follows. Section 2 briefly describes the multivariate BE of the WRF-3DVAR and the data assimilation scheme, and Section 3 outlines numerical experiments, including data, sea fog cases, model configurations, and numerical experiments. Results, analysis, and evaluation are presented in detail in Section 4. Finally, summary and conclusions are provided in Section 5.

Multivariate BE.
e cost function for 3DVAR is defined by where x and x b are, respectively, the model state and the background field, y is the observation, and H is the observation operator matrix that transforms data from model space to observation space. In equation (1), analysis increment (x-x b ) and observation innovation (y-Hx) are, respectively, weighed by BE (B) and OE (R). OE is determined by instrument errors and representation errors caused by the transformation of data from model space to observation space; BE is determined by model forecast errors. e analysis x a represents a posterior maximum likelihood estimate of the true state of the atmosphere, and it can be obtained if the cost function J(x) is minimized. e WRF-2 Advances in Meteorology 3DVAR system uses Quasi-Newton method [45] to realize the minimization process. Assuming that the model has m degrees of freedom and the assimilated observation number is p, x and x b are m-dimensional vectors, y is an p-dimensional vector, H is a p × m-dimensional matrix, B is a m × m-dimensional matrix, and R is a p × p-dimensional matrix that is diagonal because observational errors are usually supposed to be uncorrelated. In practice, the amount of calculation of J b is much larger than that of J o , since R is diagonal and p is far less than m.
Usually, a numerical model has m degrees of freedom with a typical value of 10 7 ; thus it is almost impossible to explicitly solve J b , which requires ∼O(m 2 ) calculations. However, in the model space, the correlations between different variables or between the same variables at different grid points are available only within a limited range, and thus BE is a sparse matrix. As a result, by decomposing BE, a lot of computational cost can be greatly reduced.
is process is known as control variable transform; that is, the model variables x are transformed to control variables v, and J b can then be calculated by v instead of by Note that U is the transform matrix defined in a way that satisfies B � UU T . Using the increment formulation [46], equation (1) can be rewritten as e WRF-3DVAR system takes three steps to complete the control variable transform: the horizontal transform U h , the vertical transform U v , and the physical transform U p . erefore, U � U p U v U h . U h uses recursive filter [47,48] to calculate the correlations between data on different grid points at the same level, U v implements empirical orthogonal function (EOF) decomposition to obtain the correlations between data at different model levels, and U p does linear regression to get the correlations between different variables.
In the physical transform x ′ � U p v p , the increments of model variables, that is, zonal wind u, meridional wind v, temperature T, pressure P, and mixing ratio Q, are first converted to a new set of variables, that is, stream function ψ, velocity potential χ, temperature T, surface pressure P s , and pseudo-relative humidity RH s . ese new variables are then divided into balanced and unbalanced components: where φ can be any variable. e balanced part φ b is related to other control variables μ by linear regression α μφ ; that is, φ b � α μφ μ. e unbalanced part φ u is independent and serves as input to the vertical transform. By the above procedures, the correlated model variables are transformed into uncorrelated variables.
In the WRF-3DVAR, the correlations are available only within a subset of the variables. CV5 has three basic coefficients: α ψχ , α ψT , and α ψP s , and these coefficients link χ, T, and P s to ψ: where i and j denote the horizontal dimension index, k and l indicate the vertical levels, n is the total number of vertical levels, and the regression coefficients are usually supposed to be horizontally homogeneous. Compared to CV5, the multivariate BE (CV6) contains six extra coefficients, that is, α χ u T , α χ u P s , α ψRH s , α χ u RH s , α T u RH s and α P su RH s . RH s also consists of balanced and unbalanced components. e control variable transform can be expressed by Formula (4d) indicates that an increment of four other variables can lead to an increment of RH s in CV6. As mentioned in Introduction, moisture information in the initial condition is vital for sea fog modeling. Note that the analysis of moisture is retrieved from the analysis of RH s , background temperature, and pressure together. e moisture increment is directly proportional to RH s increment because background data are recognized as constant.

Data Assimilation Scheme.
A cycling-3DVAR scheme [21] is employed to do data assimilation (DA) in this study. e scheme is shown in Figure 1. CV5/6 presents CV5 or CV6, and "obs" stands for observations. Individual assimilation cycles are connected by the WRF integration (i.e., wrf.exe), and the final analysis x a is taken as the initial condition for the WRF simulation. In this study, all WRF Advances in Meteorology simulations are designed with two domains in two-way nesting.
As shown in Figure 1, the cycling-3DVAR scheme includes three 3DVAR updates. Each update is only performed on the outer WRF model domain with its BE, and initial condition of the inner domain is obtained by interpolation of analysis field for the outer domain after its 3rd update. In the case study, the data assimilation period is 6 hours with a DA interval of 3 hours. e NMC method [32] is used to generate CV5 and CV6, assuming that the BEs are well approximated by the averaged differences between 12-hour forecasts and 24-hour forecasts, which are expressed by where x 12 and x 24 are 12-hour and 24-hour forecasts, respectively. ese forecasts can be obtained from a period of time (e.g., 15−30 days) or can be got from an ensemble prediction system. e former way is adopted here, and thus equation (5) is rewritten as where T n is the time period for statistics. Here, the forecasts are collected within 15 days centered at the initial day of the case simulation.

Numerical Experiments
3.1. Data. Initial and lateral boundary conditions for numerical modeling were derived from the NCEP Final Analysis (FNL) Data (1°× 1°; 6 hourly). Bottom boundary conditions were extracted from the daily SST dataset (0.25°× 0.25°) of the North-East Asian Regional Global Ocean Observing System (NEAR-GOOS). Observational data for the 3DVAR data assimilation include conventional surface and upper air measurements from the Global Telecommunication System (GTS), buoy and island measurements, and AIRS (Atmospheric Infrared Sounder) retrieved temperature and humidity profiles. Surface synoptic maps were downloaded from the Korea Meteorological Administration (KMA) to show weather situations of the sea fog cases in the present study. e observed fog areas were empirically retrieved following the approach of Wang et al. [22], which uses the multichannel data of Japanese Multifunctional Transport Satellite (MTSAT), including albedo, infrared, and visible cloud imageries. e NCEP ship observations in PrepBUFR format, which can be downloaded from the NCAR Computational and Information Systems Laboratory (CISL) Research Data Archive, were used to evaluate model results.

Sea Fog Cases.
Two typical sea fog cases over the Yellow Sea, one spreading narrowly over the eastern Yellow Sea along the Korean Peninsula coast (Figure 2(b)) and the other occupying a broad area over the southern Yellow Sea (Figure 3 Figure 3(a)), both cases formed over the areas controlled by a high-pressure system. However, the SST conditions are different for the two cases. Case-A occurred along the west coast of the Korean Peninsula, where the SST was about 5°C (Figure 2(c)). e fog area was under the control of a weak high-pressure system with the center located at the Korean Peninsula (Figure 2(a)). Warm-moist air mass was transported from the central Yellow Sea to the cold SST area near the coast by the southwesterly wind ( Figure 2(c)). Case-B occurred over an ocean area, where a warm SST tongue and a cold SST tongue (Figure 3(c)) coexisted with large SST gradients. Such an SST distribution is usually favorable for the formation of sea fog [49]. In addition, easterly winds developed over the southern Yellow Sea (Figure 3(c)) due to the influence of the high pressure ( Figure 3(a)), transporting huge amounts of warm-moist air mass to the fog area.

Model and Its Configuration.
e Advanced Research core of the WRF (ARW, version 3.8.1) and its corresponding WRF-3DVAR were employed in the present study for numerical simulations and data assimilations. Based on the previous work related to sea fog simulation over the Yellow Sea fog [22][23][24]42], the Yonsei University (YSU) planetary boundary layer (PBL) scheme and the Lin microphysics scheme (LIN) were selected for the present study. Appropriate, vertical levels were specified as well. Detailed model configurations are listed in Table 1. Due to the distinctly different occurrence locations and areas between the two sea fog cases, the model domains and resolutions were separately designed for the two cases ( Figure 4). Two-way nesting was implemented in the simulations.  Figure 1: Flowchart of the 3DVAR data assimilation scheme. See details in the text. experiments include six experiments, which were conducted to explore the impacts of CV6 on the 3DVAR assimilation. e evaluation experiments were conducted for the purpose of evaluating the effect of CV6 application. Table 2 lists the case study experiments for Case-A and Case-B. e initial times of these simulations are presented in Table 1, and all the simulations ran out to 24 hours. e assimilation of observations collected at different platforms may have complicated superposition effects, which is not favorable for clear analysis on the assimilation impact. erefore, only the GTS routine observations were used for assimilation. For the coastal sea fog case, only observations over the southern Korean Peninsula (Figure 4(a)) were assimilated and the data collected at other far away observational sites were excluded. For the widespread sea fog case, all the GTS observations in D1 were assimilated (Figure 4(b)).

Experimental
A suite of experiments were conducted for 10 additional sea fog cases that occurred over the Yellow Sea during 2007-2014. In order to produce better initial conditions for the WRF model, the assimilation window was extended from 6 hours for the case study to 12 hours for these cases, and the observations assimilated include the GTS routine observations, buoy and ship observations and island measurements, and AIRS retrieved data. As with the BEs of those experiments listed in Table 2, CV5 and CV6 were generated in advance separately for each case using the method described in Section 2.2. e model configurations are the same as those used for the widespread sea fog case (see the description for Case-B in Table 1 and the domains in (Figure 4(b)). e experiments consist of two groups: Group-CV5 running with CV5 and Group-CV6 running with CV6, respectively.  relationships between them are decided by the regression coefficients. In order to investigate how close these variables are related, contribution rates of other control variables to RH s were calculated from CV6 generated in advance via the NMC method for Case-A.

Single-Observation
e results are shown in Figure 5(a). At low levels of the model, the balanced component accounts for 30-40% of the total RH s , among which the term T u contributes up to 20%-30% of the total.
Other individual terms can only make a contribution of less than 5%. As a result, the term T u accounts for more than    Advances in Meteorology ∼70% of the balanced component of RH s , indicating that the temperature increment is much more important for the humidity increment than the wind and pressure increments. Similar results are also found for Case-B (figure not shown). Figure 5(b) presents the distribution of regression coefficient α T u RH s (k1, k2) of T u at level k2 to calculate RH s at level k1 (see equation (4a)). e horizontal/vertical axis stands for k2/k1. e coefficients along the diagonal of Figure 5(b) are all negative, and their absolute values are much larger than other data on the same row. As a result, the humidity increment at any specific point is basically decided by temperature increment at the same point, and a negative increment of T u always causes a positive increment of RH s .

Contribution
Rate of T u to RH s . A series of singleobservation experiments were designed based on Case-A with CV5/CV6, and the single observation is marked by the red dot shown in Figure 4(a). Each experiment was conducted at different vertical model level with a temperature observation of 1 K lower than the background. ese experiments are denoted by Exp-Sk, where k (from 1 to 20) stands for the vertical level number.
Since α T u RH s (1, 1) is negative, decrease of temperature with a maximum value of about −0.3 K (Figure 6(a)) produces an increase of mixing ratio with a maximum value of about 0.06 g/kg (Figure 6(c)).
e maximum values of temperature and mixing ratio increments at levels 1-20 from Exp-S1-20 are illustrated in Figures 6(b) and 6(d). rough vertical transform, the negative temperature increment spreads to several vertical levels and causes increases of mixing ratio there. Figure 6(d) has a very similar pattern to that of Figure 6(b), which again shows that the humidity change mainly results from the change of temperature at the same level, though there exist some slight influences from neighboring levels. Compared with CV6, although similar patterns and a little bit weaker intensities of temperature assimilation result are obtained for CV5 (cf. Figure 6(e) with Figure 6(f ) and Figure 6(a) with Figure 6(b), respectively), there is not obviously any gain for mixing ratio (Figures 6(g) and 6(h)) due to the fact that there is no relationship between moisture control variable and any other variable.
Overall, for CV6, regression coefficients α T u RH s of the same level are all negative and their absolute values are usually much larger than those of neighboring levels, leading to an increase of mixing ratio when a lower temperature observation is assimilated. For 3DVAR with CV5, assimilation of temperature observations cannot change humidity analysis, since RH s is an independent variable in CV5.

Verification.
Simulated sea fog is defined to those areas where the following criteria are met [22]: cloud liquid water (CLW) at the model's lowest level is ≥ 0.016 g/kg or cloud top is ≤ 400 m. CLW ≥ 016 g/kg is equivalent to a visibility ≤ 1 km. e cloud top is calculated based on the threshold of CLW ≥ 0.016 g/kg; and observations indicate that advection sea fogs are deeper than other types of fog but rarely exceed 400 m [56]. e fog occurred over a narrow ocean area along the western coast of the Korean Peninsula (Figures 7(a)-7(e)). As mentioned in Section 3.4, only observations over the southern Korean Peninsula Figure 4(a) were assimilated, because we want to explore the impact of assimilating observations that are close to the fog area. e fog areas forecasted by Exp-A_CV5 and Exp-A_CV6 are shown in Advances in Meteorology Figure 7. It is apparent that Exp-A_CV6 performs better than Exp-A_CV5, because the former can well capture the sea fog evolution, while the latter completely fails to reproduce the sea fog (cf. Figures 7(f )-7(j) and 7(k)-7(o) with 7(a)-7(e), respectively). Note that the difference of 925 hPa wind speeds between these two experiments is less than 1 m/ s (figure not shown), which indicates that CV6 makes a trivial change for wind field. e agreement between the initial conditions and the radiosonde observations (see the big dots in Figure 4(a) for their locations) over the Korean Peninsula was examined. Figure 8 shows the average vertical profiles of bias and root mean square error (RMSE) of temperature and mixing ratio for Exp-A_CV5 and Exp-A_CV6. Below the height of 600 m, the initial temperatures of both experiments are higher than observations (Figure 8(a)), while most of the initial mixing ratios are smaller than observations (Figure 8(b)). Nevertheless, the mixing ratio bias of Exp-A_CV6 is much smaller than that of Exp-A_CV5. e same is true for the RMSE. It means that the initial moisture condition of Exp-A_CV6 is much wetter than that of Exp-A_CV5, which results in the better forecasted sea fog area by Exp-A_CV6.
In addition, ship measurements (not assimilated; see their sites in Figure 9(a) were used to verify the moisture condition near the sea surface. Biases between model output and the measurements were calculated and results are shown in Figure 9(b). It is obvious that the initial moisture conditions of both Exp-A_CV5 and Exp-A_CV6 are drier than the observed ones. However, the dry condition is significantly improved in Exp-A_CV6 compared to that in Exp-A_CV5. Note that the sites S6, S7, and S8 are close to the Korean Peninsula. Compared to Exp-A_CV5, the improvements at these three sites in Exp-A_CV6 are, respectively, 48.7%, 76.7%, and 100%, and the corresponding biases are reduced by 0.39, 0.23, and 0.47 g/kg, respectively.

Analysis Increments.
Exp-A_CV5 performs poorly. As shown in Figures 7(f )-7(j), it produces almost no sea fog during the whole forecast period. In contrast, Exp-A_CV6 successfully simulates the sea fog event (Figures 7(k)-7(o)). From Figure 9, it is found that the initial condition of Exp-A_CV5 is quite dry. For example, the moisture bias at the site of S6 is about −0.8 g/kg. Note that this dry bias is corrected in Exp-A_CV6 and the bias is reduced to 0.4 g/kg. e reason why Exp-A_CV6 can gain such a significant moisture improvement is discussed below.
During the analysis update of the WRF-3DVAR, the minimization of the cost function (equation (1)) produces the analysis field x a . x a is determined by the balance of analysis increment (x-x b ) and observation innovation (y-Hx) under the joint constraint of OE and BE. Because OE is usually prescribed regardless of a specific sea fog case, BE plays the key role. ere are three updates for x a in the designed data assimilation scheme (Figure 1), and the 3rd update yields the initial condition for the experiments. Figure 10 demonstrates analysis fields of temperature, wind, and mixing ratio at the lowest model level for three DA updates in Exp-A_CV5 (upper row) and Exp-A_CV6 (lower row). It is clearly seen that the analysis fields of both temperature and wind are almost the same between Exp-A_CV5 and Exp-A_CV6. e analysis fields of mixing ratio  Figure 10(c) with Figure 10(f ); see bold lines with the value of 5.5 g/kg). Water vapor is much less near surface over sea area west of Korean Peninsula in Exp-A_CV5 than that in Exp-A_CV6, which leads to its failure. Analysis fields at upper levels are also checked and significant difference of mixing ratio between Exp-A_CV5 and Exp-A_CV6 exists below ∼600 m (not shown). Next, we try to investigate in detail the difference between the observation innovations of Exp-A_CV5 and Exp-A_CV6 at the 3rd update, as well as the consequence of the difference on analysis increment. Figure 11 illustrates the observation innovations and analysis increments at the lowest model level for Exp-A_CV5 and Exp-A_CV6. For the temperature innovation, there is nearly no difference between Exp-A_CV5 and Exp-A_CV6 (cf. Figure 11(a) with Figure 11(e)). Meanwhile, the difference in moisture innovation is very small, and the moisture innovation in the area denoted by a red-dashed frame (Figure 11(f )) is a little bit larger in Exp-A_CV6 than in Exp-A_CV5 (Figure 11(b)). However, Exp-A_CV6 obtains much more moisture increment than Exp-A_CV5 (cf. Figure 11(c) with Figure 11(g)). Does this moisture increment still result from temperature innovation through the way revealed in the single-observation experiments of Section 4.1? To answer this question, two extra experiments were conducted, and the value of T u (see equations (3a)-(3c) and (4a)-(4d)) is set to 0 in the 3rd cycle of the WRF-3DVAR; namely, the temperature information is not used in the assimilation process. Results of the two extra experiments show that the difference in moisture increment between Exp-A_CV6 and Exp-A_CV5 is very small (cf. Figure 11(d) with Figure 11(h)), which confirms that the larger moisture increment shown in Figure 11(g) is due to the negative temperature innovation shown in Figure 11(  increments in Figure 11(g) only covers the fog area (see Figure 7). erefore, the success of Exp-A_CV6 is undoubtedly the consequence of the moisture improvement in its initial condition, and CV6 is vital to guarantee this improvement.

10
Advances in Meteorology e temperature and moisture structures of the MABL in the initial conditions of Exp-B_CV5 and Exp-B_CV6 were compared with six radiosonde soundings (red dots in Figure 4(b). Vertical profiles of RMSE and bias of temperature and mixing ratio for Exp-A_CV5 and Exp-A_CV6 are presented in Figure 13. Similar to Figure 8, the initial temperatures of both experiments are higher than observations (Figure 13(a)). However, the initial mixing ratios are smaller than observations below 300 m see the biases in (Figure 13(b)). Compared to Exp-B_CV5, the drier condition is alleviated by 0.2 g/kg in Exp-B_CV6, which results in better performance. Ship measurements (not assimilated; see their sites in (Figure 14(a)) were also employed to check the moisture condition near the sea surface. e biases between model outputs and the measurements are displayed in (Figure 14(b)), which shows clearly that Exp-B_CV6 outperforms Exp-B_CV5 at six out of eight sites and the maximum improvement is up to 36.7%.

Differences in Initial Condition. Since Exp-B_CV5 and
Exp-B_CV6 share the same model configuration, their different performances are attributed to differences between their initial conditions. Figure 15 shows the differences in initial condition near the surface between the two experiments (Exp-B_CV6 minus Exp-B_CV5). e differences in pressure and wind are not presented because they are very small. e difference in moisture is relatively larger than that in temperature. e differences in moisture and temperature over the sea fog area (see the sea fog in Figures 12(a), 12(f), and 12(k)) are quite small. For example, the difference in temperature is within the range of about −0.2-+0.2°C, and the differences in mixing ratio and relative humidity are within the ranges of about −0.2-+0.2 g/kg and 0-5%, respectively. Note that there are two areas with large positive and negative differences in mixing ratio, respectively (Figure 15(b)). e area of positive difference is located over the northwestern Yellow Sea (see the red shading around Shandong Peninsula marked in Figure 2(b)), and the area of negative difference is located over the southern Korea Peninsula. e area of positive mixing ratio difference stretches to the southeast, while positive difference in relative humidity occurs over almost the entire Yellow Sea (Figure 15(c)). Different from the positive analysis increment of mixing ratio over the Yellow Sea in Exp-B_CV6,  there is almost no positive gain of mixing ratio in Exp-B_CV5 (cf. Figures 16(a) and 16(b)). e above results demonstrate that the initial condition of Exp-B_CV6 is wetter than that of Exp-B_CV5 in the MABL, which leads to better performance.
For Case-A described in Section 4.2, the assimilation of the observations over the Korea Peninsula (see their sites in Figure 4(a) can noticeably improve the modeling result of Exp-A_CV6, because positive analysis increment of moisture is produced by negative temperature observation innovation (see Figure 11). On the contrary, for Case-B, both Exp-B_CV5 and Exp-B_CV6 get negative analysis increments of moisture over Area-K (shown in Figure 16), and the increment for Exp-B_CV6 is larger than that for Exp-B_CV5. is can be explained by the observation innovations of moisture and temperature over Area-K. e moisture innovation is negative and the temperature innovation is positive over Area-K for both Exp-B_CV5 and Exp-B_CV6 (figures not shown). For Exp-B_CV5, the negative increment of moisture is produced because RH s is an independent variable in CV5. For Exp-B_CV6, however, since there exist correlations between RH s and other variables, especially temperature (equation (4a)), negative moisture increment and positive temperature increment jointly produce a large negative increment of moisture over Area-K. In contrast, Exp-B_CV6 gets positive moisture increment (see the yellow shading in Figure 16(b)) over other areas of the Yellow Sea due to the assimilation of other observations beyond Area-K.
Forced by the sea surface wind field (see the vectors in Figure 12), the southwesterly advection pulls drier air mass from Area-K to the southern Yellow Sea, which is not conducive to the formation and development of sea fog there. Despite this unfavorable situation, Exp-B_CV6 still performs better than Exp-B_CV5, which is attributed to the fact that its initial condition contains positive moisture increment over the Yellow Sea caused by the correlation between moisture and other control variables in CV6. To further explore the impact of CV6, two sensitivity experiments-Exp-B_CV5e and Exp-B_CV6e-were, respectively, conducted based on Exp-B_CV5 and Exp-B_CV6 but without assimilation of the observations over the southern Korea Peninsula (the domain framed by dashed line in Figure 4(b). Figure 17 compares the analysis increments of mixing ratio near the surface between Exp-B_CV5 and Exp-B_CV5e and initial difference inmixing ratio between them. It can be seen clearly that the assimilation of the observations over the southern Korea Peninsula results in an obvious decrease of mixing ratio in the MABL in the initial condition. Comparative results of the forecasted sea fog area between Exp-B_CV5/6 and Exp-B_CV5/6e are shown in Figure 18. Compared to Exp-B_CV5, the fog area simulated by Exp-B_CV5e is significantly improved by 29.0-45.4%. However, the difference between results of Exp-B_CV6e and Exp-B_CV6 is no more than ∼5%. e above results indicate that the assimilation effect of CV6 is more reliable than CV5.  e results of the above case studies show that CV6 outstrips CV5 in the both cases. It seems that the performance of CV6 is much better in Case-A than in Case-B (see Figures 7 and 12). In Case-A, the sea fog patch is small and very near the western coast of the Korean Peninsula. e assimilation of temperature observations over the Korean Peninsula can affect the area where the sea fog occurred (Figure 11(g)), which leads to relatively great improvement in Case-A. In contrast, the sea fog patch in Case-B is much larger than that in Case-A, and most of the sea fog patch in Case-B is far away from the coast; the assimilation information of observations over land is difficult to spread over the area where sea fog occurred. is results in a relatively poor performance in Case-B.
Due to the negative correlation between temperature and moisture in CV6, assimilation of negative observational innovations of temperature can definitely get positive moisture increments, which leads to improvement of sea fog simulation.
However, degree of improvement (i.e., development of sea fog area) is primarily dominated by the localization scales in BE, which limits the spatial spread of observational information. e case studies indicate that the impact of CV6 depends on weather and real-time condition. Actually, this is directly due to the individual localization scales in BEs for Case-A and Case-B. e localization scales are calculated when BE is generated by the NMC method. As mentioned in Section 2.2, forecast differences during 15 days are used to produce BE in the NMC method. us, BE comes from the temporal average of forecast differences, which means that BE does not vary with time and sometimes the localization scales are not suitable for a specific real-time condition, especially that the weather system changes rapidly. If ensemble forecast members are available from an ensemble forecast system, ensemble-based perturbations can be used instead of using forecast differences to generate CV6 in the NMC method. Perhaps this kind of CV6 becomes more effective, because some flow-dependent information could be introduced into the BE.

Quantitative Evaluation.
e evaluation method of the simulated sea fog area in case study is a so-called eyeball method, which visually compares forecasted fog and satellite observed fog. is is a qualitative and subjective method, which cannot provide a quantitative statistical result for multiple cases. According to Zhou and Du [56], a more objective and comprehensive method is to remap both the observed and forecasted fog areas onto the same grids, in which point-topoint comparisons can be conducted. Fog area prediction can be regarded as a binary event (yes or no; 1 or 0). Statistical scores, including the probability of detection (POD), false alarm ratio (FAR), bias, and equitable threat score (ETS), are used for evaluation. ey are defined as follows:  14 Advances in Meteorology bias � Here, H, F, and O refer to the numbers of correctly forecast points (hits), forecast points, and observed points; R � F(O/N) is a random hit penalty, and N is the total number of grid points. e observed fog areas are retrieved from the MTSAT using the method designed by [22], and D2 is taken as the verification domain with a mesh resolution of 0.1°. Grids covered by high-altitude clouds are excluded over both the land and the sea areas. Verification was carried out for all model outputs at 1-hour intervals. Quantitative evaluation was conducted using the statistical scores defined above for each case, and the scores were then temporally averaged. e evaluation results are presented in Table 3 Figure 17: (a) Analysis increments of mixing ratio near the surface for Exp-B_CV5 (a) and Exp-B_CV5e (b) and differences in the initial mixing ratio between them (Exp-B_CV5e minus Exp-B_CV5) (c).  Figure 15: Initial differences of temperature (a), mixing ratio (b), and relative humidity (c) near the surface between Exp-B_CV6 and Exp-B_CV5 (the former minus the latter).  Because ETS is a comprehensive scoring indicator, the result of ETS is analyzed first. It is clearly seen from Table 3 that the average ETS of Group-CV6 is much higher than that of Croup-CV5 with an improvement of 21%, and only one case (Case-5 in Table 3) gets a worse ETS among the ten cases. It demonstrates that implementing CV-6 BE in the WRF-3DVAR can obviously improve sea fog modeling. Further case-by-case analyses of POD, FAR, bias, and ETS indicate that the improvement of ETS is mainly attributed to the increase of POD. Compared to that for Group-CV5, POD for Group-CV6 improved by 27.6% on average. Meanwhile, the average decreases in FAR and bias were −5.5% and −22.9%, respectively. Similar to the cases studied in the present study, CV6 results in more water vapor in the initial condition compared to CV5 for almost every case here and hence promotes sea fog development and larger fog area (see POD improvements in Table 3). Although the results of FAR and bias show that the sea fog areas are a bit overforecasted in Group-CV6, this is acceptable because the average deterioration of FAR is only 5.5%.

Summary
Data assimilation is absolutely necessary for numerical modeling of sea fog, because modeling result is highly sensitive to initial condition, especially to initial moisture status. e WRF-ARW model and WRF-3DVAR module have been widely employed for sea fog modeling, and the WRF-3DVAR has already been proven to be able to effectively improve the initial conditions for sea fog simulations  Figure 18: Boundaries of the forecasted sea fog area for Exp-B_CV5 (red lines) and Exp-B_CV5e (blue lines) (upper row) and for Exp-B_CV6 (red lines) and Exp-B_CV6e (blue lines) (lower row). e number in each panel is the enlargement rate (%) of sea fog area produced by Exp-B_CV#e compared to that by Exp-B_CV#.
over the Yellow Sea. However, the default domain-dependent BE of the WRF-3DVAR is CV5, which does not consider the relationship between moisture and other control variables. As a result, the improvement of moisture in the initial condition has to be only dependent on the direct assimilation of moisture information. is approach causes some disadvantageous issues; for example, other control variables like temperature can hardly contribute to the improvement of moisture status. By theoretically comparing CV5 with CV6 (multivariate BE), it is found that CV6 includes the correlation between moisture and temperature, which is absent in CV5. In order to explore the impact of multivariate BE on the 3DVAR assimilation effect over the Yellow Sea fog modeling, two sea fog cases that differ greatly were selected for case study. e advantage of CV6 compared to CV5 has been primarily explained by single-observation experiments and further revealed in detail based on analysis of the modeling results from the cases study. In addition, the impact of CV6 application in the WRF-3DVAR has been evaluated by a series of experiments of extra ten sea fog events over the Yellow Sea. e major conclusions are as follows: (1) e performance of the WRF-3DVAR assimilation with CV6 is obviously better than that with CV5. e assimilation with CV6 can significantly improve the forecasted sea fog area, which is clearly demonstrated in the study of the selected cases. Particularly, for the case of spreading narrowly along the coast, the sea fog evolution is successfully reproduced when the assimilation with CV6 is implemented. In contrast, the model completely fails to capture this fog event when using the assimilation with CV5.
(2) ere exists a dominant negative correlation between temperature and moisture in CV6 near the sea surface, which makes it possible to improve initial moisture condition in the marine atmospheric boundary layer by assimilating temperature observations near the sea surface. However, the primary way to improve moisture in the WRF-3DVAR with CV5 is to assimilate moisture observations.
(3) Experimental results of the 10 additional sea fog cases clearly indicate that CV6 is more suitable than CV5 in the WRF-3DVAR assimilation for sea fog modeling. Compared to that with CV5, the average POD and ETS of forecasted sea fog area using the WRF-3DVAR with CV6 can be improved by 27.6% and 21.0%, respectively, and positive impacts are found in 9 out of the 10 cases in the present study.
e present study has demonstrated the importance of improving moisture condition in the initial condition of sea fog modeling and the necessity of using CV6 in the WRF-3DVAR to achieve this improvement. is is of reference value for the real-time operational forecast of sea fog when using the WRF-3DVAR. Despite encouraging results that have been achieved in this study, further in-depth studies are still needed. CV6 is essentially a type of climatological background error covariance, which usually varies on monthly basis and may also have a feature of rapid variation in the low-level atmosphere due to transient weather systems. us, when ensemble forecast members are available, it is better to use ensemble-based perturbations instead of using forecast differences in the NMC method. is is because some flow-dependent information might be introduced into the BE. In addition, the localizations of horizontal and vertical correlation in CV6 represent an issue that needs to be addressed.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.