Overview and Prospect of Data Assimilation in Numerical Weather Prediction

1.   Introduction

As an initial value problem, numerical weather prediction (NWP) can be achieved by advancing a numerical model, given the current state of the atmosphere, and associated lateral boundary, and top and bottom boundary conditions. More accurate initial conditions can lead to improved numerical weather forecasts, while the initial conditions are the best possible estimate of the atmospheric state using all available information (Talagrand, 1997). Data assimilation (DA) is a sophisticated process to combine the noisy observations with uncertain short-term forecasts, resulting in the optimal estimate of the atmospheric state (Daley, 1991; Kalnay, 2003; Zou, 2025).

In the early stages of DA for weather forecasting, empirical methods were developed, such as the interpolation schemes (Panofsky, 1949; Gilchrist and Cressman, 1954), successive correction method (Cressman, 1959), and Newton relaxation method (Hoke and Anthes, 1976). At the same time, Chinese scholars innovatively proposed to use recent historical weather data for future state forecasting, which transformed the NWP initial value problem into an extrapolation forecast process based on historical weather evolutions (Koo, 1958a, b). Chou (1974) further proposed using functional extreme value problem with multi-time historical observations, as an equivalent way to solve the differential equations that approximately describe the atmospheric processes, which is the primary concept of the variational methods later developed and widely used in NWP. As numerical models advanced, Chinese scholars deepened the theoretical basis of DA, pointing out that NWP is an initial value problem and also an inverse problem (Chou, 2007), through which the initial conditions, boundary conditions, and model parameters can be optimally estimated. But challenges remain due to the ill-posed problems, which could be solved by introducing regularization from inverse problems to DA and adding stabilization functionals to the objective functional (Huang et al., 2003). The China Meteorological Administration (CMA) operational system initially relied on the imported DA technology. Since 2000, CMA has been dedicated to developing variational DA methods and systems, gradually establishing the operational regional three-dimensional variational DA system (GRAPES-MESO, now CMA-MESO) and global four-dimensional variational DA system (GRAPES, now CMA-GFS) (Xue and Chen, 2008; Zhang et al., 2019; Shen et al., 2020), which significantly improves the accuracy of global operational forecasts (Fig. 1).

Download

Larger image

PowerPoint slide

Fig  1.  Progress of global NWP by CMA. Monthly mean evolution of anomaly correlation coefficients (ACC) at 500-hPa geopotential height for days 3, 5, and 7, respectively, from January 2010 to August 2024 (solid lines for northern hemisphere, and dashed lines for southern hemisphere).

In addition to the fundamental DA theories and methods, the effective and efficient integration of observation data into DA systems is crucial for NWP. The assimilation of satellite retrievals into numerical models began in the 1970s, marking a significant contribution to the field of numerical forecasting (Smith et al., 1970). Since the 1990s, the direct assimilation of satellite radiances has been made possible through advancements in fast radiative transfer models and variational assimilation methods (Saunders et al., 2018; Weng et al., 2020; Johnson et al., 2023). This development has further enhanced the role of satellite data in NWP (Eyre et al., 2020). In 2009, China’s global/regional integrated numerical forecast model (GRAPES) was quasi-operationally implemented, significantly advancing the utilization of satellite data, particularly from the Fengyun (FY) series. In 2015, FY-2D cloud motion vectors were operationally assimilated into GRAPES, followed by the operational assimilation of FY-3C microwave radiances and occultation data in 2016 (Li and Liu, 2016; Li G. et al., 2016; Li J. et al., 2016). Since then, a variety of satellite observations from FY-4A/B and FY-3D/E have been operationally assimilated, marking a significant milestone in the quantitative application of FY satellite data. Moreover, radar data has played an important role for monitoring and forecasting the convective-scale weather systems. Chinese researchers have demonstrated the effectiveness to assimilate radar data in improving the accuracy of forecasts for convective-scale weather systems (Chen M. et al., 2014; Shao et al., 2016; Chen Y. D. et al., 2018; Sun et al., 2020b).

The concept of target observation (i.e., adaptive observation) strategy has refined DA by collecting and assimilating additional high-quality observations in limited “sensitive areas”, in order to improve the initial conditions and then the subsequent forecasts of high-impact weather events. To identify the sensitive areas for target observations, Chinese scholars considered the nonlinear nature of atmospheric and oceanic motions, and proposed the conditional nonlinear optimal perturbation (CNOP) method. The additionally collected observations based on the CNOP-identified sensitive areas have been proved to more effectively improve the forecast of high-impact weather events, compared to those collected upon sensitive areas identified by traditional linear approximation methods.

This paper reviews the development of DA theories and methods, the progress in assimilating multi-source observations, the development of targeted observation strategies, and advances in China’s operational DA systems for NWP. The challenges and opportunities for DA in the context of rapid development of numerical models, observing systems, data science, and artificial intelligence are also discussed.

2.   Theories and methods of DA

Along with the development of numerical models, observing systems, and computational science, the DA theories and methods have evolved from early empirical objective analysis to analysis theories based on statistics, and then to assimilation methods that incorporate atmospheric dynamics. This section focuses on the DA methods that are supported by statistical theories and have been operationally used in NWP. The nonlinear DA and coupled DA methods that are yet to be operationally applied are discussed in the prospect.

2.1   Variational methods

Analysis methods with statistical foundations begun to develop around the 1980s. Eliassen et al. (1960) first derived the multivariate optimal interpolation equations based on observations and background fields. Subsequently, the optimal interpolation (OI) method was proposed, which seeks the optimal weight matrix in the physical space, such as at grid points (McPherson et al., 1979) or over finite volume elements (Lorenc, 1981). OI utilizes estimated background error covariances based on the differences between short-term forecasts and radiosonde observations (Hollingsworth et al., 1986; Thiébaux and Pedder, 1987). Gandin (1963) independently derived the multivariate OI equations and applied them to objective analysis in the Soviet Union. OI became the operational analysis method from the 1980s to the early 1990s.

Different from the OI that locally updates the weights given an “influence radius”, the three-dimensional variational method (3DVar; Sasaki, 1970) uses global optimization algorithms to minimize the cost function and directly obtains the minimum of control variables. 3DVar also has an observation-space form known as the physical-space statistical analysis system (PSAS; Da Silva et al., 1995), which seeks the minimum of the cost function in the observation (physical) space. Although 3DVar, PSAS, and OI are equivalent in terms of their solutions (Lorenc, 1986), 3DVar and PSAS use more general and global background error covariances than OI, such as those based on forecast differences at the same forecast time (Parrish and Derber, 1992; Rabier et al., 1998).

As an extension of 3DVar, four-dimensional variational method (4DVar; Lewis and Derber, 1985; Courtier and Talagrand, 1990) considers temporal distributions of observations within an assimilation window (Daley, 1991) and can implicitly account for temporal evolutions of background error covariances (Thépaut et al., 1993). To efficiently obtain the optimal solution, 4DVar can be expressed in an incremental form (Courtier et al., 1994; Lorenc, 1997), by which the optimal perturbation relative to a reference state rather than the whole optimal state is solved, and the solving process can be accelerated through the “preconditioning” (Parrish and Derber, 1992; Derber and Bouttier, 1999). Strong-constrain 4DVar assumes the model being perfect (Sasaki, 1970), and the perfect-model assumption can be relaxed by model error corrections (Derber, 1989; Zupanski, 1993) or model error representation using weak constraints (Bennett, 1992; Egbert et al., 1994; Bennett et al., 1996). When the model is perfect and the background error covariances at the initial time are accurate, the analysis of 4DVar at the end of the assimilation window is equivalent to that of the generalized Kalman filter (Lorenc, 1986). Since the mid-1990s, 3DVar and 4DVar have become mainstream operational DA methods.

2.2   Ensemble Kalman filters

Both OI and 3DVar methods use static background error covariances, but “errors of the day” reveal the importance of flow-dependent background error covariances (Kalnay et al., 1997). Kalman filter (KF, Kalman, 1960; Kalman and Bucy, 1961) uses background error covariances that evolve with the numerical model over time, which establishes the mathematical framework for four-dimensional DA. Extended Kalman filter (EKF, Ghil et al., 1981; Daley, 1995), as an extension of KF to nonlinear models, can provide the best linear unbiased estimate (BLUE) for the state and its error covariances. EKF is considered the “gold standard” of DA, but it requires massive computations due to the update of background error covariances based on linear model matrices.

Ensemble Kalman filter (EnKF, Evensen, 1994; Houtekamer et al., 1996) uses the Monte Carlo method to estimate the background error covariances based on samples of ensemble forecasts, which can be seen as a simplified EKF. EnKF can approximate the KF’s analysis solution and is suitable for high-dimensional dynamical systems, while providing ensemble initial conditions for subsequent ensemble forecasts (Houtekamer et al., 2005, 2014). The original EnKF is a stochastic one, by use of perturbed observations for each ensemble member to achieve consistent analysis and associated error covariances (Burgers et al., 1998; Houtekamer and Mitchell, 1998). To avoid sampling errors caused by perturbing observations, deterministic EnKFs have been proposed, such as the ensemble adjustment Kalman filter (EAKF, Anderson, 2001), ensemble square-root filter (EnSRF, Whitaker and Hamill, 2002), and local ensemble transform Kalman filter (LETKF, Bishop et al., 2001; Hunt et al., 2007). Deterministic EnKFs solve for optimal Kalman gains using the analysis error covariances and achieve equivalent solutions without covariance localization (Tippett et al., 2003). As an extension of EnKF, ensemble Kalman smoother (EnKS, Evensen and van Leeuwen, 2000) further assimilates future observations to update the analysis of EnKF, based on temporal sample correlations.

