The Regularized WSM6 Microphysical Scheme and Its Validation in WRF 4D-Var

Sen YANG, Deqin LI*, Liqiang CHEN, Zhiquan LIU, Xiang-Yu HUANG, and Xiao PAN

1Institute of Atmospheric Environment, China Meteorological Administration, Shenyang 100166, China

2National Center for Atmospheric Research, Boulder, CO 80301, USA

3Institute of Urban Meteorology, China Meteorological Administration, Beijing 100089, China

ABSTRACT A cold cloud assimilation scheme was developed that fully considers the water substances, i.e., water vapor, cloud water, rain, ice, snow, and graupel, based on the single-moment WSM6 microphysical scheme and four-dimensional variational (4D-Var) data assimilation in the Weather Research and Forecasting data assimilation (WRFDA) system. The verification of the regularized WSM6 and its tangent linearity model (TLM) and adjoint mode model (ADM) was proven successful. Two groups of single observation and real sounding data assimilation experiments were set up to further verify the correctness of the assimilation scheme. The results showed that the consideration of ice, snow, and graupel in the assimilation system of the 4D-Var, as opposed to their omission in the warm rain Kessler scheme, allowed the water substances to be reasonably updated, further improving the forecast. Before it can be further applied in the assimilation of observational data, radar reflectivities, and satellite radiances, the cold cloud assimilation scheme needs additional verification, including using conventional ground and sounding observations in the 4D-Var assimilation system.

Key words: 4D-Var, data assimilation, linearization, numerical weather prediction, WSM6

1. Introduction

Since four-dimensional variational (4D-Var) data assimilation was introduced for use in atmospheric research and operational centers, it has become one of the most advanced and widely used data assimilation techniques (Lewis and Derber, 1985; Le Dimet and Talagrand, 1986; Zou et al., 1992;Courtier et al., 1994; Huang et al., 2009; Zhang et al., 2014;Liu et al., 2020). Compared to the three-dimensional variational (3D-Var) scheme, the ensemble Kalman filter (EnKF),the hybrid ensemble variational, and other data assimilation techniques (Rakesh and Kutty, 2021), the 4D-Var scheme has three main useful features (Rabier et al., 2000; Huang et al., 2009; Sun et al., 2021):

(1) Implicit definition of flow-dependent forecast error covariance.

(2) The NWP model and physics process, which builds a strong constraint to maintain a dynamic balance in the minimization.

(3) The ability to use observations at the time of their measurement within the assimilation time windows.

4D-Var data assimilation was first introduced into meteorology in simple models (Lewis and Derber, 1985; Le Dimet and Talagrand, 1986), and then complex numerical weather predictions (NWPs) were conducted (Navon et al.,1992; Zou et al., 1995; Rabier et al., 1998). The European Centre for Medium-Range Weather Forecasts (ECMWF)was the first center to operationally implement the 4D-Var(Rabier et al., 2000). Following the success of the ECMWF,other 4D-Var systems based on operational NWP models have been developed and used in global NWP centers, including in France (Janisková et al., 1999; Gauthier and Thépaut,2001), the United Kingdom (Rawlins et al., 2007), Canada(Gauthier et al., 2007), Japan (Kadowaki, 2005), the United States Navy (Xu et al., 2005; Rosmond and Xu, 2006), and China (Liu et al., 2018; Zhang et al., 2019). 4D-Var has also been developed and used in regional models by some NWP centers (Honda et al., 2005) and research institutes (Sun and Crook, 1997; Huang et al., 2002; Huang et al., 2009; Liu et al., 2020; Rakesh and Kutty, 2021).

Most of the operational 4D-Var implementations have adopted an incremental formulation (Courtier et al., 1994;Lorenc, 2003; Huang et al., 2009), which involves running a high-resolution nonlinear model in the outer loop to calculate the departure of the model trajectory from observations and low-resolution tangent linear and adjoint mode models(TLM and ADM, respectively) in the inner loop. This is also called a multi-resolution incremental 4D-Var (MRI-4DVar) (Zhang et al., 2014; Liu et al., 2020). The Weather Research and Forecasting (WRF) model is one of the most advanced NWP models, and it is freely available and widely used to develop operational regional numerical weather forecast systems around the world. The 4D-Var in WRF data assimilation (WRFDA) was developed by extending the 3DVar to the time dimension based on its basic framework(Huang et al., 2009). To decrease the computational cost of 4D-Var minimization, Zhang et al. (2014) and Liu et al.(2020) developed an MRI-4DVar version that uses multiple horizontal resolutions for the inner loop and a three-stage procedure (observer, minimizer, and regridder) in the outer loop.

In 4D-Var, the NWP model is a strong constraint; however, it usually employs simple or simplified physics processes, which are especially important in high-resolution models (Mahfouf and Rabier, 2000; Rabier et al., 2000). An accurate representation of microphysical processes is critical for improving cloud and precipitation forecasts in high-resolution models. To better understand the predictability of clouds and precipitation, microphysical parameterization schemes (MPSs) are used to represent the realism of cloud microphysical processes. Most MPSs are derived from the early bulk microphysics schemes, e.g., the Lin scheme and the Rutledge-Hobbs scheme (Lin et al., 1983; Rutledge and Hobbs, 1984). In MPSs, the hydrometeor size spectrum is assumed to follow an exponential (Kessler, 1969) or gamma distribution (Walko et al., 1995), which is related to the particle size distribution, parameterized mass-diameter, and fall speed-diameter relationships. Simple schemes usually use a small number of cloud species (e.g., cloud droplets, rain,and ice). More complex schemes involve additional finer cloud species (graupel, hail, ice, and snow) and use more prognostics for each cloud species (Milbrandt and Yau,2005; Morrison et al., 2005; Hong et al., 2006; Thompson et al., 2008; Furtado et al., 2018).

In variational assimilation, the TLM and ADM of the moist physical processes link the changes in the 4D-Var control variables with the changes in cloud and precipitation fields (Janisková and Lopez, 2013). Due to the highly nonlinear characteristics and complexity of MPSs, developing a linearized version is challenging. Although the current NWP model uses more advanced MPSs in its forecast model, its data assimilation system usually employs relatively simple MPSs, which typically lack ice-phase hydrometeors, as do its TLM and ADM. In WRFDA, a warm rain Kessler scheme (Dudhia, 1989) is developed and used to constrain the relationships among rainwater, cloud water, moisture,and temperature. The Kessler scheme includes the condensation process of water vapor into the cloud, accretion of the cloud by rain, automatic conversion of the cloud to rain, and the evaporation of rain to water vapor (Xiao et al., 2007).However, the Kessler scheme only builds a constraint among rainwater, cloud water, moisture, and temperature but lacks ice-phase microphysical processes (Wang et al.,2013). Because moisture processes usually have discontinuities and strong nonlinearity, it is difficult to develop a TLM and ADM for most ice-inclusive MPSs, which may cause a variational assimilation to fail (Geer et al., 2017). Although ignoring ice-phase water substances in the current version of WRFDA does not cause serious forecast deviations for assimilated humidity variables (Wang et al., 2013), it is still desirable to develop a more complex mixed-phase MPS to assimilate humidity and to work towards including a cloud and precipitation assimilation.

Therefore, to fully consider ice-phase hydrometeors in 4D-Var, this study developed a regularized model of WSM6 and its TLM and ADM for 4D-Var in WRFDA. The remainder of this article is organized as follows. Section 2 introduces the 4D-Var in WRFDA and the regularized WSM6 MPS, as well as the development and verification of the TLM and ADM. Sections 3 and 4 summarize the single observation and real-sounding data assimilation experiments conducted to verify the correctness of the assimilation scheme, and the conclusions are presented in section 5.

2. Construction of a cold cloud assimilation scheme based on WSM6 in WRF 4D-Var

2.1. 4D-Var data assimilation

The 4D-Var in WRFDA is the temporal extension of 3D-Var and uses the adjoint model’s backward integration to calculate the gradient of the cost function over an assimilation window (Huang et al., 2009; Zhang et al., 2014; Liu et al., 2020). With the model as a strong constraint during minimization, the cost function in 4D-Var is expressed as:

where N is the number of discrete subwindows within a single assimilation time window, i as an individual subwindow; x and xbare the analysis and background fields of the model variables that are valid at the beginning of the assimilation window, respectively; yiis the observation vector at time i;B and R represent the background and observation error covariance matrixes, respectively; H is the observation operator; Midenotes the nonlinear model integration from time 0 to time i; and the superscripts “—1” and “T” denote the inverse of a matrix and the transpose of a column vector,respectively.

Using xgas a first guess, let δx = x — xgand δxg= xb—xg; then, the incremental form of 4D-Var is expressed as follows (Huang et al., 2009; Wang et al., 2013; Liu et al.,2020):

where di=yi-Hi[Mi(xg)] is the first guess of the departure of the trajectory from observations, which is calculated using the nonlinear WRF model and observation operator.H and M are the tangent linear operators of H and M, respectively.

Because B is a huge matrix for the original model state variables, the control variable transform U is introduced to reduce the condition number and to accelerate the minimization algorithm. Letting B = UUT, δx = Uv, and δxg= Uvg;the control variable form of Eq. (3) becomes:

and the gradient of the cost function, with respect to v is:

where HTand MTare the adjoints of H and M.

After the minimization procedure converges at the specified number of iterations and convergence conditions, the analysis of model state variables can be updated. Further details of the calculation of the control variable transforms,U and UT, can be found in Huang et al. (2009) and Liu et al.(2020).

2.2. WSM6 microphysics scheme

This study developed a linearized version of the WSM6 and its TLM and ADM in WRF 4D-Var. There are two main factors to consider when choosing WSM6. One is that it is easier to linearize a single-moment MPS than a more complex scheme, and it is also easier to develop its TLM and ADM. Another is that WSM6 is one of the most widely used and well-behaved schemes in the northeast Asia area in the meso-scale NWP system (Teng et al., 2020).