EnKF faces the challenges of filter divergence, especially when it is applied to high-dimensional dynamical systems, due to limited ensemble sizes, model errors, and linear correlations. One way to combat the filter divergence is covariance localization (Hamill, 2001; Anderson, 2012), which is typically a function of the distance between observations and model variables (Gaspari and Cohn, 1999). Covariance localization can be implemented through the background error covariance matrix (Houtekamer and Mitchell, 2001; Lei et al., 2018) or the observation error covariance matrix (Hunt et al., 2007). The localization function varies with observation types (Zhang et al., 2009a; Lei and Anderson, 2014b), state variable kinds (Kang et al., 2011; Lei et al., 2015), and assimilation times (Anderson, 2007; Chen and Oliver, 2010). Thus, adaptive localization methods have been developed (Bishop and Hodyss, 2009; Lei and Anderson, 2014a; Zhen and Zhang, 2014; Flowerdew, 2015; Lei et al., 2020). Another approach to address the filter divergence is covariance inflation (Anderson and Anderson, 1999; Houtekamer and Mitchell, 2005). It can be implemented by empirically or adaptively multiplying the ensemble perturbations (Anderson, 2009; Miyoshi, 2011; El Gharamti, 2018), increasing posterior ensemble perturbations relative to prior ensemble perturbations or prior ensemble spread (Zhang et al., 2004; Whitaker and Hamill, 2012; Ying and Zhang, 2015), or augmenting additive noises (Wang X. G. et al., 2013; Yang et al., 2015), particularly accounting for model uncertainties (Buizza et al., 1999; Berner et al., 2009; Ha et al., 2015; Zeng et al., 2020). Since the early 2000s, EnKF has been operationally used at Environmental Canada and National Centers for Environmental Prediction (NCEP) (Houtekamer and Mitchell, 2005; Whitaker et al., 2008).

2.3   Hybrid ensemble-variational methods

The static background error covariance matrix used by the variational methods is full rank, but unable to capture the “errors of the day”. On the other hand, EnKF constructs flow-dependent background error covariance matrix using short-term ensemble forecasts, but the flow-dependent one is often rank-deficient and affected by sampling errors and model errors. Thus, hybrid ensemble-variational methods that combines the advantages of variational methods and EnKFs have been developed. The hybrid ensemble-variational method can directly combine the static and flow-dependent background error covariance matrices (Hamill and Snyder, 2000), to mitigate the impact of rank deficiency and sampling errors in the flow-dependent background error covariances (Wang X. G. et al., 2008, 2013; Zhang et al., 2009b; Kleist and Ide, 2015a, b). The static background error covariances can also be incorporated into the EnKF to represent model errors (Meng and Zhang, 2008). Ensemble-4DVar (En4DVar; Lorenc, 2003; Bonavita et al., 2012) embeds the flow-dependent background error covariances into the cost function of the variational method through the alpha control vector, which is equivalent to directly combining the static and flow-dependent background error covariances (Wang et al., 2007). En4DVar outperforms either standalone 4DVar or EnKF (Zhang et al., 2009b; Buehner et al., 2010a, b).

Compared to En4DVar, the 4D ensemble-variational method (4DEnVar, Liu et al., 2008) captures the temporal evolution of error covariances based on ensemble forecasts, eliminating the need for tangent-linear and adjoint models. Similar to 4DEnVar, the dimension-reduced projection 4DVar (DRP-4DVar, Wang et al., 2010; He et al., 2017) and the nonlinear least-squares ensemble 4DVar (NLS-En4DVar, Tian and Feng, 2015; Tian et al., 2018) have been proposed. However, 4DEnVar cannot account for the evolution of static background error covariances within the assimilation window (Wang and Lei, 2014) and struggles to handle time-varying localization (Bishop and Hodyss, 2009). Consequently, 4DEnVar is inferior to En4DVar (Lorenc et al., 2015; Poterjoy and Zhang, 2015, 2016).

Unlike the hybridization of static and flow-dependent background error covariances, the hybrid gain approach (Penny, 2014) mixes the analyses from the variational method and EnKF. The hybrid gain approach also outperforms the standalone EnKF and 4DVar with static background error covariances (Bonavita et al., 2015). Since the early 2010s, the hybrid ensemble 4DVar has been operational at centers such as European Centre for Medium-Range Weather Forecast (ECMWF) and the Met Office (Bonavita et al., 2012; Clayton et al., 2013). 4DEnVar has later been implemented at operational centers of Canada, the U.S. and so on (Buehner et al., 2015; Caron et al., 2015; Kleist and Ide, 2015b).

Hybrid ensemble-variational methods require a variational system and an EnKF system, but inconsistencies between the two systems could result in suboptimal analyses. The ensemble variational integrated localized method (EVIL; Auligné et al., 2016) constructs the ensemble analyses using the analysis error covariances from the variational solution, to avoid the need for an ensemble framework. The integrated hybrid ensemble-variational method (IHEnKF; Lei et al., 2021) approximates the static background error covariances by a large size of climatological perturbations, and achieves the solutions of hybrid ensemble-variational and hybrid gain methods within a pure ensemble framework. Moreover, IHEnKF can update the ensemble perturbations by the hybrid background error covariances, which leads to superior ensemble analyses than the traditional hybrid ensemble-variational methods.

2.4   Multi-scale and balanced DA

The atmospheric state and its evolution span multiple spatial and temporal scales, and multi-source observations also capture information cross different scales. Thus, multiscale DA methods have been proposed, to effectively use the multi-source observations to constrain the multi-scale atmospheric state. Multiscale DA methods can be iteratively and sequentially implemented. Observations representing large scales (e.g., conventional observations) are first assimilated using a broad localization lengthscale; then small-scale observations (e.g., radar data) are assimilated with a tight localization lenghscale, aiming to effectively extract information from observations at different scales (Zhang et al., 2009a; Xie et al., 2011; Sodhi and Fabry, 2022). Alternatively, all observations can be assimilated using different localization lengthscales, and the final analysis is constructed by combining the analysis increments with various localization lengthscales (Miyoshi and Kondo, 2013). Different from the iterative assimilation methods, multiscale DA can be performed in a single step within the ensemble-variational framework or pure ensemble assimilation framework. Different localization lengthscales are applied to the background error covariances at different scales, and then the multi-scale state variables are simultaneously updated by the multi-source observations across all resolvable scales (Buehner, 2012; Buehner and Shlyaeva, 2015; Wang X. G. et al., 2021; Wang and Wang, 2023).

When the numerical models use the analyses produced by DA methods as initial conditions to advance, they face “insertion noises” or “initialization shocks”. Thus, balanced DA is required to minimize the insertion noises caused by unbalanced initial conditions, which could results in spurious gravity waves and adversely affect subsequent forecasts (Temperton and Roch, 1991). It is straightforward for 4DVar to include balance constraints in the cost function. However, intermittent EnKF faces the issue of imbalances, and then continuous EnKF that transforms the intermittent EnKF to a continuous form has been proposed (Bergemann and Reich, 2010; Lei et al., 2012). Moreover, to mitigate the initialization shocks, EnKF can utilize shorter assimilation windows (He et al., 2020; Slivinski et al., 2022), or leverage additional current observations to obtain more accurate past initial conditions (Kalnay and Yang, 2010). There have been general initialization methods, such as the digital filter that eliminates rapid oscillations (Lynch and Huang, 1992), the incremental analysis update that preserves large-scale analysis increments (Bloom et al., 1996), and the four-dimensional incremental analysis update that retains the evolution of both large- and small-scale increments within the assimilation window (Lorenc et al., 2015; Lei and Whitaker, 2016).

3.   Multi-source atmospheric observations

3.1   Satellite data

Over the past five decades, the role of satellite DA in NWP has grown significantly. Following the launch of the “Nimbus 3” satellite in April 1969, which carried the first temperature detector, the retrieval products from the Satellite Infrared Spectrometer (SIRS) were first “incorporated” into the objective analysis of the U.S. National Meteorological Center and had a significant impact on the analyses in the Pacific region and the forecasts across U.S. (Smith et al., 1970). Subsequently, a series of international experiments on satellite product assimilation were conducted (Atkins and Jones, 1975; Desmarais et al., 1978; Druyan et al., 1978; Kelly et al., 1978; Gilchrist, 1982; Uppala et al., 1984). However, due to the of low vertical resolution of satellite data and relatively large temperature retrieval errors of 2–3 K, the assimilation experiments during this period generally had a neutral impact on forecasts, with a more pronounced effect in the southern hemisphere (Ohring, 1979). Eyre and Lorenc (1989) pioneered the direct assimilation of the radiative brightness temperature observed by satellites in NWP and successfully integrating the radiative information from the Television Infrared Observation Satellite Program (TIROS) vertical sounding through one-dimensional variational assimilation (Eyre et al., 1993). In October 1995, the NCEP (Derber and Wu, 1998) and in January 1996, the ECMWF (Andersson et al., 1994; McNally and Vesperini, 1996; Saunders et al., 1997) led the way in directly assimilating satellite radiances using 3DVar. ECMWF further advanced this by adopting direct satellite DA in 4DVar in November 1997. Subsequently, other operational centers successively directly assimilated radiances in 3DVar and 4DVar systems (Chouinard et al., 2002; Joo and Lee, 2002; Okamoto et al., 2002; Li J. et al., 2016).

3.1.1 Fast radiative transfer model

The fast radiative transfer model quickly maps the atmospheric state variables to the observed quantities, serving as an observation forward operator. It mainly consists of an atmospheric gas absorption module, a particle scattering module, a surface emissivity module, a radiative transfer solution module, and the corresponding tangent linear and adjoint modules. Currently, three fast radiative transfer models are widely used in satellite DA for NWP: the Radiative Transfer for the TIROS Operational Vertical Sounder (RTTOV) model developed by the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) (Saunders et al., 2018), the Community Radiative Transfer Model (CRTM) model developed by the National Oceanic and Atmospheric Administration (NOAA) of the U.S. (Johnson et al., 2023), and the Advanced Radiative Transfer Modeling System (ARMS) model developed by the CMA (Weng et al., 2020; Yang et al., 2020). ARMS, a fast radiative transfer model independently developed by China, has replaced RTTOV in CMA-GFS since 2023. Its innovations are mainly reflected in the following aspects. (1) A calculation scheme for atmospheric transmittance coupled with the real spectral response function, effectively enhancing the simulation accuracy of microwave channels (Kan et al., 2024). (2) A scattering database for non-spherical cloud particles and aerosol particles based on the T-Matrix and the Discrete Dipole Approximation (DDA) methods, supporting all-sky satellite DA (Yang et al., 2020). (3) A polarized Bidirectional Reflectance Distribution Function (pBRDF) based on a two-scale ocean roughness model, improving the physical mechanism of the bidirectional reflection at the sea–air interface and enhancing the simulation of the satellite radiative transfer with integrated active and passive sensors (He and Weng, 2023). (4) An improved physical model of microwave land surface emissivity (LandEM) and the development of the Chen–Weng rough surface reflectivity model, increasing the estimation accuracy of complex surface emissivity (Liu et al., 2024). (5) Refinement of the discrete ordinate radiative transfer theory, overcoming assumptions and dependences on the azimuthal symmetry property in atmospheric scattering and surface reflection process, and developing a general vector radiative transfer solution scheme (VDISORT) for simulating Stokes vector radiation of satellites across the full spectrum range (Zhu et al., 2024).

3.1.2 Infrared radiance DA

Due to the spectral limitations and challenges in assimilating radiative quantities in cloud areas (Li et al., 2022a), infrared radiance DA mainly focuses on the direct assimilation of clear-sky radiances or radiances with partial cloud cover. For clear-sky radiances, precise cloud detection is essential. The clear-sky channel cloud detection scheme sorts the channels to determine cloud top height and assimilates channels above the cloud top, improving utilization of the satellite data (McNally and Watts, 2003). For assimilating infrared sounding with partial cloud cover, Li et al. (2005) proposed the “optimal cloud clearing” technique, converting partially cloud-covered infrared sounding into equivalent clear-sky radiative quantities, effectively improving the utilization rate of infrared sounding data in rain and cloud areas and improving tropical cyclones forecasts (Wang P. et al., 2014, 2017). To address challenges in directly assimilating the infrared radiances in rain and cloud areas, Jones et al. (2013) and Chen et al. (2015) successfully assimilated satellite cloud water and cloud ice path products retrieved from visible and near-infrared soundings by using the EnKF and variational methods, respectively. Meng D. M. et al. (2019) introduced hydrometeors into the extended control variables and achieved the hybrid assimilation of infrared retrievals in cloudy areas.

In addition, techniques such as channel correlation processing and principal component analysis have been developed to optimized channel assimilation, addressing the computational costs and spectral-related observational errors of infrared hyperspectral data (Rabier et al., 2002; Collard, 2007; Matricardi and McNally, 2014; Zhou et al., 2024). Typically, about 200 channels per instrument are assimilated, which also avoids spectral regions affected by trace gases like ozone. CMA-GFS currently has the ability to assimilate the infrared hyperspectral data of China’s polar-orbiting and geostationary satellites, such as FY-3D/E HIRAS (Liu and Xue, 2014), FY-4A/B GIIRS (Yin et al., 2020, 2021; Han et al., 2023), and is also capable of assimilating the infrared hyperspectral data such as METOP-B/C IASI (Li G. et al., 2016) and NOAA 20 CRIS in real time. In addition, infrared imager data, including FY-2 VISSR, FY-4A/B AGRI (Wang et al., 2018), H8/H9 AHI and GOES-18 ABI data, etc., have also achieved operational application in CMA-GFS.

3.1.3 Microwave radiance DA

Among numerous satellite instruments, microwave sounding data can penetrate the cloud and rain, and provide information on the vertical distribution of atmospheric temperature and water vapor over the whole sky and the entire surface, significantly improving the forecasting accuracy (Bormann et al., 2019; Li et al., 2024; Luo et al., 2025). In 2009, ECMWF implemented the world’s first operational assimilation system for all-sky satellite data (Bauer et al., 2010). Subsequently, satellite DA techniques in rain and cloud areas were successively applied to the operational models of the Japan Meteorological Agency and NCEP (Okamoto, 2014; Zhu et al., 2014). Initially, ECMWF adopted an indirect assimilation strategy of 1D–4DVar for rain and cloud areas. For satellite observations affected by rain and clouds, the temperature and humidity in the background field were used as the first guess values. Through 1DVar, the total water vapor content (TWPC) corresponding to satellite data was retrieved, and then TWPC was taken as a virtual observation and assimilated by 4DVar (Geer et al., 2008). Since 1DVar can retrieve the atmospheric state matching the rain and cloud conditions of satellite observations, it can avoid the mismatch between the background and observations under the cloudy and rainy conditions. Meanwhile, it can also conduct all-sky quality control through 1DVar without the need to adopt a complex cloud detection scheme. However, the TWPC retrieved by this method has already implied the humidity information of the background field, and false observation increments will be generated when it is subsequently brought into 4DVar to update the background (Geer et al., 2010). Therefore, this method was later replaced by the direct all-sky assimilation at ECMWF (Bauer et al., 2010; Geer et al., 2010). CMA-GFS has also successfully carried out the all-sky assimilation of FY satellite microwave imager, significantly improving global water vapor analyses and forecasts (Xie et al., 2023). In the past decade, the most important progress in satellite DA in NWP has been the DA techniques under the influence of clouds and precipitation (Geer et al., 2018).

Microwave radiance assimilation affected by surface conditions has also gained attention. Accurate surface emissivity estimation is crucial, with methods including physical model method, statistical model method, and dynamic inversion method (Tian et al., 2015). The physical model method aims to establish a numerical model based on the physical relationship between surface emissivity and various surface parameters (such as surface type, vegetation parameters, soil parameters, etc.). However, its calculation accuracy depends on the accurate estimation of a large amount of input information of surface parameters, which is difficult to obtain, hindering the wide application of this method (Weng et al., 2001). The statistical model method refers to using the historical data set of surface emissivity as the empirical or semi-empirical estimated value of the actual instantaneous emissivity, such as the TELSEM data set (Aires et al., 2011) and the CNRM data set (Karbou et al., 2010), etc. This type of method is easy to use, but the historical statistical data set is difficult to represent the current instantaneous emissivity state and cannot estimate the temporal changes of surface features or spatial changes of surface features at the sub-pixel scale. The dynamic inversion method refers to obtaining the surface emissivity given observed brightness temperature by calculating the upward and downward radiation of the atmosphere, as well as the atmospheric transmittance and surface temperature (Karbou et al., 2006). This method not only has strong usability but also can consider the dynamic changes of emissivity under complex conditions, so it has been widely used in the assimilation of land surface microwave data (Krzeminski et al., 2009; Baordo and Geer, 2016; Xiao et al., 2023b). Based on the dynamic inversion method, the CMA-GFS model has successfully achieved the operational assimilation of the near-surface channels of the AMSU data on land, significantly improving global lower atmosphere analyses and forecasts, particularly in the northern hemisphere (Xiao et al., 2023b).

3.1.4 From dual-satellite constellation to three-satellite constellation

Joo et al. (2013) found that 64% of the reduction in numerical forecast errors of NWP was contributed by satellite observations, with polar-orbiting meteorological satellites contributing approximately 90% of this reduction. Before 2010, the international polar-orbiting meteorological satellites operated in a dual-satellite constellation (“AM” satellite and “PM” satellite). Recognizing the limitations of the dual-satellite system could not provide complete global coverage within the 6-h assimilation window of the global NWP model, the World Meteorological Organization (WMO) proposed a three-satellite constellation model (“Dawn”, “AM”, and “PM”) in 2009. In 2014, the CMA clearly stated the plan to launch the “Dawn” satellite in the feasibility study report for the third batch of FY-3 satellites. In 2021, FY-3E was successfully launched, and the Chinese researchers realized that the observation system of the three-satellite constellation could effectively make up the observation gap of polar-orbiting satellites within the 6-h assimilation window (Zhang et al., 2022). The observation data of FY-3E have been applied not only in China’s NWP operational system (Li et al., 2024), but also in the operational models of many organizations such as the ECMWF, the Met Office, the Japan Meteorological Agency, and the Korea Meteorological Administration ensuring global observation needs for NWP (Zhang et al., 2024).

3.2   Radar data

Radar data has high temporal and spatial resolutions, allowing it to capture fine information for convective-scale weather systems. Proper utilization of radar data can significantly improve the dynamical and microphysical characteristics of convective weather systems in the initial conditions. Thus, the effective assimilation of radar observations is one of the key factors in improving convective-scale NWP (Wan et al., 2005; Sun et al., 2014; Sun et al., 2020b).

3.2.1 Radar observation operators