The WSM6 scheme is a bulk microphysics scheme that includes ice-phase water substances (cloud ice, snow, and graupel). Most of the processes in WSM6 are consistent with the microphysical scheme developed by Rutledge and Hobbs (1983) and Dudhia (1989), although its cloud-rain auto-conversion adopts the scheme developed by Tripoli and Cotton (1980). The WSM6 scheme was developed by adding additional processes related to graupel species into the WSM5 scheme (Hong and Lim, 2006). With an increase in the number of different water substances, the amount of rainfall in WSM6 increases, and the local maximum becomes stronger than in the WSM5 in high-resolution forecast experiments (Hong and Lim, 2006). In WSM6, six classes of prognostic water substances are considered: the mixing ratios of water vapor (Qvapor), cloud liquid water(Qcloud), cloud ice (Qice), snow (Qsnow), rain (Qrain), and graupel (Qgraupel) (Hong and Lim, 2006). There are 34 production terms for six water substance transformation processes in WSM6, and the detailed computational procedures are given by Hong and Lim (2006). Hong et al. (2009) also explained the differences between the LIN and WSM6 schemes in detail. This study developed a linearized version of the WSM6 and incorporated it into the WRF TLM and ADM (WRFPLUS) code.

2.3. Regularized WSM6 for the construction of a cold cloud assimilation scheme

Processes such as condensation, convection, and vertical diffusion, in most MPSs are highly nonlinear and often discontinuous (and therefore not differentiable), which may cause convergence issues during 4D-Var minimization when directly linearizing them. Discontinuous processes,such as those related to the presence of threshold or on/off processes (Zou, 1997), need to be regularized before developing their linearized version (Janisková et al., 1999). Regularization is used to develop a smooth and regular physical package (avoiding problems of convergence), as long as the parameterizations can represent the order of magnitude in intensity and the general feedback of physical phenomena present in the atmosphere (Janisková et al., 1999). The purpose of regularization is to make the physical parameterizations more continuous and to give the 4D-Var a better physical process.The derivation of a tangent-linear physics parameterization by making a simplified, regularized version of the nonlinear model can be used to remove or smooth over the worst of the nonlinearities (Janisková and Lopez, 2013). The tangentlinear equivalent is then created, and careful testing and further regularization may be necessary if the initial models show excessively strong instabilities (Geer et al., 2017). Various technical modifications, such as recording the on/off switches, spline fitting, and vertical smoothing, are used to replace the original parameterizations of the moisture processes in NWPs (Zou et al., 1993; Zupanski, 1993; Zou,1997; Janisková et al., 1999).In WSM6, the distinction between gaseous, solid, and liquid precipitation involves discontinuities related to the temperature and water substances (Hong et al., 2006). For example,in the homogeneous freezing of cloud water, when the temperature is lower than —40°C, a discontinuity will occur as cloud water is transformed into cloud ice. In the auto-conversion process from cloud to rain, a discontinuity will occur when the mixing ratio of cloud water is larger than 5.0265 ×10—4g kg—1.In this approach, an S-shaped curve logistic function is used to smooth the discontinuities in WSM6. The equation is as follows:

where c is the value of the sigmoid midpoint and k is the logistic growth rate or steepness of the curve. The larger the value of k used by the function, the closer it is to the IF statement (e.g., Heaviside step function). A smaller k corresponds to a smoother function, as does its derivative. Because the hard-coded minimum threshold of the water substances in WSM6 is 10—15, k = 1017can guarantee f (—10—15, 0) = 0 and f (10—15, 0) = 1 in Eq. (5). For the temperature, the hardcoded minimum threshold of temperature in WSM6 is 10—9,k = 1011can guarantee f (—10—9, 0) = 0 and f (10—9, 0) = 1 in Eq. (5). Therefore, k = 1011was used in the smoothing of the temperature-related IF statement and k = 1017was used in the smoothing of the water substance-related process.

In developing the TLM and ADM of WSM6, it is difficult to find the errors in the consistency check because the water substances and temperature are updated after all conversion processes are completed in the original version of WSM6.Therefore, the water substances and temperature are updated upon completion of each physical process in the regularized WSM6.

The magnitude and precision of different variables differ among MPSs. It is easy to ignore a small variation in magnitude, but this could cause the rounding-off of digits when calculating the sum of two variables. For example, in the calculation of the production rate for the deposition—sublimation rate of ice (PIDEP), there is a need to update the mixing ratio of water vapor, mixing ratio of ice crystals, and temperature because the summation will result in a rounding error,leading to further errors in the calculation of the TLM and ADM. In addition, using intermediate variables will also cause error accumulation. To eliminate error accumulation,all of the production terms of the physical processes in WSM6 were rewritten, and the intermediate variables were deleted, particularly the densities of the different water substances. After all of the above regularization steps were completed, the regularized version of WSM6 was named RWSM6. A detailed comparative analysis of the water substances simulated from WSM6 and RWSM6 is discussed in section 4.2.

2.4. Construction of the TLM and ADM of RWSM6

The TLM and ADM of the existing programming model require a manual derivation of the model equation and computer differential processing; however, this becomes almost impossible as the model becomes more complex. Automatic differentiation is a collection of techniques used to obtain the analysis derivatives of differentiable functions, which is provided in the form of a computer program(Hascoet and Pascual, 2013). In recent years, significant efforts have been made to develop a 4D-Var system based on the WRF model and to couple the variational data assimilation algorithm (WRFVAR) with the WRFPLUS. Xiao et al. (2008) simplified the parametrization routines of the WRF model and then used the Transformation of Algorithm (TAF)automatic differentiation tool (Giering and Kaminski, 2003)to generate the TLM and ADM of the WRF model. Manual modifications and tests were conducted to verify its correctness. TAPENADE is one of the most effective automatic differentiation tools and can be used to automatically develop adjoint models given the source of an original program (Hascoet and Pascual, 2013). Zhang et al. (2013) used TAPANADE to develop a new TLM and ADM version of the WRF model with manual intervention and found that it showed efficient parallel performance and scalability.In the present study, TAPENADE was also used to generate the TLM and ADM of RWSM6, and the interface of the microphysical process in WRFPLUS was added. The correctness of the TLM and ADM of RWSM6 was tested within WRFPLUS, and the detailed results are discussed in section 4.3. Figure 1 presents a conceptual diagram of the regularization, development, and verification of the TLM and ADM of WSM6.