Traditional Doppler weather radar detect variables mainly including the radial wind and reflectivity. Radial wind represents important dynamical features of the internal structure of convective-scale weather systems, while reflectivity contains microphysical information about these systems. Radar radial wind has been widely used in various DA systems, including the VDRAS 4DVar system (Sun and Crook, 1997), ARPS 3DVar system (Gao et al., 1999), WRFDA system (Xiao et al., 2005), and GRAPES 3DVar system (Liu et al., 2010). But conventional radial wind observation forward operators only introduce information of radial wind, and could be insufficient to analyze tangential wind. Based on the assumption of a uniform wind field within the regional azimuth, Luo et al. (2014) incorporated radar radial wind speed and its spatial variations into the radial wind forward operator, enabling the analysis of tangential wind information during radial wind assimilation. This radial wind forward operator was introduced into both the GSI (Chen et al., 2017) and GRAPES (Ma et al., 2016) assimilation systems.

Compared to radial wind forward operators, reflectivity observation forward operators are more complex. Early studies established reflectivity forward operators based on the empirical relationship between reflectivity and rainfall (Sun and Crook, 1998; Xiao and Sun, 2007). Tong and Xue (2004) and Gao and Stensrud (2012) further incorporated ice-phase particles such as snow and hail, developing reflectivity forward operators based on the Lin microphysical parameterization scheme, and implemented direct assimilation of reflectivity within the EnKF and variational frameworks. Similarly, Hawkness-Smith and Simonin (2021) constructed reflectivity forward operators to achieve direct assimilation of reflectivity within the Met Office assimilation system. Liu et al. (2022) introduced raindrop number concentration based on the Thompson microphysics scheme, developing a reflectivity forward operator for 2-moment hydrometeors. However, reflectivity observation forward operators remain somewhat empirical, and the associated errors are relatively large.

Jung et al. (2008b) estimated the backscatter cross-section parameters of various precipitation particles using the T-matrix algorithm, developing a more accurate forward reflectivity forward operator and achieving successfully assimilated reflectivity in an EnKF. Following Jung et al. (2008a), Wang and Liu (2019) developed the tangent linear and adjoint operators to perform direct reflectivity assimilation within a variational framework. Zeng et al. (2013; 2014) and Jerger et al. (2013) incorporated different physical aspects into the reflectivity forward operator, developing a three-dimensional reflectivity forward operator. Wang and Liu (2019) further developed a forward reflectivity observation operator for ice-phase particles based on Jung et al.(2008a), parameterizing reflectivity as a rapid polynomial relationship for the mixed ratio. These more complex and accurate reflectivity forward operators significantly reduce the errors from forward operators, leading to better assimilation of radar reflectivity.

In addition, many studies have focused on developing radar observation forward operators. Jung et al. (2008a) implemented a polarization radar data simulator, which uses spherical particles to represent hydrometeors and calculates spectral properties with either online Rayleigh approximations or offline lookup tables (Mishchenko and Travis, 1994). This polarization radar data simulator has been applied to low-frequency S-band, C-band, and X-band radar. Wolfensberger and Berne (2018) developed a cross-platform polarimetric radar observation forward operator, which includes hail particle radar simulations and represents all hydrometeor particles as uniform spherical bodies. Oue et al. (2020) developed a cloud-resolving model radar simulator that can simulate both polarimetric radar and lidar observations, though it is currently limited to ground platforms and does not explicitly handle melting particles. Zhejiang University (ZJU-AERO, Xie et al., 2024) designed an accurate and efficient radar observation forward operator that incorporates scattering calculations for hydrometeors and constructs an optical property database, allowing it to handle non-spherical and inhomogeneous hydrometeor particles in the atmosphere.

3.2.2 Conventional radar DA

Radial wind observations contain important dynamical features of convective-scale weather systems and provide crucial tools for monitoring and studying convective-scale weather systems (Xu, 2003; Liang, 2007; Yang et al., 2008). Assimilation of radar radial wind is relatively mature and can significantly improve analyses and forecasts of convective-weather systems (Gao et al., 2004; Li et al., 2012; Zhu et al., 2013; Chen et al., 2014; Shao et al., 2016; Chen et al., 2019; Mu et al., 2019; Chen et al., 2025).

Comparing to assimilation of radial winds, the assimilation of reflectivity is more complex. Current methods of radar reflectivity assimilation are generally divided into two categories: direct and indirect assimilation. Direct assimilation projects model state variables into the observation space, directly comparing the background with the observed reflectivity. The innovation and associated uncertainties are used for assimilation and leads to the analysis. Direct assimilation of radar reflectivity has been applied effectively (Sun and Crook, 1997; Tong et al., 2005; Sheng et al., 2006).

However, direct assimilation of radar reflectivity using variational methods also faces challenges. When the background hydrometeor content is low, the gradient of the observation term in the cost function can be large, which often prevents effective convergence of the minimization. The nonlinear reflectivity forward operator is difficult to construct, often resulting in unrealistic hydrometeor analyses (Wang et al., 2013a, b; Liu et al., 2019). Incremental variational assimilation improves the nonlinearity of the reflectivity forward operator through multiple outer loops, but with increased computational cost. Compared to variational methods, EnKFs can use nonlinear forward operators of radar reflectivity (Lan et al., 2010a, b; Yussouf and Stensrud, 2010), effectively use reflectivity observations with complex microphysical processes. However, nonlinear forward operators do not satisfy the assumption of Gaussian error distributions of EnKFs, and thus, suboptimal solutions are obtained, especially for dense and strongly nonlinear reflectivity observations. Moreover, the EnKF is also affected by sampling errors and model errors (Liu et al., 2020).

Background error covariances play a crucial role in convective-scale DA, and appropriate background error covariances can lead to more coherent analyses (Chen et al., 2013; Chen et al., 2022）; Zheng et al., 2023). Wang and Wang (2021) developed static background error covariances including hydrometeor control variables, which is incorporated into a hybrid ensemble-variational framework, leading to improved supercell predictions than ensemble-based background error covariances. Furthermore, assimilating radar reflectivity with hydrometeor-included background error covariances can significantly improve thermodynamic conditions and heavy rainfall forecasts, due to the vertical and multivariable correlations (Zheng et al., 2023).

To avoid linearization errors caused by the nonlinear reflectivity forward operator during direct assimilation, many studies and operational systems often use indirect assimilation for radar reflectivity. In indirect assimilation, the radar reflectivity is first inverted into model variables during the assimilation process, where the type and proportion of hydrometeors are determined using the prior temperature, and then these inverted model variables are assimilated (Wang et al., 2013a). Some studies have proposed a new background-dependent hydrometeor inversion method, which updates the type and proportion of hydrometeors in real time based on the background characteristics. This improves reflectivity assimilation and enhances weather forecasts (Chen et al., 2020, 2021; Huang et al., 2022). The indirect assimilation method avoids constructing tangent linear and adjoint operators for reflectivity observation operators, improving the mathematical conditions for solving the cost function, while having relatively lower computational costs compared to direct assimilation. As a result, it is widely used in various studies and operational systems (Fan et al., 2013; Lai et al., 2020; Zhang et al., 2019).

3.2.3 Dual-polarization radar data assimilation

In recent years, several countries, including China, have started upgrading their dual-polarization radar networks (Wu et al., 2018). Dual-polarization radar can improve the identification the phase state characteristics of hydrometeors, and the effective use of dual-polarization radar data can improve microphysical initial conditions, such as hydrometeor content (Zhao et al., 2019). Assimilation of dual-polarization radar observations has made progress. Chinese researchers have developed dual-polarization observation forward operators based on single- and double-moment microphysical schemes and directly assimilated simulated dual-polarization radar observations using EnKFs (Jung et al., 2008a, b; Xue et al., 2010). Further, Putnam et al. (2019) used and EnKF to assimilate real dual-polarization radar observations, leading to improved analyses and forecasts of polarimetric quantities, although their assimilation was limited to horizontal reflectivity factor and differential reflectivity below 2 km.

The EnKF is limited by the finite ensemble size, which makes it difficult to properly estimate the background error covariances with rank deficiencies, faces the challenges of imbalances in the analyses and model errors. Consequently, research on dual-polarization radar assimilation based on variational methods has also been carried out. Li et al. (2017) conducted case studies of single-station dual-polarization radar DA using variational methods, and their results showed that additional assimilation of differential reflectivity and specific differential phase could further improve reflectivity analyses and forecasts. Variational methods require the construction of tangent linear and adjoint operators for dual-polarization observations. To establish more reasonable tangent linear and adjoint operators, Kawabata et al. (2018) constructed the tangent linear and adjoint operators for dual-polarization observations based on liquid-phase particles, and Wang et al. (2019) developed the tangent linear and adjoint operators for horizontal/vertical reflectivity, including ice-phase particles.

3.3   Other types of observations

The radio occultation technique is regarded as one of the most promising means in current atmospheric detection. It can provide information on the neutral atmosphere and ionosphere with global distribution over the whole sky. In particular, the “Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC)” implemented in 2006 created conditions for the operational application of occultation DA. For the assimilation of occultation data, the most suitable assimilation quantities are the bending angle and refractivity and there are two kinds of forward operator, namely one-dimensional or two-dimensional. Currently, the advanced global operational centers nearly all have assimilated occultation data, indicating the importance of occultation data for NWP (Healy et al., 2005; Poli et al., 2009; Liu and Xue, 2014). The CMA-GFS model started to operationally assimilate occultation data in 2009. The assimilated quantity is the refractivity with height range from 1 to 50 km. The CMA-GFS achieved the assimilation of FY-3D GNOS occultation refractivity relatively early (Wang et al., 2020). Occultation data have always played an important role in the numerical prediction system of CMA, with the contribution to the 24-h forecast error always ranked among the top three in CMA-GFS.

Atmospheric motion vectors (AMV) are one of the satellite data that were among the earliest usage in DA. Even with a large amount of satellite data assimilated by the DA system, the role of AMV in improving NWP remains non-negligible (Forsythe et al., 2007). In the past decade, satellite wind retrieval algorithms have achieved remarkable progress in aspects such as the selection of tracer representative pixels and height assignment (Xu, 2020). In particular, due to the improvements of the selection of motion representative pixels and the estimation of translucent cloud height for FY satellite cloud motion winds (Zhang et al., 2017a, b), the cloud motion wind data of FY-2E had reached the level of similar international products as early as 2011 (Salonen and Bormann, 2015). Compared to other observation means, the errors of AMV are still relatively large. In particular, the error in specifying the height of clouds makes the height of the atmospheric motion represented by AMV highly uncertain. The thousands of channels in the vertical of the FY-4A GIIRS have brought new opportunities for the height assignment of the three-dimensional wind field. Studies have shown that the three-dimensional horizontal wind field can be effectively retrieved in clear-sky and partially cloudy areas (Ma et al., 2021; Li et al., 2022b), and reasonable assimilation of three-dimensional dynamic information has a positive effect on the track and intensity forecasts of typhoons (Meng et al., 2024).

A scatterometer is a spaceborne radar that measures the backscatter of the sea surface from multiple directions, from which the wind direction and wind speed of the sea surface can be derived (Stoffelen and Anderson, 1997). Satellites measure a set of backscatter values in different directions for the same sea surface area, but these can correspond to multiple different wind directions, thus bringing new problems to DA. Currently, when assimilating the scatterometer data, one wind direction closest to the background is generally selected from several possible wind directions. There are also cases where the wind speeds of multiple different wind directions are assimilated simultaneously, allowing the assimilation system to provide adaptive weights. In 2022, the assimilation of the scatterometer wind data of China’s HY-2B satellite was achieved in the CMA-GFS 4DVar, significantly improving the analyses in the lower troposphere over the ocean surface (Wang et al., 2023).

In 2018, the European Space Agency (ESA) successfully launched the world’s first spaceborne wind lidar satellite, ADM-Aeolus. It can provide high spatial and temporal resolutions, near-real-time global radial wind speed information with vertical resolutions ranging from 0.25 to 2 km from the ground to 30 km. In 2022, the operational application of Aeolus data was achieved for the first time in CMA-GFS. The assimilation of Aeolus data can significantly reduce the wind analyses errors in the tropics and the southern hemisphere. In the tropics, the reduction in the average error (compared with ERA5) can reach 10%. The forecast improvements for the first three days in the northern hemisphere, southern hemisphere, and tropics are relatively significant, and the prediction contribution in the east Asian region is neutral.

Additionally, China’s ground-based automatic weather station network has developed since 2008, reaching approximately 70,000 stations by 2023. High spatiotemporal resolution ground-based observations have become one of the key observation types of CMA-MESO. However, due to China’s complex terrain, including the “Roof of the World” with an average elevation above 4,000 m, plateaus such as the Yunnan–Guizhou Plateau, the Loess Plateau, and the Sichuan Basin with average elevations between 1,000 and 2,000 m, as well as hills and plains below 1,000 m, there are significant height differences between the relatively smooth model terrain and the actual terrain of observation stations. Studies by Xu Z. F. et al. (2006, 2007, 2009) have shown that failing to effectively resolve the height differences between the model and stations can negatively impact the assimilation of surface data. Therefore, Xu et al. (2021, 2023) developed DA schemes for 2-m relative humidity and temperature observations in complex terrain within the CMA-MESO 3DVar. These schemes improved both the quantity and quality of assimilated surface data and enhanced forecasts of surface state variables. Moreover, a DA scheme for the surface pressure based on the hydrostatic equation was implemented to replace the surface pressure extrapolation assimilation scheme by Lian and Xue (2010). This new approach not only increased the utilization of surface observations but also resolved the issue of deteriorating precipitation forecasts with increased number of ground observations.

4.   The target observation strategies

Large uncertainties often occur in high-impact weather forecasts, such as heavy rainfalls, typhoons, etc. Not enough data with either quantity or quality is an important obstacle limits the forecast skill for high-impact weather events. Hence, Snyder (1996) proposed the concept of target observation strategy, which adds a few observations with high quality in sensitive areas to improve the forecasts in concern. Extensive researches have demonstrated that additional observations in sensitive areas are greatly helpful for accurate forecasts of high-impact weather events such as typhoons (Aderson, 2010). For example, the Observing System Research and Predictability Experiment (THORPEX) displayed the importance of target observations in improving the track forecasts of typhoons (Shapiro and Thorpe, 2004). In recent years, China has made significant advancements in the research and application of targeted observations for forecasting high-impact weather events. (Duan et al., 2023).

4.1   Concept of target observations

THORPEX, under the auspices of WMO, is a world weather research programme accelerating improvements in the accuracy of one day to two week high-impact weather forecasts for the benefit of society, economy and environment. As an important part of THORPEX, target observations were highly promoted and refer to the augmentation of the regular observing networks with additional observations in some sensitive but data-sparse areas, aiming to reduce the initial condition errors and improve the subsequent forecasts.

The field campaigns of target observations have advanced rapidly under THORPEX. The Atlantic-THORPEX Regional Campaign (A-TReC) started during the autumn of 2003 for the northern hemisphere. A large quantity of in-situ and remotely sensed observations was collected, targeted at 1- to 3-d forecasts of potential high-impact weather events over Europe (Rabier et al., 2008). The African Monsoon Multidisciplinary Analysis campaign, utilizing the rawinsonde and driftsonde balloons, was aimed at improving short-range forecasts of western African rainfall and easterly waves that may lead to tropical cyclogenesis (Agustí-Panareda et al., 2010). Several campaigns over Europe, aimed partially at improving short-range forecasts of specific high-impact weather events such as winter flow distortion past Greenland, summer rainfall in Central Europe or autumn heavy precipitation events in the Mediterranean region have taken place between 2007−2009 (Wulfmeyer et al., 2008；Prates et al., 2009；Jansa et al., 2011). The THORPEX Pacific Asian Regional Campaign (T-PARC) possessed a broader scope than aforementioned experiments and focused on the large northern Pacific basin. The summer phase in 2008 was aimed at investigating a wide variety of issues related to the science and predictability of the life cycle of typhoons from formation through recurvature and extratropical transition, including the impact on the flow far downstream in the mid-latitude storm track (Elsberry and Harr, 2008). In the winter phase, the primary purpose was to investigate the potential for targeted aircraft and rawinsonde observations to improve forecasts of weather systems over north America beyond the 1- to 3-d ranges. In order to improve forecasts over Scandinavia and the use of data from polar-orbiting satellites over the Antarctic, field campaigns of target observations also covered the polar areas (IPY; Irvine et al., 2011).

4.2   Target observation strategies

The target observation strategy requires additional observations in some key sensitive areas to reduce the initial condition errors and further improve the forecast skills of high-impact weather events. The methods to identify the sensitive areas can be classified into two types, one is based on the analysis sensitivity, and the other is based on the observation sensitivity.

Analysis sensitivity captures the dependence of forecast uncertainty on initial perturbations at different sites. The higher the sensitivity, the larger forecast errors are induced by the initial errors at these sites. Hence, these sites are identified as the sensitive areas for target observations. The typical analysis sensitivity methods include the adjoint sensitivity (Bergot, 1999; Wu et al., 2007), singular vectors (SVs; Palmer et al., 1998), ensemble transform technique (Bishop and Toth, 1999), etc. The former two methods have been widely utilized in the field campaigns of target observations as T-PARC and IPY.

Observation sensitivity introduces observations and DA techniques, which identifies the sensitive areas by evaluating the reduction of the forecast errors led by the simulated observations at different sites. The sites where the observations can bring the maximum reduction of forecast errors are identified as the sensitive areas. The typical observation sensitivity methods include the Hessian SVs (Barkmeijer et al., 1998) and ensemble transform Kalman filter (ETKF; Bishop et al., 2001). The latter has been utilized in many field campaigns as A-TReC, T-PARC, and IPY. Results from the T-PARC summer phase showed improved track forecasts of two long-lived typhoons, with 20%–40% error reduction in numerical models of NCEP and Korea Meteorological Administration, although little track forecast improvements in the models of ECMWF and Japan Meteorological Agency with forecast lead times longer than 72 h (Weissmann et al., 2011).

All the methods mentioned above utilize linear approximation to some degrees; however, the atmospheric state and its evolutions are characterized by nonlinear nature. Mu et al. (2009) proposed to identify the sensitive areas of target observations according to the structures and locations of initial perturbations through the CNOP method (Mu et al., 2003). Plenty of observation system simulated experiments demonstrated that assimilating observations in the sensitive areas identified by the CNOP method can improve the track forecasts of typhoons, which is more prominent than that by traditional SVs (Qin and Mu, 2012; Chen et al., 2013; Feng et al., 2022; Chan et al., 2023; Qin et al., 2023). In recent years, the CNOP method has also been used in theoretical research and field campaigns of offshore ocean environment forecasts, which greatly improved the ocean forecast skills (Liu et al., 2021). As an effective method to identify the sensitive areas of target observations for both atmosphere and ocean (Mu et al., 2017; Duan et al., 2018; Jiang et al., 2022, 2024; Yang et al., 2022, 2023), the CNOP method is expected to be further applied in real-time operational forecasts and improve the NWP forecasts (Fig. 2).

Download

Larger image

PowerPoint slide

Fig  2.  The CNOP method and its application in forecasts of extreme weather and climate events.

4.3   Practices and experiences in China