2.5. Cloud control variables and background error statistics

In WRFDA, the total water mixing ratio is used as a control variable to assimilate most of the observations. Xiao et al. (2007) used a warm rain scheme and the total water mixing ratio to assimilate radar reflectivity. The results showed that the short-range quantitative precipitation forecast (QPF)skills improved after the assimilation of Doppler radar data.When humidity-related observations were assimilated, the forward warm rain process and its backward adjoint distributed this information to the increases in other variables. This also means that multivariate correlations were obtained in 4DVar through the NWP model (Xiao et al., 2002, 2007).Unlike the warm-rain scheme, WSM6 includes the mixing ratios of water vapor, cloud water, rain, cloud ice, snow,and graupel, which can correspond to the water substance control variables in the 4D-Var of WRFDA. The types of water substance variables contained in WSM6 were also more suitable for convective synoptic systems than the Kessler scheme. The moisture control variable is the pseudo relative humidity (defined as Q/Qb,s, where Q and Qb,sare the specific humidity and saturated specific humidity from the background field), cloud water, rainwater, cloud ice, snow, and graupel. Both horizontal and vertical correlations are considered. The background error statistics of the cloud variables need to be included in the be.dat file. The generalized software package (GEN_BE), developed by Descombes et al. (2015),was used to generate the background error statistics for control variables. The background error statistics were computed using the National Meteorological Center method (Parrish and Derber, 1992), which uses pairs of the differences between the 12- and 24-h forecasts. A total of 31 daily forecasts in July 2018 were used to generate the background error covariance. Background error covariance (BEC) implicitly evolves within the time window through a linearized model. Although static B (BEC at the beginning of the time window) is typically the same for all analysis cycles, the 4DVar sets up the multivariate correlations through its backward adjoint integration and distributes the increment information to the other variables.

3. Validation of the RWSM6 TLM and ADM

3.1. Case of Typhoon Ampil (2018)