The national Landfalling Tropical Cyclone Research Project (LTCRP) in China was conceived and funded in 2009. The main objectives of the project are to investigate the characteristics of structure and intensity changes during TC landfall and associated physical mechanisms, and then to develop new technology to advance the forecast skills for landfalling TCs (Duan et al., 2019). The national LTCRP has been ongoing for about 10 years from 2008 to 2018 and includes three main components: field campaigns, scientific research, and technical developments. Field campaigns for 24 TCs under the national LTCRP have obtained a large amount of observations, which helps understanding the boundary structure of TCs (Zhang J. A. et al., 2011, 2015; Ming et al., 2014; Bi et al., 2015; Tang et al., 2015; Zhao et al., 2015; Wang M. J. et al., 2016, 2018; Zhao et al., 2017; Ming and Zhang, 2018; Wen et al., 2018; Wu et al., 2018) and promotes the developments of new DA technique, nowcast system, and assessment system (Cha and Wang, 2013; Wen et al., 2017; Li et al., 2018; Liu et al., 2018; Lu et al., 2018; Chen et al., 2019; Bao et al., 2023).

The operational GRAPES model achieved identifying sensitive areas using the SVs in 2013 (Liu et al., 2013; Li and Liu, 2019) and started to conduct field campaigns of target observations for TC forecasts. The forecasts in GRPAES-SVs usually concern a wide area (10°–35°N, 105°–125°E) of southern China and near sea (Liu et al., 2019, Zhang et al., 2019). A few years later, several operational and scientific research departments have started cooperation in 2020 to utilize the CNOP method to identify the sensitive areas of target observations for TCs (Duan and Qin, 2022; Duan et al., 2023). Field campaigns were conducted for TCs Higos (2020), Maysak (2020), Chanhom (2020), Conson (2021), Chanthu (2021), and Mulan (2022), and the forecasts have been improved (Feng et al., 2022; Chan et al., 2023; Qin et al., 2023).

Field campaigns of target observations for TCs have been going on after the national LTCRP. The Hong Kong Observatory has conducted field campaigns for the TCs over South China Sea using dropsondes since 2006 (Chan et al., 2018). In 2020, Meteorological Observation Centre of CMA organized a field campaign for TC Sinlaku (2020) using an unmanned aerial vehicle. The collected data help understand the formation and effect of helicall roll in the TC boundary layer (Chen et al., 2021; Tang et al., 2021) and improve the forecast skills of TC track and intensity. Taking the TC Atsani (2020) for example, assimilating several dropsonde data in the sensitive areas identified by the CNOP method has obtained comparable track forecasts as assimilating all available dropsonde data.

The successful launched FY-4A in December 2006 supplements a new and powerful observing technique for field campaigns of target observations (Han et al., 2023). In 2018−2021, FY-4A conducted nine field campaigns of target observations, of which eight were aimed at TCs (Lei et al., 2019; Meng D. M. et al., 2019; Han et al.,2025). Taking the TC Maria (2018) for example, FY-4A scanned it at every 15-m. Such high frequent data help profile the atmospheric structure around the TC (Yin et al., 2021). Feng et al. (2022) showed that the data obtained by FY-4A improved the TC track forecasts with lead time longer than 2.5 d, especially with a 50% reduction of track forecast errors in 3- to 3.5-d forecasts. In particular, assimilation of the FY-4A data corrected the wrong landfall forecast of TC Chanthu in Taiwan Island.

FY-4B was successfully launched in June 2021 and rapidly conducted target observations three times from 2021 to 2023, including two for TCs. In 2022, assimilating the FY-4B data improved forecast skills for TC Mulan (2022) in track, intensity, and precipitation (Chan et al., 2023). Without the FY-4B target observations, TC Mulan (2022) was forecasted to make landfall in Leizhou peninsula. However, TC Mulan was forecasted to move westward and then turn northward with the FY-4B target observations assimilated, which is much closer to the fact.

In summary, China has made breakthrough in target observations for TCs, promoted interactions between forecasts and observations, fulfilled a combination of nonlinear method in identifying sensitive areas and field campaigns of target observations. All these achievements have provided theories and techniques for field campaigns of target observations for TCs and other high-impact weather event forecasts.

5.   China’s operational DA systems for NWP

5.1   Development of operational DA systems for NWP

DA systems built on theoretical and methodological foundations, are a critical component of operational NWP systems. Operational DA systems typically include real-time observation acquisition modules, observation pre-processing modules, and assimilation modules, along with the core methods of observation quality control and DA. Since the 1980s, as China’s operational NWP systems have advanced, the DA systems have evolved from simple to complex, from primarily imported systems to independently developed ones (Fig. 3).

Download

Larger image

PowerPoint slide

Fig  3.  Journey of the independently developed NWP operational DA systems at CMA.

China’s earliest operational NWP models were the three-layer primitive equation model (Model A), starting in July 1980, and the five-layer northern hemisphere primitive equation model (Model B), operational since February 1982. At that time, assimilation, also called objective analysis, employed the relatively simple SCM method, primarily assimilating conventional observations (Wang et al., 1984).

From the late 1980s, the National Meteorological Centre of China gradually established and developed the limited area forecast system model (LAFS), operational in 1991, and the high-resolution LAFS model (HLAFS), operational in May 1996, and meanwhile regional DA was also achieved. The assimilation scheme was an improved version of OI from the U.S. National Meteorological Center, performing 3D multivariate analysis for height and wind fields and univariate analysis for relative humidity, along with nonlinear normal mode initialization (Xue et al., 1992), to assimilate conventional and non-conventional observations, including the satellite-derived moisture, temperature profiles, and cloud-tracked winds (Guo et al., 1995). During the “Ninth Five-Year Plan” period, the National Meteorological Centre adopted the Mesoscale Model 5 (MM5) from the National Center for Atmospheric Research (NCAR) and updated the SCM method to a dynamical relaxation method for assimilating conventional observations (Jiao, 2010).

Later on, CMA imported the ECMWF spectral model to establish a global medium-range NWP operational system, which included the T42L9 system (operational in June 1991), the T63L16 system (operational on the domestically developed Galaxy-II supercomputer in October 1993), the T106L19 system (operational on CRAY C92 in July 1997), and the T213L31 system (operational in September 2002) (Jiao, 2010). In these T-series global medium-range forecasting systems, a global DA system based on OI was established, with nonlinear normal mode initialization to ensure more coherent alignment between the analyses and forecasts. Assimilated observations included the weather reports received via domestic communication lines and GTS (Li and Qiu, 1992; Li, 1994). The T639L60 system, operational right before the 2009 flood season, introduced NCEP’s GSI variational assimilation system, marking a transition from OI to 3DVar, with the ability of assimilating microwave sounding data from polar-orbiting satellites (Guan et al., 2008).

By the late 20th and early 21st centuries, CMA shifted its strategy from importing systems to independently developing its own systems. This decision led to the development of the first-generation multi-scale and unified DA and NWP system—GRAPES (Xue and Chen, 2008; Shen et al., 2020). In July 2006, the GRAPES regional model system (GRAPES-MESO V2.0) with a 30-km horizontal resolution became operational, featuring a regional isobaric 3DVar assimilating conventional observations, with twice-daily cold starts. Before the 2008 flood season, GRAPES-MESO was upgraded to V2.5 with a 15-km resolution and isobaric 3DVar system. In 2010, the GRAPES-RAFS (Rapid Analysis and Forecast System) began quasi-operational use, offering 3-h cycles of analyses and forecasts (Xu et al., 2013). In July 2014, GRAPES-MESO V4.0 was launched with a 10-km resolution, a new cloud analysis module, and the capability to assimilate GPS/PW, FY-2E cloud drift winds, and GNSS/RO data (Huang et al., 2017; Shen et al., 2020). In June 2020, CMA-MESO V5.0 became operational, upgrading the DA system to a 3DVar with 3-km resolution and 3-h cycles of analyses and forecasts (Huang et al., 2022).

The first formal version of the GRAPES global model, GRAPES-GFS V1.0, began quasi-operational use in 2009, with initial conditions provided by a global isobaric 3DVar system. In 2016, GRAPES-GFS V2.0 with a 25-km horizontal resolution became operational, upgrading to a model-level 3DVar to reduce errors introduced by spatial interpolation and variable transformation along with improved background error covariances (Wang J. C. et al., 2014). DA for the satellite and occultation data were also enhanced (Wang J. C. et al., 2015, 2016; Han and Bormman, 2016). GRAPES-GFS V2.2, operational in July 2018, upgraded the assimilation system from 3DVar to 4DVar (Zhang et al., 2019). This marked China’s entry into the international forefront of operational 4DVar systems, making CMA one of the few operational centers with self-developed and operational 4DVar systems (Shen et al., 2020). In 2021, GRAPES-GFS and GRAPES-MESO were renamed CMA-GFS and CMA-MESO. In May 2023, CMA-GFS V4.0 became operational, increasing the global 4DVar system’s resolution from 25 to 12.5 km and replacing the RTTOV radiative transfer model with the domestic ARMS model, a significant step toward independent control of core technologies in operational DA systems. According to the WMO model verification standards, forecast improvements can be measured by the historical evolution of the 500-hPa geopotential height anomaly correlation coefficient (ACC) in the northern and southern hemispheres (Fig. 1). For the northern hemisphere, the 5-d ACC at 500 hPa improved from 0.695 in September 2010 to 0.768 in September 2015, further reaching 0.846 in September 2023.

In June 2024, CMA-MESO V6.0 with nationwide 1-km resolution passed operational evaluation and began real-time parallel operational trials, upgrading the 3DVar system with 1-km horizontal resolution and 1-h cycle of analyses and forecasts. Meanwhile, CMA-GFS V4.2 upgraded the global 4DVar system to a global En4DVar system, with a year-long retrospective test showing systematically better analyses and forecasts than the global 4DVar operational system, now the global En4DVar has been operationally implemented since 31 December 2024.

5.2   Global 4DVar operational system

4DVar is an extension of 3DVar in the time dimension, and thus 3DVar framework is the foundation for the 4DVar framework. The CMA-GFS global 3DVar and 4DVar systems use stream function, unbalanced velocity potential, unbalanced nondimensional pressure, and specific humidity as the analysis variables. The balanced components of velocity potential and nondimensional pressure are calculated using a combination of dynamical and statistical methods. Since the analysis variables are independent, the corresponding background error covariance matrix is diagonal. A second-order autoregressive correlation function is used for the horizontal correlation of univariate variable, computed via spectral filtering, while the lengthscales of horizontal and vertical correlations are statistically derived from ensemble samples.

Additionally, the CMA-GFS global 3DVar and 4DVar systems employ an incremental scheme. In this incremental scheme, the high-resolution forecast model is integrated to calculate the observation increments, while low-resolution tangent linear and adjoint models are used for the minimization process, significantly reducing the computational cost and improving the efficiency of functional minimization.

The global tangent linear model and adjoint model are the core components of the global 4DVar system. The CMA-GFS global tangent linear and adjoint models are the first non-hydrostatic ones applied in an operational 4DVar system internationally. Moreover, the computational time for the dynamical frameworks of the tangent linear and adjoint models is only about three times that of the dynamical framework of the global forecast model, demonstrating outstanding computational performance. The tangent linear and adjoint models also incorporate comprehensively linearized physical processes, including vertical diffusion, subgrid-scale topography parameterization, large-scale condensation, and convective parameterization (Gong et al., 2019; Liu et al., 2019).

The CMA-GFS global 4DVar operational system developed preconditioned Lanczos-CG and L-BFGS algorithms. By default, the Lanczos-CG algorithm is used for its faster convergence and smoother process. However, the L-BFGS algorithm has better fault tolerance, and the system automatically switches to L-BFGS when the Lanczos-CG algorithm fails to converge.

The basic configuration of the CMA-GFS global 4DVar operational system (Version 4.0) includes a horizontal resolution of 0.125°/0.75° (outer loop/inner loop), 87 vertical layers, model integration time step of 300 s/900 s (outer loop/inner loop), a 6-h assimilation window, observation profiling interval of 30 m, and a maximum of 50 iterations for minimization. To balance the tradeoff between the qualify of analyses and forecasts and timeliness, the CMA-GFS global 4DVar operational system runs one analysis–forecast cycling system and one analysis–forecast system. The CMA-GFS global 4DVar operational system is performed four times daily, producing model initial conditions for standard time points that are then used for subsequent10-d forecasts.

The CMA-GFS global assimilation system has developed critical technologies for satellite DA, including quality control, cloud detection, and bias correction. It has also established a domestically developed fast radiative transfer model ARMS replacing the previously used RTTOV. CMA-GFS global 4DVar has successively incorporated satellite data from FY polar-orbiting microwave temperature and humidity sounders (Xiao et al., 2023a, b), microwave imagers (Xiao et al., 2020), occultation data (Wang et al., 2020), FY geostationary infrared hyperspectral (Yin et al., 2020, 2021; Han et al., 2023), infrared imagers (Wang and Han, 2018), HY-2B microwave imager SMR (Li and Han, 2024), scatterometer ocean winds (Wang J. C. et al., 2023), etc. Currently, the CMA-GFS global 4DVar operational system assimilates observations including radiosonde, surface, aircraft reports, cloud drift winds, occultation data, scatterometer winds, GPS precipitable water, GNSS reflectometry winds, NOAA, METOP, FY-3 microwave temperature and humidity sounders, infrared hyperspectral data, GCOM microwave imagers, and FY2 imagers. Satellite observations accounts for approximately 80% of all observations, playing a critical role in global DA.

5.3   The regional 3DVar operational system

The regional numerical forecasting system implemented in CMA is CMA-MESO (formerly GRAPES-MESO), which is mainly used to improve the prediction of hazardous weather and lower tropospheric phenomena. In order to rapidly update the model trajectories, CMA-MESO uses a 3DVar assimilation system and a cloud analysis system to assimilate spatiotemporally dense observations. In June 2020, the CMA-MESO v5.0 system with 3-km horizontal resolution and 3-h update was operationally implemented (Shen et al., 2020; Huang et al., 2022). Since October 2024, the CMA-MESO v6.0 system with 1-km horizontal resolution and 1-h update has been operationally implemented.

The CMA-MESO kilometer-scale 3DVar is based on the GRAPES unified variational data assimilation framework, and has been further developed for convective-scale weather systems. In terms of the analysis framework, new minimization control variables are constructed based on the dynamical characteristics of the convective-scale system, and a simplified weak constraint of the continuum equation is introduced to get better balanced small- and medium- scale analyses (Wang et al., 2024). Meanwhile, multi-scale analysis schemes are also developed to capture the large-scale patterns, which include a horizontal correlation model using multiple Gaussian scale superposition, a weak constrain with large-scale information, and a blending scheme that merges small- and medium- scale information to the global large-scale information (Zhuang et al., 2020; Wang R. C. et al., 2021). The kilometer-scale 3DVar mainly updates wind, temperature, pressure, and humidity variables, while the cloud fields are diagnostically updated by the cloud analysis system and introduced into the model trajectory by nudging (Zhu et al., 2017). Based on the operational 3DVar system, a prototype En3DVar system has also been developed to obtain flow-dependent analyses. In addition, cloud control variables have also been added into the En3DVar framework to provide a basis for direct assimilation of radar reflectivity.

In terms of the application of observations, CMA-MESO kilometer-scale 3DVar takes full advantages of shared observation operators of the unified variational framework and achieves the direct assimilation of conventional and satellite observations. On this basis, the CMA-MESO system further develops assimilation algorithms for spatiotemporally dense observations. For radar data, radial wind and wind profile radar data are directly assimilated in the kilometer-scale 3DVar, and reflectivity data is assimilated in the cloud analysis system. The observations of China’s new generation of geostationary satellites, FY4A and FY4B, can be directly assimilated in the kilometer-scale 3DVar, while the FY2G brightness temperature and total cloud cover data are assimilated in the cloud analysis system. For automatic surface observations, kilometer-scale 3DVar can assimilate multiple observations such as 10-m wind, 2-m temperature, 2-m humidity, and surface pressure. In addition to conventional surface variables, GNSS/MET humility data is also assimilated in the kilometer-scale 3DVar. Currently, a total of 17 types of observations are assimilated in the CMA-MESO system, with radar data contributed to the highest percentage.

The assimilation of radar reflectivity data can improve the TS scores of heavy precipitation forecasts by 5%–18%, and the assimilation of radial winds further improves the quality of low-level winds. The assimilation of wind profile radar observations can also improve the track and intensity forecasts of typhoons (Wang et al., 2019). Since the cloud analysis system relies on empirical relations and has limitations for the convective-scale numerical simulations, the development of direct assimilation techniques for radar reflectivity (including dual polarization quantities) is being carried out based on the CMA-MESO 1-km variational assimilation system. In addition, research on the assimilation of new types of ground-based remote sensing observations, including the X-band weather radar data, microwave radiometer temperature and humidity profile products, and cloud radar products, is also under way.

5.4   DA of TC observations

Since 2004, a set of TC initialization schemes (Qu et al., 2009), which consists of initial vortex formation, vortex relocation, and vortex adjustment, has been successively developed and adopted in operation on T213 global spectral model of CMA for TC forecasts. In 2014, the bogus vortex embedding in the initial vortex formation scheme was replaced by DA of vortex wind and pressure data, which were developed and applied to T639 Global Spectral Model (Qu et al., 2016). In 2018, with the continuous improvement of CMA’s self-developed global modeling system (CMA-GFS), a new initialization scheme that assimilates the evolution trend of the TC position and central pressure profile based on the 4DVar assimilation system, was produced (Qu et al., 2022). The operational application shows that the new initialization scheme can significantly improve the TC track and intensity forecasts of CMA-GFS. There are two prospects for future advancements of improving the analyses and forecasts of TCs, one to develop TC initial perturbation techniques, aiming to improve the flow-dependent background error statistics for TCs, and the other to develop En4DVar with effective and efficient assimilation of observation information for TCs.

6.   Prospects for DA

6.1   Integrating DA for high resolutions and long lead times

One development aspect for NWP is persistently improving the model resolution, which leads to better resolved small-scale, nonlinear, and non-Gaussian processes. Meanwhile, advances in observing systems produce more indirect types of observations and increased observations with higher spatiotemporal resolutions (Fig. 4). Consequently, the current mainstream DA methods that often assume Gaussian error distributions and linear relationships may no longer be optimal (Yano et al., 2018). The iterative EnKF (IEnKF; Sakov et al., 2012) replaces linear regression with a transform matrix from the previous iteration, better capturing the nonlinear error growth. Other nonlinear filtering methods have been proposed, including the Gaussian mixture filters (Bengtsson et al., 2003), maximum likelihood ensemble filters (Zupanski, 2005), rank histogram filters (Anderson, 2010), rank-matching filters (Lei and Bickel, 2011), quantile-conserving ensemble filters (Anderson, 2022, 2023).

Download

Larger image

PowerPoint slide

Fig  4.  Prospect for the development of DA in NWP.