The case of Typhoon Ampil (2018) in the western Pacific in 2018 was used to validate the correctness of the regularized RWSM6, and the results of the 4D-Var assimilation were used to forecast sensitivity experiments. Figure 2 shows the path of Typhoon Ampil (2018), and the best path of the typhoon is from the Japan Meteorological Agency(https://www.jma.go.jp/). At 2000 UTC 18 July, Ampil was generated over the northwest Pacific Ocean and then moved towards the northeast. On 19 July at 2000 UTC, Ampil turned and moved steadily north, skirting the coast of Chongming Island, Shanghai, with the intensity of a strong tropical storm (central maximum wind speed: 28 m s—1; lowest central pressure: 982 hPa). Ampil passed through Jiangsu, Shandong, and Hebei provinces, then turned northeast on 24 July at 0600. It entered Liaoning on 24 July at 1700 UTC and moved away at 2200 UTC. The maximum precipitation in western Liaoning was 146 mm over 24 h. After Ampil weakened and merged into a westerly wind trough, precipitation occurred again in Liaoning. The observed maximum precipita-tion at the auto weather station in southwestern Liaoning was 182.5 mm between 0500 and 2000 UTC on 25 July,and the rainfall rates at the Huludao and Dalian sites reached 86 and 105 mm h—1, respectively. The main feature of this case is that the precipitation in the second stage affected by Ampil was more serious than that in the first stage. In this study, we mainly focused on the second stage of precipitation in Liaoning Province.

3.2. Verification of the regularized WSM6

A comparison experiment was conducted to verify the correctness of the linearized RWSM6 scheme. The WRFv4.2 was run with the original WSM6 and RWSM6,respectively, from 0000 UTC 25 July 2018 to 1200 UTC 25 July 2018, with the 6-hour forecast at 1800 UTC 24 July 2018 as the initial and boundary conditions of the National Centers for Environmental Prediction Global Forecast System data. The horizontal grid spacing was 3 km in a single domain, with 401 × 351 grid points. The time step was 18 s.The number of vertical levels was 50, and the model top was at 10 hPa. The other physical parameterizations selected for use were the Rapid Radiative Transfer Model(RRTM) longwave radiation scheme (Mlawer et al., 1997),Dudhia shortwave radiation scheme (Dudhia, 1989), Yonsei University (YSU) planetary boundary layer scheme (Hong et al., 2006), and Noah land surface model (Chen and Dudhia, 2001). The cumulus parameterization was turned off in this 3 km setting.

The water substances produced by WSM6 and RWSM6 were water vapor, cloud water, rain, cloud ice,snow, and graupel. The magnitude of water vapor was larger than that of other water substances, and there were also differences in their spatial distribution. The spatial distributions of water vapor predicted by WSM6 and RWSM6 were similar (not shown). Figure 3 shows the mixing ratios of cloud water (Qcloud) and rain (Qrain) at 850 hPa and ice(Qice), snow (Qsnow), and graupel (Qgraupel) at 300 hPa were obtained with WSM6 and RWSM6. The figure also shows that RWSM6 could capture the spatial distribution characteristics of hydrometeors from the original WSM6. There were some differences in the simulation of water substances by WSM6 and RWSM6, especially for the extreme values of Qsnowand Qgraupelfrom RWSM6, which were slightly lower than those from WSM6. In recalling the regularization process, the strict thresholds of water substances and temperature were carefully chosen by k in the logistic function, which was not the main source of differences in the ice-phase hydrometeors simulated by WSM6 and RWSM6. One possible source of the differences was that the water substances and temperature were updated immediately after each physical process was completed. This could have changed the triggering and calling conditions for some of the physical processes.

Precipitation is sensitive to microphysical parameterization (Hong and Pan, 1998). The 12 h precipitation obtained from the dense automatic weather stations (DAWS) in China from 0000 UTC 25 July 2018 to 1200 UTC 25 July 2018 was used to evaluate the difference in the precipitation forecast between WSM6 and RWSM6 (Fig. 4). Both WSM6 and RWSM6 captured the observed spatial distribution of precipitation in Fig. 4. As shown in Fig. 3, the differences in water substances between WSM6 and RWSM6 further contributed to the precipitation, and there were some differences in the details, as shown in Fig. 4.

Fig. 1. Concept diagram of the regularization of WSM6, and the development and verification of TLM and ADM.

Fig. 2. The track of Typhoon Ampil (2018). The dotted line is the model domain, and the best track of the typhoon was obtained from the Japan Meteorological Agency(https://www.jma.go.jp/).

Fig. 3. (Continued).

Fig. 3. Comparison of the (a, b) Qcloud and the (c, d) Qrain forecasts at 850 hPa, and (e, f) Qice, (g, h) Qsnow,and (i, j) Qgraupel forecasts at 300 hPa with (a, c, e, g, i) the original WSM6 and (b, d, f, h, j) the regularized RWSM6 at 0000 UTC on 25 July 2018.

3.3. Validation of the TLM and ADM of RWSM6

According to Thépaut and Courtier (1991) and Navon et al. (1992), the TLM and ADM can be validated using the TLM and ADM test procedures. In WRFDA, WRFPLUS contains the TLM and ADM based on a simplified WRF model.WRFPLUS includes a few simplified physics packages, and it also includes tools for TLM and ADM tests following the method of Navon et al. (1992). The test results and sensitivities obtained using WRFDA are reported in Xiao et al.(2008). Based on WRFDA, the interfacing of RWSM6 to WRFPLUS was completed successfully, and the TLM and ADM test of RWSM6 also worked well.

As in Navon et al. (1992) and Zhang et al. (2014), let NLM and TLM denote the nonlinear model and tangent linearity, where x and Δx are the column vectors of model-state variables and perturbations of state variables, respectively.The correctness of the TLM can be determined as follows:

where ＜, ＞ is the inner product of the two vectors. For values of α that are small but not close to zero, the ratio Φ(α)should be close to 1. For a tangent linear check using the 4DVar in the WRFDA system, Table 1 shows the results of Φ(α) for values of α ranging from 10—1to 10—11with an integration time of 60 min, and the numerator and denominator are also shown separately. To ensure that WSM6 is called,the increments of humidity in the upper layer (200 hPa) and at the surface were added as 1 and 10 g kg—1, respectively.Due to the highly nonlinear relationship of microphysical processes, Φ(α) was not close to 1 when n = 7, 8, and 9. Similar results were found in Navon et al., (1992) and Wang et al.(2013). The effect of the length of integration time (5, 15,30, 60, and 180 min) was also tested. Although there were differences among the test results, it could be guaranteed that Φ(α) was close to 1 for n = 11. Therefore, the results from Table 1 confirmed the correctness of the TLM of RWSM6.

Table 1. Results of the WRF TLM correctness check with 60 min integration time.

As in Zhang et al. (2014) and Wang et al. (2013), the accuracy of the ADM with the RWSM6 scheme was verified against the TLM based on the following definition:

where the left side involves only the TLM, while the right side involves the ADM. It is necessary for the values of the left and right sides to be exactly the same.

The precision of the numbers on both sides of Eq. (5) is widely used for verifying the correctness of the ADM(Wang et al., 2013; Zhang et al., 2014). According to Eq.(5), an integration for 1 h gives the left-side a value of 0.42361302578090×1012, while the right-side value is 0.42361302578096×1012, representing a 13-digit accuracy on 64-bit machines, which has similar accuracy to that reported in Wang et al. (2013). This confirmed the mathematical accuracy of the WRFPLUS with RWSM6.

4. Results

4.1. Single observation experiments

A single observation experiment was conducted to illustrate the impact of the temporal evolution of the BEC matrix and the TLM and ADM integration (Huang et al., 2009;Wang and Lei, 2014, Chen et al., 2021; Sun et al., 2021).The model played a role in propagating the observed information in the 4D-Var analysis, and a single observation test clearly showed the flow-dependent covariance structure(Huang et al., 2009; Wang et al., 2013).

To illustrate the correctness of the TLM and ADM of RWSM6, the results obtained using the linearized Kessler scheme were also introduced as a reference. Two sets of single observation experiments (Sing_Kessler and Sing_WSM6)were conducted to compare the difference when using the warm rain Kessler scheme as opposed to the cold cloud WSM6 scheme. The model setting of the two single observation experiments was the same as that described in section 4.2, and the MPS in the WRF model was WSM6. The forecast field after 6 h of model integration was used as the initial field. A single specific humidity observation at 0000 UTC on 25 July was placed upstream of the weather system at 119.81°E, 40.73°N, 500 hPa (z = 18) and was valid at the end of the 3-h assimilation window. The value of the specific humidity observation increment was set at 1 g kg—1, and the standard deviation of the observation error was set at 1 g kg—1.

Because the Qvaporwas larger than that of other water substances, Fig. 5 shows the increment of the initial field of the Qvaporfrom two sets of single-point experiments, as well as the forecasts at 1, 2, and 3 h within the 3-h assimilation window. Due to the effect of the weather system, after adding a single point of humidity observation, the increment from the two groups of experiments continued to move toward the northeast with the integration of the WRF model. The water vapor in the two single-observation experiments existed at 500 hPa, where the corresponding temperature was well below freezing. Despite identical initializations of the humidity field in both experiments, the difference in the predicted humidity by the models gradually increased over the 3-h assimilation window. Due to the lack of ice-phase physical processes in the Kessler scheme, Qvapordid not change much over the entire 4D-Var time window in the Sing_Kessler experiment, and the Qvaporfrom the Sing_Kessler experiment was stronger than that of the Sing_WSM6 experiment in all time windows. This can be understood due to a more reasonable distribution of water substances when using WSM6 rather than the Kessler scheme.

Fig. 4. Accumulated precipitation (mm) from 00 UTC on July 25 to 12 UTC on July 25 for (a) the observation and (b) the forecasts from the experiments with the original WSM6 and (c) the regularized RWSM6.

Fig. 5. The increment in Qvapor (g kg—1) at 500 hPa at (a, e) 0 h, (b, f) 1 h, (c, g) 2 h, and (d, h) 3 h after assimilating a single specific humidity observation at the end of the 3-h assimilation window in the single observation experiments of (a—d)Sing_Kessler and (e—h) Sing_WSM6. The red lines from A (39.00°N, 117.51°E) to B (40.96°N, 120.09°E) in panels (a) and (e)were used for the vertical cross-sections shown in Figs. 4 and 5. The blue + indicates the location of observation point C.

In the Sing_Kessler experiment, the increments of only three water substances (Qvapor, Qcloud, and Qrain) were updated with the assimilation of a single humidity observation(not shown). Figures 6a—c show cross-sections of the analysis increments of Qvapor, Qcloud, and Qrainin Sing_Kessler, and Figs. 6d—i show cross-sections of the analysis increments of Qvapor, Qcloud, Qrain, Qsnow, Qice, and Qgraupelin Sing_WSM6.In comparing the increments of the water substances, the Qvapor, Qcloud, and Qrainobtained by the two sets of experiments were similar, but there was more Qrainpresent in Sing_WSM6 than in Sing_Kessler. Because the Qvaporwas much larger than the other water substances, the Qcloud, Qrain,Qsnow, Qice, and Qgraupel, and their icons are enlarged by 1000 times in Fig. 6. There were larger amounts of Qgraupelthan Qsnowand Qicebecause the microphysical processes responsible for graupel accounted for most of the snow in WSM6 based on WSM5 (Fig. 6) (Hong and Lim, 2006).Based on the vertical sections of Qsnow, Qice, and Qgraupel, the heights of the three types of water substances were significantly higher than those of the other three water substances,between 600 and 200 hPa.

Figure 7 shows cross-sections of the initial Qvaporand the forecasts at 1, 2, and 3 h within the 3-h assimilation window along the two endpoints of A (39.59°N, 118.03°E) to B(40.73°N, 119.81°E) in Figs. 5a and 5e. Even with the same setting in the forward WRF model, the vertical structure of Qvaporvaried greatly between the two experiments. As shown in Fig. 5, with the warm rain Kessler scheme, the intensity and height of the water vapor were nearly unchanged in the 3-h time window because of the lack of a process to convert water vapor to other water substances. The positive increments of Qvaporabove 700 hPa in Sing_Kessler were stronger than those in Sing_WSM6. The results also showed that the increment of Qvaportriggered from the convection system gradually replaced the single-point information with the forward integration of the WRF model, and the differences in the incremental Qvaporgradually increased as the two single-observation experiments were run.

4.2. Real sounding data assimilation experiment

Real-sounding data were also used to verify the correctness of the TLM and ADM of WSM6 as shown in this section. Two real-sounding data assimilation experiments using Kessler (Real-Kessler) and WSM6 (Real-WSM6) were conducted for precipitation during the second stage of precipitation from Typhoon Ampil (2018). Each experiment performed 12-h forecasts starting from 0000 UTC to 1200 UTC 25 July.

Because sounding observation data were available only twice a day at 0000 and 1200 UTC, and the two experiments were conducted mainly to test the correctness of the TLM and ADM of WSM6, the time window of the 4D-Var was set to 12 min. The number of outer loops in the 4D-Var was set to nine. There were 18 sounding stations (Fig. 8) in the model domain used in the real-sounding data assimilation experiments. The experiments adopted the same physics options as those used in the previous single observation experiment.

WSM6 is more complex than Kessler; therefore, more computing resources (82%) were required for the assimilation analysis in 4D-Var in this case. As shown in Fig. 9, the values of the cost functions for the two experiments were consistent. The cost functions converged well in Real WSM6 with four iterations, and they did not increase compared with Real_Kessler. This was partly due to the small amount of sounding observation data used in this test. The cost function was significantly reduced after the first three iterations. In addition, a jump in the value of the cost function occurred at each new iteration, which was explained by the start of a new outer loop (Wang et al., 2013). Four outer loops were sufficient to produce a reasonable value of the cost function,and the number of first and second internal iterations was high.

Fig. 6. The cross-sections of the increments of (a) Qvapor, (b) Qcloud, and (c) Qrain for the Sing_Kessler experiment at 0 h along A to B in Figs. 3a and 3e at 0 h and the increments of (d) Qvapor, (e) Qcloud, (f) Qrain, (g) Qice, (h) Qsnow, and (i) Qgraupel in the Sing_WSM6 experiment.

Fig. 7. The cross-sections of the increment of Qvapor (g kg—1) at (a, e) 0 h, (b, f) 1 h, (c, d) 2 h, and (d, h) 3 h after assimilating a single specific humidity observation that was valid at the end of the 3-h assimilation window along A to B in Figs. 3a and 3e in the single observation experiments of (a—d) Sing_Kessler and (e—h) Sing_WSM6.

Fig. 8. The locations of the 18 sounding stations.

Fig. 9. Cost function with respect to the number of total iterations accumulated by the nine outer loops for the Real_Kessler and Real_WSM6 experiments.

The water substances in the initial field of the model were adjusted after assimilating the sounding humidity observation data. Figure 10 shows the average profile of the initial Qcloud, Qrain, Qice, Qsnow, and Qgraupelof the two experiments and again after a 6-h forecast in the model domain. Because the magnitude of Qvaporwas larger than that of the other hydrometeors, Figs. 10c and 10f also show the difference in the vertical profiles of Qvaporfrom the Real-WSM6S and Real-Kessler experiments. In the warm-rain scheme, at the initial time after assimilating the sounding data, only three water substances were updated, with small increments in the Real_Kessler experiment. In the Real_WSM6 experiment,six water substances were updated. The differences in the initial ice-phase water substance of the two experiments were maintained for almost 1 h before being gradually reduced by the model integration. In this way, the vertical distribution of the predicted water substances after 6 h tended to be consistent (Figs. 8d—e). This was because WSM6 was used in both sets of assimilation experiments, and the subsequent integration of the model resulted in an increase in the content of water substance, which remained consistent overall. The difference in water vapor content (corresponding to Qvapor) was mainly concentrated below 700 hPa. From a comparison of Qvaporin the Real_Kessler and Real_WSM6 experiments,with assimilated sounding data, the increment of Qvaporwas mainly above 500 hPa in the Real_WSM6 experiment, and the increment of Qvaporwas mainly below 500 hPa in the Real_Kessler experiment under the initial model conditions.As with the integration of the WRF model, this difference was mainly apparent below 600 hPa. These results indicated that the spatial distribution of the initial water substance(more types and spatial distributions) could be more reasonably updated with a more complex MPS. The ice-phase hydrometeor in the initial field was redistributed with model integration of 1 h because the distribution and conversion of ice-phase water substances were strongly influenced by Qvapor, and the magnitude of Qvaporwas larger than that of the other water substances.

4.3. Precipitation forecasts in the real sounding data assimilation experiment

To evaluate the performance of the WSM6 and Kessler schemes in 4D-Var, the impact of data assimilation using two MPSs on the precipitation forecast was assessed.Figure 11 shows the 3-h accumulated precipitation from observations and the corresponding forecast rainfall from the two real-sounding experiments. For the 3-h accumulated precipitation, both the Real_Kessler and Real_WSM6 experiments did not induce the precipitation noted in the observation(Fig. 3a). This may be because cloud and precipitation observations cannot be assimilated when cloud or precipitation is missing in the first guess due to the zero-gradient problem(Geer et al., 2017). After 3 h of WRF integration, the precipitation forecasts of the two experiments were relatively close to the observations, and the predicted precipitation from the Real_Kessler experiment was slightly greater than that in the observations. As the model continued to integrate, two centers of maximum precipitation developed, and the Real_Kessler experiment began to overestimate the precipitation on the left side of the precipitation centers. Considering the results in Fig. 10, the increase in water vapor in the Kessler scheme was mainly below 500 hPa, and precipitation was more likely to occur in the upstream area of the synoptic system. The precipitation from the Real_WSM6 experiment matched the observations better than the Real_Kessler experiment.

Fig. 10. The vertical profile of the average hydrometeors of Qcloud, Qrain, Qice, Qsnow, and Qgraupel in the model domain at the (a, b)initial and (d, e) 6-h forecasts for the (a, d) Real_Kessler experiment and (d, e) Real_WSM6 experiment and the difference in the Qvapor between the Real_WSM6 and Real_Kessler experiments at the (c) initial time and (f) 6-h forecasts.

Figure 12 shows the observed 12-h accumulated precipitation from 0000 UTC to 1200 UTC on 25 July 2018, compared to the forecast precipitation from the Real_Kessler and Real_WSM6 experiments. The figure indicates that the precipitation center was mainly distributed over western Liaoning and northern Hebei. The Real_Kessler experiment forecasted precipitation bands with a stronger intensity to the north of Hebei, but these bands were weaker to the west of Liaoning. The precipitation from the Real_WSM6 experiment was more consistent with the observation than the precipitation from the Real_Kessler experiment.

The threat score (TS) and bias score verification statistics were calculated for the two assimilation experiments. A TS value of 1 indicates a perfect rainfall forecast by the experiments. The bias score indicates the tendency of a model to underestimate (when the bias score is less than 1) or overestimate (when the bias score is greater than 1) an event. As shown in Fig. 13, the Real_WSM6 assimilation experiment had a higher TS score than Real_Kessler, with more than 50 mm of precipitation within a 12-h period. Meanwhile,the results in Fig. 10 showed that the increase in the amount of water vapor and ice-phase water substances in the initial field in Real_WSM6 was greater than that in Real_Kessler above 500 hPa. This indicates that using a more complex MPS can improve the upper-level water substances and further improve the precipitation forecast. In this case, both experiments slightly underestimated the amount of precipitation at the 50 mm levels.

Fig. 11. Time series of the 3-h accumulated rainfall (mm) for (a, d, g, j) the observation, (b, e, h, k) the Real_Kessler experiment,and (c, f, i, l) the Real_WSM6 experiment during the 12-h period of 0000 UTC 25 July 2018 to 1200 UTC 25 July 2018.

Fig. 12. (a) The observed 12-h accumulated rainfall (mm) from 0000 UTC to 1200 UTC on 25 July 2018 and the corresponding forecast precipitation from (b) the Real_Kessler experiment and (c) the Real_WSM6 experiment.

Fig. 13. The 12-h accumulated precipitation forecast scores for the thresholds of 1, 5, 10, 25, and 50 mm h—1 for the (a) threat score (TS) and the (b) bias score (BIAS).

5. Summary and Discussion

The main objective of this study was to develop the capability of including independent water vapor, cloud water,rain, ice, snow, and graupel water substances in a WRF incremental 4D-Var for convective-scale data assimilation. A relatively widely used single-moment MPS, WSM6, was chosen for regularization, and then its TLM and ADM were developed in WRFPLUS. Then, single-observation and real-sounding data assimilation experiments were conducted to validate its correctness.

A logistic function was used to smooth the water substances and temperature discontinuities of the on-off process in WSM6. The difference parameter k was used for temperature and moisture processing. The original WSM6 was also used to compare and validate the rationality of the linearization process. The forecast experiment with the regularized WSM6 predicted precipitation well, and the forecasted hydrometeors were consistent with the original WSM6. The TLM and ADM of RWSM6 were developed using TAPENADE, and their interfaces in the minimization algorithm were added to WRFDA and WRFPLUS. The verification results showed that the TLM and ADM of RWSM6 were correct.

A single and real-sounding data assimilation experiment was used to validate the rationality and stability of the iterative convergence in the assimilation system. As a reference, the results of the warm rain Kessler scheme are also presented for comparison. The results of the single-point experiment revealed flow-dependent characteristics in the updated WRF 4D-Var system, and the vertical structure and horizontal distribution of the assimilation of different water substances were considered reasonable. In the real-sounding data assimilation experiment, data from 18 sounding stations were assimilated in the WRF 4D-Var system. The results indicated that the updated WRF 4D-Var, with the cloud, water, rain, ice, snow,and graupel physics of RWSM6, produced reasonable analyses and forecasts and did not increase the computing resource requirement too severely in the iteration process.

The real-sounding data assimilation experiments revealed that the water substances under the initial conditions were more reasonably updated with WSM6 than with the Kessler scheme. The initial Qvaporhad a greater effect on the redistribution of water substances for both experiments with Kessler and WSM6 because its magnitude was larger than that of the other water substances. From the results of the precipitation forecast, this cold cloud assimilation scheme with water vapor, cloud water, rain, ice, snow, and graupel could further improve the prediction of precipitation in the WRF model.

This work mainly focused on the development of the cold cloud assimilation scheme in WRF 4D-Var. Single-observation and real-sounding data assimilation experiments were set up to test the correctness and rationality of updating the initial water substance in 4D-Var. A large number of assimilation experiments, types of data, and individual cases are needed to further verify the applied effect of this cold cloud scheme. The next step is to assimilate more types of observational data (e.g., radar reflectivity and satellite radiation data) at the convective scale in 4D-Var to verify the correctness and effectiveness of the cold cloud scheme.

Acknowledgements.This work was supported by the National Key R & D Program of China under Grant Nos.2018YFC1507302 and 2018YFC1506803, National Natural Science Foundation of China No. 42275171, the Liaoning Province Key R& D Program of Liaoning of China under Grant No. 2020JH2/10300091, and the Bohai Rim Meteorological Science Collaborative Innovation Fund under Grant No. QYXM201901.