A fully Bayesian, nonlinear method is the particle filter that provides particles following the Bayesian posterior distribution. Particle filter can be implemented through methods like bootstrap sampling (Gordon et al., 1993; Douc and Cappé, 2005), importance sampling (van Leeuwen, 2003; Robert and Cassela, 2004), and importance sampling with proposal (Doucet et al., 2000; Spiller et al., 2008). However, particle filter faces the “curse of dimensionality” when it applied in high-dimensional dynamical systems (Snyder et al., 2008). To mitigate the “curse of dimensionality”, techniques such as the implicit particle filter combining filtering and smoothing (Chorin and Tu, 2009), hybrid ensemble-particle filters (Santitissadeekorn and Jones, 2015), equal-weight particle filters (Ades and van Leeuwen, 2015; Zhu et al., 2016), and localized particle filters (Penny and Miyoshi, 2016; Poterjoy, 2016) have been developed.

Another trend in NWP is extending the forecast lead times, from the current forecast limit of two weeks to seasonal-to-decadal (s2d) lead times (Fig. 4). To extend the predictability, coupling atmospheric models with slowly evolving components of the Earth system, such as the ocean, land surface, and cryosphere, becomes necessary. Thus, the advance of coupled DA becomes essential. Weakly-coupled DA assimilates component-specific observations within each model component, using the coupled model to spread the observation information cross components (Zhang et al., 2007; Sugiura et al., 2008; Saha et al., 2010; Laloyaux et al., 2016). Comparatively, strongly-coupled DA uses cross-components observations to simultaneously update state variables from all components (Sluka et al., 2016; Sun et al., 2020a), providing more balanced analyses and reducing initialization shock, and strongly-coupled DA can lead to improved forecasts (Smith et al., 2015; Chen and Zhang, 2019).

6.2   Hybrid machine learning and DA

Machine learning (ML) has brought new opportunities for DA, numerical models, predictions and projections, especially in efficiently and effectively assimilating the vast satellite observations. Both DA and ML can be viewed as inverse problems under the Bayes theory (Geer, 2021). Compared to traditional DA methods, ML has advantages for capturing nonlinear features, approximating nonlinear systems, and processing massive observation datasets, and thus, hybrid ML and DA approaches have been rapidly developed (Fig. 4). Direct integration of ML and DA can use data-driven ML models to replace numerical models, providing efficient forecasts for DA cycles. FengWu-4DVar combines the multi-modal neural network-based FengWu model with 4DVar, to leverage the short-term forecasts and automatic differentiation capability of DL for efficiently solving the 4DVar analysis (Xiao et al., 2024a). Similarly, the ClimaX based on the Vision Transformer (ViT) architecture, can integrate with LETKF, enabling cyclical ensemble DA while diagnosing ML-based models (Kotsuki et al., 2024).

ML can also enhance various DA components. Observation forward operators convert model state variables to estimated observations, but the Jacobians of forward operators could be computationally expensive. ML can efficiently approximate the forward operators, their Jacobians, and even the second-order Hessian matrices, especially for highly nonlinear forward operators (Storto et al., 2021). ML-based forward operator for satellite radiances can replace the fast radiative transfer model and simultaneously perform bias correction (Liang et al., 2023). The adjoint models required by 4DVAR is difficult to construct and also computationally expensive. ML offers an alternative by simulating the physical parameterizations to directly obtain the tangent linear and adjoint models, significantly improving the computational efficiency (Hatfield et al., 2021). Covariance localization is essential for EnKF successfully applied in high-dimensional dynamical systems, but the widely used localization functions are symmetric functions of distances. ML can extract nonlinear error characteristics from data, and generate non-symmetric and nonlinear localization functions (Wang Y. M., et al., 2021).

As a major source of forecast errors, model errors need be appropriately addressed in DA. ML can learn model errors resulting from unresolved processes in numerical simulations (Rasp et al., 2018; Bolton and Zanna, 2019; Gagne II et al., 2020; Brajard et al., 2021). Neural networks and other ML architectures can learn model errors with nonlinear and multi-scale characteristics from the analyses, forecasts, and observations, which can then be represented in DA through additive terms (Bonavita and Laloyaux, 2020; Farchi et al., 2021b). Moreover, through cycling DA, ML can learn and correct model errors online using priors and posteriors (Farchi et al., 2021a; Peng et al., 2024), or simultaneously estimate errors in state variables and model parameters (Bocquet et al., 2021; Malartic et al., 2022).

To address the computational complexity of high-dimensional dynamical systems, latent-space DA has been proposed, which assimilates data in the latent space created by ML-based autoencoders, combining ML’s efficiency with DA’s optimization (Binev et al., 2017; Arcucci et al., 2019; Casas et al., 2020). Compared to the variational methods and EnKFs that often assume Gaussian error distributions, variational autoencoders can estimate non-Gaussian error distributions, and the combined variational autoencoders and variational methods outperform the traditional variational methods (Xiao et al., 2024b). ML can also integrate with nonlinear filters, like the deep Kalman filters (Krishnan et al., 2015, 2017) and Kalman variational autoencoders (Fraccaro et al., 2017). DiffDA, based on the GraphCast neural network, leverages similarities between numerical models and denoising diffusion models to directly produce the analysis (Huang et al., 2024). Furthermore, ML can facilitate reconstruction and assimilation based on incomplete or coarse-resolution observing networks (Wang et al., 2022; Howard et al., 2024).

6.3   DA for all-sky satellite radiances and dual-polarization radar data

In China’s operational systems, the assimilated satellite data primarily focuses on clear-sky radiances, with the assimilation of satellite data in rain and cloud areas has not yet been operationally implemented. The assimilation and application of new remote sensing data, including active remote sensing instruments such as precipitation radars and wind lidars, as well as ground-based payload observation data, still need to be improved. For Earth system coupled DA, the fast radiative transfer model needs to consider an ultra-high spectral atmospheric transmittance calculation model covering the full spectrum range and containing multiple atmospheric components. Additionally, constructing a satellite and new payload coupled forward operator based on artificial intelligence is essential. The assimilation of visible data will be a very important direction for the future application of satellite DA. Specifically, the all-sky assimilation of infrared and visible data will provide more accurate analyses for numerical models, particularly at the convective scale (Schröttle et al., 2020). Developing advanced radiative transfer models that account for the scattering characteristics of cloud particles will enable better assimilate of infrared radiances affected by clouds. Meanwhile, an improved description of land surface processes is worth to develop, aiming to increase the assimilation of radiance in surface-sensitive channels. Microwave sounding data contribute the most to the forecast accuracy of NWP. But currently, microwave sounders are only carried on low-earth orbit meteorological satellites, with long revisit cycles and low temporal resolutions, making it difficult to provide continuous observation for weather systems. The geostationary orbit microwave sounding being designed and constructed for China’s FY satellites is an important supplement to the existing microwave sounding system (Lu and Gu, 2016). It can not only provide high temporal resolution 3D information for the atmospheric thermodynamic structure but also generate wind field products at different heights over the whole sky using the water vapor tracking method, providing more dynamical information for NWP (Zhang et al., 2021). In addition, the current satellite DA techniques basically neglect the synergy effect among the multi-instrument observations. Fully considering the synergy and complementary effects among the multi-instrument observations, such as the joint application of imaging and sounding data (Di et al., 2024), can lead to improved assimilation outcomes. This represents a significant future development direction for satellite DA.

Developing a reasonable variational-based dual-polarimetric radar observation forward operator is crucial. Based on the dual-polarization radar observation operator developed by Kawabata et al. (2018), Zhang et al. (2024) developed a variational direct assimilation scheme for dual-polarization radar based on hydrometeor control variables. Cycling assimilation and forecast experiments for real cases demonstrated that the assimilation of dual-polarization radar data can improve the thermodynamic and microphysical characteristics of both the analyses and forecasts. Research on the development of dual-polarization radar observation forward operators and adjoint operators provides foundational support for the variational assimilation of polarimetric radar quantities. However, due to the complexity of the tangent linear and adjoint operators of polarimetric radar and the uncertainties in parameterization schemes, further research for more refined polarimetric radar observation forward operators is needed, particularly with respect to the treatment of ice-phase and mixed-phase hydrometeors.

6.4   Advanced operational DA systems

Currently, China’s independently developed global and regional operational weather forecasting systems, CMA-GFS and CMA-MESO, can effectively assimilate multi-source observations, and play an important role in daily weather forecasts, warning, and meteorological disaster prevention and mitigation. However, the existing model dynamical frameworks are insufficient to meet the demand for seamless NWP of the Earth system. Major operational centers have been developing quasi-uniform grid numerical models, and China has also completed the next-generation high-precision scalable atmospheric model (Li et al., 2020). Development of the high-precision scalable atmospheric DA system for the next-generation model is progressing intensively. The goal for the next-generation atmospheric DA system is to establish an integrated global/regional hybrid ensemble-variational assimilation framework, with the basis of global 4DVar. Research on advanced assimilation techniques for multi-source observations, including the Beidou navigation soundings, satellite data in cloud and precipitation areas, radar reflectivity, X-band radar, dual-polarization radar, phased-array radar, and ground-based vertical remote sensing, will be conducted to promote the operational use of innovative observations. Meanwhile, studies on applying artificial intelligence algorithms in various areas of multi-source DA are underway.

To address the multi-scale seamless NWP for the Earth system from weather to climate, research on DA techniques across different components of the Earth system will be conducted. Based on the high-precision scalable atmospheric DA system, the goal is to develop an ocean–land–atmosphere–ice coupled DA system to achieve effective and efficient assimilation of multi-source observations across different components of the Earth system, providing coherent, balanced, and high-quality initial conditions for the Earth system model. As illustrated in Fig. 4, the aim is to establish an operational DA system for the Earth system prediction and projection, based on the advances in DA theories and methods.