Research and Operational Development of Numerical Weather Prediction in China

Numerical weather prediction (NWP) is a core technology in weather forecast and disaster mitigation. China’s NWP research and operational applications have been attached great importance by the meteorological community. Fundamental achievements have been made in the theories, methods, and NWP model development in China, which are of certain international impacts. In this paper, the scientific and technological progress of NWP in China since 1949 is summarized. The current status and recent progress of the domestically developed NWP system—GRAPES (Global/Regional Assimilation and PrEdiction System) are presented. Through independent research and development in the past 10 years, the operational GRAPES system has been established, which includes both regional and global deterministic and ensemble prediction models, with resolutions of 3–10 km for regional and 25–50 km for global forecasts. Major improvements include establishment of a new non-hydrostatic dynamic core, setup of four-dimensional variational data assimilation, and development of associated satellite application. As members of the GRAPES system, prediction models for atmospheric chemistry and air pollution, tropical cyclones, and ocean waves have also been developed and put into operational use. The GRAPES system has been an important milestone in NWP science and technology in China.


History of NWP research and application in China
The numerical weather prediction (NWP) objectively and quantitatively calculates future weather based on mathematical and physical laws. Following almost 70 years of development since its experimental success in the 1950s, NWP has become a sophisticated, strict, interdisciplinary science and technology. By virtue of NWP, weather forecast turns into an objective, quantitative science from the statistics and weather-map based, semiempirical technology (Zeng, 2013;Benjamin et al., 2019). The history of NWP, from its proposal a century ago to its 70 years of development and application, demonstrates that progress in NWP rides on a multi-year, steady, continuous accumulation of scientific knowledge and technological advances, and is thought to be one of the most influential achievements among the various fields of physics (Ji, 2011;Bauer et al., 2015). It can be argued that NWP comprehensively reflects the level of a country's meteorological technology, and is an important component of the core technology that ensures the well-being of its economy and people's livelihoods as well as its defense security.
China started the NWP research in 1954, and is one of the countries initiating such research relatively early. From the 1950s to the late 1960s, the NWP research was highly active and made advances that had considerable international influence (Gu, 1959;Blumen and Washington, 1973). Most of the research in the 1950s focused on quasi-geostrophic barotropic models. In an article published in 1959, Gu presented a systematic summary. He noted that the essentials of NWP were our understanding to the essence of atmospheric processes as well as the roles of physical factors. He also pointed out that extensive work needed to be done to develop methods for establishing, selecting and solving models (Gu, 1959). In the early research, two-layer quasi-geostrophic equations were solved by using graphical methods. This approach produced relatively good short-range forecasts for outbreak of cold surge and 48-h forecasts of the 500-hPa geopotential height (Liao, 1956;Gu et al., 1957). It is worth noting that Chen et al. (1957) studied forecast of front using a two-layer model with potential temperature as the vertical coordinate. This is believed to be a relatively early study that employed potential temperature as the vertical coordinate. In February 1960, the Numerical Prediction Team (NPT) of the Central Meteorological Bureau of China (CMBC) collaborated with the Institute of Geology and Geophysics and the Institute of Computing Technology of the Chinese Academy of Sciences (CAS) to publish 24-and 48-h 500-hPa geopotential height forecasts for Asia and Europe using a quasi-geostrophic barotropic model (Numerical Prediction Team of the Central Meteorological Bureau of China, 1965). They presented an excellent summary in a paper titled "China's numerical prediction operations," in which they noted the importance of difference scheme design, initial wind field, and topographic effect, and further pointed out that baroclinic models should be the future direction (Numerical Prediction Team of the Central Meteorological Bureau of China, 1965). Shen and Mu (1965) systematically summarized the application of 48-h 500-hPa geopotential height maps in daily weather forecast during 1964-1965. They emphasized the importance of understanding NWP errors and its capabilities in practical applications. To this day, these conclusions are still extremely important for guiding the use of NWP.
While studying quasi-geostrophic barotropic models and using them in weather forecasting operations, Chinese scientists also studied quasi-geostrophic baroclinic models as well as primitive equation models. During his visit (1959)(1960) to the Institute of Applied Geophysics of the Academy of Sciences of the USSR, Zhu (1961Zhu ( , 1962) developed a three-layer nonlinear model in a spherical coordinate system at 1000-, 500-, and 300-hPa levels. He investigated the topographic effects on 24h forecasting. In 1961, five years before the first use of atmospheric primitive equations in USA, Zeng, for the first time, produced real-case weather forecasts using the primitive equations and implemented the relevant operation at the Moscow World Meteorological Center (Zeng, 1961(Zeng, , 2013. Zeng proposed a "semi-implicit (or semiexplicit)" difference scheme based on the theory of the adaptation process of atmospheric motions for solving the primitive equation models. This approach improved computational stability and significantly reduced computational cost by implicitly solving the fast waves (Zeng, 1961(Zeng, , 1963a. Several years later, Robert published his paper on the semi-implicit algorithm (Robert, 1969), and further proposed the semi-implicit semi-Lagrangian (SISL) scheme in 1985 (Robert et al., 1985). To this day, the SISL remains being extensively used in operational NWP models in the world. In view of the computational complexity introduced by topography, Zeng (1963b) proposed a hydrostatic subtraction (or standard stratification subtraction) method in 1963. This method greatly reduces computational errors of pressure gradient forces resulted from the subtraction of two major terms, and later played an important role in reducing computational errors due to steep topography in spectral models, as well as in semi-Lagrangian advection calculations (Simmons and Chen, 1991;Temperton et al., 2001).
In the 1950s and early 1960s, Chinese scientists also conducted extensive investigations in dynamic meteorology, atmospheric circulation, dynamic and thermodynamic effects of the Tibetan Plateau, adaptation theory of the atmospheric motions, and the atmospheric dispersion theory (Xu, 1959;Zhao, 1959;Ye, 1963). These achievements directly or indirectly promoted the development of NWP in China (Ji et al., 2005). Moreover, there were considerable excellent studies in cloud microphysics, cumulus dynamics, mesoscale dynamics, and boundary layer turbulence (Su, 1958(Su, , 1959Zhou, 1963;Chao and Zhou, 1964), providing comprehensive theoretical support to understanding the physical processes in NWP. During this period, it should be mentioned that a novel approach different from the initial value problem of NWP was proposed. This approach suggested that NWP be formulated as a problem by using the recent historical data rather than only using the initial value (Gu, 1958a, b). Later, Chou further developed this idea into the functional minimization problem using multi-time historical observations (Chou, 1974). The approach proposed by Gu and its further extension into functional minimization problem by Chou are the general framework of the present four-dimensional variational assimilation (4D-Var) in NWP.
The 1970s is an important period when NWP was officially put into operation internationally, and gradually became an indispensable scientific support for daily weather forecasting (Kalnay, 2002). The founding of the ECMWF in 1975 was a milestone during this period. Though lacking environment for systematic research at this time for Chinese scientists, fundamental studies on NWP had been conducted (Tao et al., 2003). The NPT at the Institute of Atmospheric Physics (IAP) of the CAS followed international progress in a timely fashion and introduced the then-advanced achievements in NWP made by other countries to China (Numerical Prediction Research Team, 1975). The Mathematical and Physical Basis of Numerical Weather Prediction (Volume 1) by Zeng in 1979 is a representative work of this period (Zeng, 1979). In this work, Zeng presented a sound analysis and systematic exposition on the mathematical and physical principles of NWP.
The Chinese economic reform of the late 1970s and 1980s presented new opportunities for research and operational applications of NWP in China. During that period, research conducted by Chinese scientists was concentrated primarily on computational scheme-related issues, including the theoretical problem of computational stability and the design of difference scheme (Zeng, 1978;Zeng and Ji, 1981;Ji and Zeng, 1982). Since the 1990s, Chinese scientists have also made internationally influential or leading advances in areas such as predictability, new algorithms for numerical atmospheric models, and new data assimilation methods. In addition, development of NWP models was also active following the Chinese economic reform, and played a crucial role in the establishment of modern NWP operations in China (Tao et al., 2003).
Designing a suitable computational scheme to allow discretized model equations to satisfy certain physical constraints is one of the key issues in developing numerical models. These physical constraints generally include conservation of total energy, conservation of total mass, and conservation of enstrophy. Such properties of the discretized equations, on the one hand, ensure the reasonable long-term integral performance of model; on the other hand, can also help overcome computational instability. In 1966, Arakawa (1966Arakawa ( , 1972 demonstrated that conservation of mean kinetic energy and enstrophy can be maintained by designing a finite difference Jacobian expression for the advection term, which properly represents the interaction between grid points. This is a classic work in the field of computational physics. Chinese scholars noted the problem in the Arakawa's instantaneous conservation scheme (Ji, 1981;Ji and Zeng, 1982). Ji and Zeng (1982) proposed an implicit difference scheme for advection term. This scheme is computationally stable under temporal discretization and any initial conditions, and it conserves the energy, the "generalized energy," and the "average scale" (Ji and Zeng, 1982). While the implicit total energy-conservative scheme is computationally stable, it is difficult to solve in practice due to excessive computing cost. On this basis, Wang and Ji proposed an explicit scheme with conservation of quadratic quantities of physical importance (Wang and Ji, 1990;Ji and Wang, 1991). Based on this scheme, Zhong (1992Zhong ( , 1993 further developed a third-order scheme. Some researchers also explored the possible application of fully implicit difference schemes (Chen et al., 2007). Investigations on or related to these computational schemes consistently demonstrate the importance attached by Chinese scientists to basic numerical computation methods. These achievements are believed to be the important accumulation for developing the general circulation models in IAP, CAS (Zeng et al., 1985;Wang et al., 2004).
Moreover, studies on predictability have also made advances with international influence. Mu et al. (2003) popularized the linear singular vector (LSV) to the nonlinear case, and proposed conditional nonlinear optimal perturbation approach (CNOP). LSV represents the type of initial perturbations with the highest growth rate in tangent linear models. CNOP represents the type of initial perturbations with the greatest nonlinear development at the time of prediction, which satisfy certain physical constraints. Mu et al. (2010) further extended CNOP to situations where initial and model parameter errors coexist, and referred to CNOP related only to initial perturbations as CNOP-I and to CNOP related to model parameter perturbations as CNOP-P. The CNOP has been applied to a broad range of applications in predictability research at various timescales, target observation, and ensemble forecasting.
Operational spectral and grid-point models at major NWP centers in the world all adopt spherical Gaussian or latitude-longitude grids and the SISL time integration scheme. With the development of multi-core high-performance computers, computing bottleneck becomes stand out in terms of scalability, which results from the discretization algorithms and the grid structures utilized in the current models. The scalability becomes an urgent issue in the future high-resolution conservative global model on the massively parallel computers. The scalable dynamic core with high accuracy and conservative property is one of the most important challenges for developing the next generation NWP model (Mengaldo et al., AUGUST 2020 2019). Shen and his research team developed a finitevolume discretization algorithm with multi-moment constraints for developing the next generation atmospheric model in the China Meteorological Administration (CMA). This method ensures strict numerical conservation. In addition, compared to the conventional finitevolume method, much more local freedom, flexibility to various grid structures, as well as locally intensive computing characteristics can be achieved Chen et al., 2014). This new method has received interests from international peer experts who are developing the same type of next-generation models (Smolarkiewicz et al., 2016). Moreover, toward the future NWP models at grey-zone scale, investigations on saleaware model physics have also made progress. Huang et al. (2018) developed a scale-aware cumulus convection schemes based on the concept proposed by Arakawa. Zhang et al. (2018) developed a three-dimensional (3D) scale-adaptive turbulence parameterization scheme. They all show reasonable simulations from hundred meterscale and kilometer-scale resolutions to tens of kilometer resolutions.
Advanced data assimilation techniques are widely recognized to be one of the main contributors to the significant improvement of the NWP forecast skill (Bannister, 2017). In China, systematic studies on data assimilation began relatively late but developed rapidly. In the 1980s and 1990s, optimal interpolation objective analysis and nonlinear normal mode initialization had been developed and implemented into operation (Xue et al., 1992;Zhu et al., 1992;Tu and Zhang, 1995). Since 2000, variational data assimilation technique (Var) started to be developed (Zhang et al., 2019); of which, the Var system of the Global/Regional Assimilation and PrEdiction System (GRAPES) is a representative achievement, and has been implemented into operation to serve the daily NWP for more than 10 years (Xue and Chen, 2008).
Chinese scientists also proposed new approaches to advancing or developing the 4D-Var and ensemble Var techniques. The 4D-Var algorithm based on a dimensionreduced projection (DRP4D-Var) proposed by  and the nonlinear least squares ensemble 4D-Var algorithm (NLS-En4D-Var) proposed by Tian et al. (2018) are representative works. The DRP4D-Var uses ensemble technique to replace the adjoint model in calculating the gradient of the cost function. The model variable space and the observational variable space are dimensionally reduced and projected to the sample space. The principal modes of the localized correlation function are used to expand the sample space. Finally, a 4D-Var analysis is carried out in the expanded sample space. The DRP4D-Var avoids tangent linear and adjoint models in the the minimization procedure, achieving a significant decrease in computing costs. In the NLS-En4D-Var algorithm, the advantages of 4D-Var and the ensemble Kalman filter are mutually complemented. The NLS-En4D-Var algorithm can be strictly derived to transform an ensemble 4D-Var problem into a nonlinear leastsquares optimization problem, which integrates the nonlinear least-squares optimization problem and the ensemble data assimilation algorithm together. These two new 4D-Var algorithms have demonstrated high potential for operational application.
The period from the late 1970s to 1990s is a small climax for developing the NWP models in China. In 1980, a three-layer primitive-equation adiabatic model covering the Asian and European regions was jointly developed and implemented into operation by the joint team of the CMA and IAP of the CAS. This model had been used for supporting the 48-h weather forecasting, being referred to as Model A. In 1982, the CMA operational models were upgraded to a hemispheric model and a regional five-layer grid-point model, known as model B and small model B, respectively (Tao et al., 2003). The model B series was the joint efforts by scientists from IAP, Peking University, and the National Meteorological Center (NMC) of CMA. In 1995, based on the the basic framework proposed by Yuling Zhang of Peking University, a new limited area forecast system (LAFS) was implemented in NMC, later developed into a high-resolution LAFS (HLAFS; Guo et al., 1995). These models supported daily weather forecasting as well as disaster prevention and reduction in China in the 1980s and early 1990s (Chen and Xue, 2004). The HLAFS had been in operation until 2006 when it was replaced by the GRAPES mesoscale NWP system (GRAPES_Meso).
In addition, models were also independently developed in Shanghai, Guangzhou, and Lanzhou for their respective local applications, e.g., the typhoon model of Shanghai (Zhu and Yin, 1987), the tropical weather forecast model of Guangzhou (Xue et al., 1988), and the model suitable for complex terrain of Lanzhou (Yan, 1987). In view of China's complex terrain and rainstorm forecasting, based on Mesinger's work (Mesinger, 1984), Yu developed a prediction model that accounts for steep topography by introducing η-coordinates in 1986, which was later referred to as the regional eta-coordinate model (REM; Yu et al., 1994). By updating and improving the physical processes over the years, the REM has further evolved into an advanced REM (AREM). The AREM was once used in scientific research and operations, such as in the Wuhan Institute of Heavy Rain (Yu et al., 2004). During this period, apart from the operation-oriented NWP models, there were several individually de-

678
Journal of Meteorological Research Volume 34 veloped research models in the meteorological research community. Zheng (1980Zheng ( , 1989 developed a seven-level hemispherical spectral model with a spherical function as the basis function. This model utilized the full-form vorticity and divergence equations, and successfully produced a batch of real-case forecasts. Zheng's work was approximately concurrent with international research on spectral models. In Zheng's model, the nonlinear terms were calculated by using the classical "interaction coefficient approach" rather than the currently popular "transformation method." Moreover, investigations on non-hydrostatic models were also conducted relatively early in China, such as the non-hydrostatic regional model by Wang and Zhou (1996) and the 3D non-hydrostatic model by Hu and Zou (1992). At that time, these non-hydrostatic models were limited to research only. It is worth noting that the present paper focuses on reviewing the efforts and achievements in NWP, not including long-term integrated atmospheric circulation and climate models. See elsewhere for the independent and innovative development of atmospheric circulation models [Zeng et al. (1985) for the IAP AGCM-I; Wang et al. (2004)

Development of the modern NWP operational system in China
As mentioned above, while Chinese scientists had carried out exploratory experiments on NWP operations in the 1950-1960s, it was difficult to build operational capability due to the then conditions. Operational capability for short-range forecast was established after 1980 with the applications of the models A and B. The development of global medium-range NWP operational system began in the mid-1980s, and was approved and supported by the Ministry of Science and Technology of China (MOST) in 1985 as a National Key Science and Technology Project during the "Seventh Five-Year Plan" period (Li, 2010). In view of the large gap between China and advanced countries in terms of research and development capacity and supporting conditions at that time, a technology roadmap involving the importation of full prediction models from advanced countries was approved. A cross-department joint task team consisting of nearly 300 scientists from 17 organizations/institutes was established and led by Zechun Li (Li, 2010 ). Based on the global spectral model introduced from the ECMWF, a T42L9 global prediction system was implemented into operation in 1991 on the Chinese-manufactured Galaxy supercomputer. This system had been upgraded several times in the 1990s and 2000s to the T63L16, T106L19, T213L31, and T L 639L60 successively (later referred to as the T series). It should be pointed out that the T series had played a seminal role in nationwide weather forecasting operations and services, and laid an important foundation for subsequent GRAPES development.
At the turn of the 21st century, the CMA made an important decision to shift the technology roadmap from primarily importing technologies to independent development, and thereby started developing the new-generation NWP system in China. In 2001, the project entitled "Innovative Research on Meteorological Numerical Prediction System in China" was approved by the MOST Key National Science and Technology Project of the "Tenth Five-Year Plan." This five-year project succeeded in developing a preliminary version of the unified GRAPES (Xue and Chen, 2008). Later, under the continuous support from the MOST and CMA, more than 10 years was taken to establish an integrated GRAPES operational system, including deterministic and ensemble NWP models on horizontal resolutions of 3-10 km for regional and 25-50 km for global forecasting. Especially, a specialist team has been trained and growing in various fields of NWP. During the development of GRAPES, significant progress has been made in several core technologies of NWP, such as the non-hydrostatic dynamic core, 4D-Var, satellite data assimilation, scalable numerical algorithms of high-accuracy for atmospheric models, and the cloud microphysics schemes. These prove to be a solid accumulation in line with scientific and technological advancement, facilitating future development of NWP in China.
Earlier research and development of NWP in China at various stages have been summarized in Chen and Xue (2004), Li (2010), and Shen et al. (2010). The present paper reviews and discusses in a comprehensive manner the major scientific and technological progress, operational applications, as well as the future roadmap of NWP in China, with an emphasis on the achievements in developing the GRAPES since 2010.

Research and operation of GRAPES
The GRAPES consists of the regional and global unified models and associated assimilation systems. At its development stage in early 2000, there was a trend internationally to develop the non-hydrostatic grid-point atmospheric models. It should be noted that the route to select the grid-point model for designing the GRAPES was scientific and reasonable at that time. In addition, the AUGUST 2020 strategy adopted to develop assimilation systems, which involved the use of 3D-Var first and a gradual transition to 4D-Var, also conformed to the pattern of scientific development.

Milestones of GRAPES development and operation
By the end of the first 5-yr project (2001)(2002)(2003)(2004)(2005), an SISL non-hydrostatic dynamic core, physical processes suitable for mesoscale NWP, regional 3D-Var at the isobaric surface, and parallel computational supporting software had been achieved. Based on these, a mesoscale NWP system, GRAPES_Meso 2.0, was established and put into operation in July 2006. The horizontal resolution was 30 km. This system replaced the HLAFS as the operational mesoscale prediction system at NMC of the CMA as mentioned above.
The GRAPES_Meso 2.0 was the first operational application of the GRAPES achievement. In 2007, this system was upgraded to the GRAPES_Meso 2.5 by increasing its horizontal resolution from 30 to 15 km. Systematic problems, such as scattered forecast rainbands and underestimation of precipitation intensity, were discovered during the operational application. In view of these problems, a high-accurate conservative scalar advection scheme based on piecewise rational method (PRM) and a monotonic horizontal diffusion scheme that accounts for influence of steep topography were developed. In addition, effective topography was taken into account according to the model effective resolution. These improvements significantly improved the simulation accuracy for water substances, and much better precipitation forecast was obtained (Shen et al., 2011;Shen and Zhou, 2013). The PRM advection scheme also significantly improved the typhoon intensity forecast in the GRAPES Regional Typhoon Prediction Model (GRAPES_TYM). These achievements resulted in the first important upgrade of the GRAPES_Meso from version 2.5 to 3.0 in February 2010.
The mesoscale data assimilation system is another key issue that affects the forecasting performance of the GRAPES_Meso. Previously, the geopotential height, stream function, velocity potential, and water vapor (specific or relative humidity) on the isobaric surface were used as the analysis variables in GRAPES_Meso 3D-Var. This was incompatible with the prediction model. Accordingly, an additional interpolating process was required, regardless of whether the prediction model provided a background field or the assimilation system provided analysis to the model as an initial field. This inevitably introduced additional errors (Xue and Liu, 2012). This problem also existed during the early stage of the GRAPES_GFS development, and will be mentioned later. To address the aforementioned problem, the 3D-Var system was rewritten and reconstructed at the model grid space for both global and regional uses (Xue and Liu, 2012). This was also an important step that led to the subsequent success of the GRAPES global 4D-Var. Based on this new 3D-Var system, Wang R. C. et al. (2018) further updated the zonal and meridional winds as control variables and achieved direct assimilation of temperature and surface pressure. Moreover, they rendered it possible to apply the 3D-Var to the kilometer-scale NWP by optimizing the background error covariance (BC), and introducing new balance constraints and horizontal correlation models as well as large-scale constraints. Meanwhile, a cloud analysis system was developed based on the scheme in the ARPS (Advanced Regional Prediction System), which is capable of using minute-level sounding data, ground observations, Feng-yun (FY) geostationary meteorological satellite cloud products, and 3D Doppler weather radar mosaic data (Zhu et al., 2017). These new developments were important scientific and technological advancements since the GRAPES_Meso 3.0.
The above work significantly facilitated continuous improvement of the GRAPES_Meso forecast skill and led to its diverse applications. A rapid assimilation and forecast system (GRAPES_RAFS) and a regional typhoon prediction system (GRAPES_TYM) were brought into operation in 2012. In 2014, the resolution of the GRAPES_Meso was upgraded from 15 to 10 km , and the vertical resolution was also increased from 30 to 50 levels (especially in the boundary layer). In 2016, a GRAPES_Meso quasi-operational system with 3-km horizontal resolution was implemented for better forecasting the severe weather in eastern China. This system was further expanded to cover the whole country in 2019.
The GRAPES_GFS was started to develop almost from scratch since July 2007. On the basis of the SISL non-hydrostatic dynamic core and the isobaric-surface regional 3D-Var, the GRAPES global model and global assimilation system were developed. The satellite data assimilation, physics optimization, as well as the computational stability of the system were the main challenges. In March 2009, a testing system in the NMC operational environment was built up, and soon became the quasi-operational one (Shen et al., 2009). This is the first version of GRAPES_GFS, and referred to as GRAPES_GFS 1.0, in which the model horizontal resolution is 0.5° × 0.5°a nd there are 36 levels in the vertical. Main achievement during this period is not only the knowledge accumula-tion of global NWP, but also the establishment of an NWP-oriented operational receiving and processing procedure for global multi-satellite platform data in CMA. This is an excellent example of how the independent and innovative development of NWP drove full-chain progress in meteorological operations.
It is noted that the classic two-time-level SISL algorithm was used in developing the GRAPES (Gospodinov et al., 2001;Temperton et al., 2001). For the purpose of maintaining the computational stability, a relatively large semi-implicit weight coefficient was adopted in the model, so the accuracy of time discretization was only at the first order. In addition, time extrapolation was required when calculating departure points and the nonlinear terms on the right-hand side of the equations. This led to computational noise, and further resulting in computational instability in some cases. Moreover, due to the use of non-conservative equation set and the semi-Lagrangian algorithm, the conservation of model global mass was a challenging problem, which was poorly accounted for in designing the algorithm. Data assimilation was another issue that significantly affected global forecast skill. As mentioned previously, the assimilation system of the GRAPES_GFS 1.0 is a 3D-Var on the standard isobaric surface, with geopotential height, stream function, velocity potential, and relative humidity as analysis variables on the Arakawa A-grid. The model initial value was derived through variable transformation, hydrostatic relationship, and interpolation. In particular, interpolation errors inevitably limited the accuracy of initial value. The above issues affected the GRAPES_GFS forecasting performance to a large extent since its quasioperation.
Three problem-tackling stages followed the quasi-operation of the GRAPES_GFS 1.0. The first stage was represented by the improvements to the physical processes, including the introduction of the Rapid Radiative Transfer Model (RRTM) for General Circulation Models (RRTMG), the Common Land Model (CoLM), and the subgrid orographic gravity wave drag parameterization, as well as the optimization of cumulus convection schemes. These improvements equipped the GRAPES global model with a set of physical processes suitable for global forecasting. During the second stage, the focus was placed on tackling the GRAPES 3D-Var for avoiding errors caused by interpolation due to the different grid spaces used in Var and the model (Xue and Liu, 2012). In addition, activities in better use of satellite data were also carried out so that much more satellite data can be assimilated. The third stage was represented by diagnosing and optimizing the cyclic system of assimilation and prediction. During this stage, several key issues were solved, which greatly promoted the operation of GRAPES_GFS. These key issues include problems existing in dynamic core and physical processes, as well as several detailed items in the assimilation system. These issues affected the stability of the cyclic system and analysis quality to a certain extent. More details can be found in the overview paper (Shen et al., 2017).
The above improvements promoted the upgrade of GRAPES_GFS and its official operation. A new version with a horizontal resolution of 0.25° × 0.25° and 60 layers in the vertical was developed and approved for operation at the end of 2015 (Shen and Wang, 2015), hereafter referred to as GRAPES_GFS 2.0. The operational performance of GRAPES_GFS 2.0 so far demonstrates that it has overtaken the imported T L 639 system in terms of the statistical indices and precipitation threat scores at various thresholds. However, it should be noted that the scores of geopotential height exhibit some gap compared with the advanced centers such as ECMWF and NCEP.
After the successful operation of GRAPES_GFS, a momentous update of this system is the implementation of the GRAPES 4D-Var. The details of the assimilation framework were further addressed, including the computational efficiency, time slot of observations, multi-outerloop configurations, wind-pressure balance relationship, BC, etc. In particular, linearized physics, including vertical diffusion, deep cumulus convection, subgrid orographic blocking drag, and large-scale condensation, were developed. For improving the satellite radiance assimilation, a dynamic bias correction technique was implemented. See elsewhere for the scientific and technical details of the GRAPES 4D-Var system (Zhang et al., 2019). In July 2018, the GRAPES_GFS 2.4 with the 4D-Var system was put into operation.
Ensemble forecast is an important part of the integrated operational NWP system. The ensemble prediction systems (EPSs) for the medium range and short range are quite important to address forecast uncertainties at synoptic and meso scales. In 2010, the Numerical Weather Prediction Center of the CMA (NWPC) began to develop the GRAPES mesoscale and global ensemble prediction systems, hereafter referred to as the GRAPES Regional Ensemble Prediction System (GRAPES_REPS) and GRAPES Global Ensemble Prediction System (GRAPES_GEPS). For the GRAPES_REPS, initial perturbations are generated by using an ensemble transform Kalman filter (ETKF). The ETKF provides a set of perturbations, which are added to the global analysis to provide the initial values for ensemble members. Model uncertainties are represented by two particular types of AUGUST 2020 stochastic physics schemes, including the Stochastically Perturbed Physics Tendency (SPPT) Scheme and Stochastic Kinetic Energy Backscatter (SKEB). The GRAPES_REPS on the 15-km horizontal resolution was brought into operation in 2015. In 2019, the GRAPES_ REPS was further upgraded to 10-km horizontal resolution and 50 layers in the vertical. For the GRAPES_ GEPS, the singular vector (SV)-based initial perturbation method was adopted based on the GRAPES global tangent and adjoint models, and model uncertainties were represented by using the SPPT. The GRAPES_GEPS with a resolution of 50 km and 60 vertical layers was put into operation in 2018.
On the basis of the above achievements, an integrated NWP system for deterministic and ensemble forecast with resolutions in the range 3-10 km for regional forecasting and 25-50 km for global forecasting was established by the end of 2018. In addition, a full-chain support system was also built up.

Main scientific and technical characteristics
The imported T series had been used in weather forecasting as well as disaster prevention and reduction in China for over two decades since 1991, and was deeply influential to weather forecasters nationwide. Therefore, it is not easy to establish a brand new operational system based on the GRAPES to replace all functions of the T series for forecasting operations. The description of the GRAPES milestones in the previous section shows that long-term team efforts and accumulation of knowledge and experience are very important for research transferring to operation and for continuously improving its operational capability. This is also the case at the advanced NWP centers (Bauer et al., 2015).
Over the 13 years of continuous improvement of the GRAPES, significant advances have been made in various aspects, such as the dynamic core, physical processes, assimilation, and applications of various observational data. This section devotes to summarizing these new progresses, while the earlier research and development achievements can be found in the book Scientific Design and Application of the GRAPES Numerical Prediction System by Xue and Chen (2008). Table 1 summarizes the progress in development of the GRAPES dynamic core and model physics. The comparison with the preliminary version is also given. An iterative predictor-corrector SISL algorithm, similar to that currently implemented in the operational NWP models of the UK Meteorological Office (UKMO) and the Meteorological Service of Canada (MSC), has been developed, together with the replacement of the isothermal reference by a 3D state. These revisions greatly improved the computational accuracy and stability of GRAPES, and largely mitigated the computational noise problem caused by classic two-time-level SISL algorithm. Moreover, marked progress has also been achieved in the research and development of the physical processes. Compared with GRAPES 1.0, the physics package was completely updated both in the completeness of physics and new development. Also, as mentioned in the previous section, a lot of optimization activities have been carried out. In particular, a double-moment microphysics and a prognostic cloud cover scheme were developed and implemented successfully (Ma et al., 2018). These efforts in physical processes have become one of the important factors in enhancing the regional and global forecast skills.  Table 2 is given to show the advances in data assimilation system. Compared with the isobaric-surface 3D-Var in the early-stage GRAPES Var system, the GRAPES global 4D-Var and a regional 3D-Var suitable for kilometer-scale NWP have been established. The significant progress on Var systems in the past 10 years is the same big contributor as the model improvement to the forecast skill. It should be pointed out that the operation of global GRAPES 4D-Var is a big step for further advancing the GRAPES NWP system.
Tables 3 and 4 summarize data application in the current GRAPES global and regional prediction systems, respectively. Careful use of observations in data assimilation system is quite critical in improving the NWP performance (Magnusson and Kallen, 2013). As shown in Tables 3, 4, the early GRAPES_Meso operational system used only conventional data. In comparison, the current GRAPES system is capable of assimilating more non-conventional data, including nationwide radar re-flectivity and radial wind data, ground-based Global Positioning System (GPS) precipitable water data, domestic and foreign geostationary and polar-orbiting satellite multi-platform multi-channel radiance data, and GPS occultation data. The increasing capability to process and assimilate various types of observations is one of the important achievements in the past 10 years. See Han et al. (2010), Li and Zou (2014), Liu and Xue (2014), Wang et al. (2016), and Wan et al. (2017) for details. 2.2.2 Overview and performance of the GRAPES operational system By the end of 2018, based on the GRAPES achievements, an integrated NWP operational system with appropriate configurations of high-and low-resolutions has been established, which includes the global and regional systems for both deterministic and ensemble forecasts. Compared with the T series based on imported technologies, it is revealed that development of a totally new operational system (GRAPES series) independently needed Table 2. Comparison of the current and early development-stage versions of the GRAPES in terms of the Var system GRAPES 1.0 regional isobaric-surface 3D-Var Operational system regional GRAPES 3D-Var  not only extensive operation-oriented transformation of research and development achievements both in models and assimilation, but also reconstruction and development of various sub-systems, from the data retrieval, preprocessing, and monitoring of observational data to the post-processing, verification, and visualization of model results. Apart from the marked progress in model, assimilation, and application of observations, substantial achievements have also been made in the development and application of support systems, such as database and management information system, testing platforms, and analytical and diagnostic tools. The operational GRAPES system consists of a global deterministic prediction system with a horizontal resolution of 0.25° × 0.25° and 60 vertical layers, a global ensemble prediction system of 0.5° × 0.5° resolution and 60 vertical layers, a mesoscale deterministic prediction system of 0.09° × 0.09° resolution and 68 vertical layers for the Asia-Pacific region, a high-resolution (0.03° × 0.03°) NWP system, and a mesoscale 0.1° × 0.1° ensemble prediction system with 50 vertical layers that covers the Chinese mainland. In addition, a GRAPES-driven global/ regional ocean wave prediction system, a sandstorm prediction system for East Asia, and a nuclear pollution emergency prediction system are also important components of the GRAPES system. Tables 5 and 6 summarize the configuration of the GRAPES deterministic and ensemble prediction systems, respectively. The details of model configurations, operational schedules, forecast length, resolutions, and ensemble approaches are presented in the two tables.
It is worth noting that there has been continual and rapid development of NWP in the world, and many countries have strengthened their investment in this technology in the past 20 years. While China has made great advances, there remains a certain gap between the GRAPES and advanced international operational NWP systems. This gap exists in several aspects of the NWP system, such as the much advanced data assimilation techniques including cloudy-and rainy-area satellite data assimilation, convection-permitting scale ensemble system, kilometer-scale data assimilation techniques, land surface data assimilation, etc. These important issues need to be addressed in the future. Figure 1 shows the year-on-year improvement of the precipitation forecast threat score (TS) by the GRAPES_Meso. Clearly, capability of precipitation forecast exhibits gradually better skills. This progress is considered mainly due to the resolution increasing, dynamic core and physics improvement, as well as effective application of radar reflectivity data through cloud analysis as mentioned above. The GRAPES_Meso is becoming an important, even indispensable, support to daily weather forecasting at NMC and other provincial as well as local weather offices. In June 2019, the GRAPES_Meso with 3-km horizontal resolution covering the Chinese mainland was put into operation. This system exhibits much better forecast skill for heavy rainfall events. Here, the average equitable TS (ETS) for 3-h rainstorm forecasts produced by the Shanghai Meteorological Service (SMS-WARMS based on WRF) and the ECMWF are used for comparison, as plotted in Fig. 2. The reason for such a comparison is that these two forecast products are the most popular products among forecasters. As demonstrated in Fig. 2, the GRAPES_Meso-3 km was significantly advantageous in the 36-h rainstorm forecasts, especially in the first 20-h forecast range. Figure 3 shows temporal evolution of the 5-day hindcast errors in typhoon track and intensity by GRAPES_ TYM. Figures 3a, b, and c are for the track, minimum central pressure, and maximum wind speed of the typhoons over the western Pacific during 2016-2018, respectively. For comparison, forecast errors during the same period given by the ECMWF IFS (Integrated Forecasting System) and NCEP GFS (Global Forecast System) are also included. These two global models are often used by forecasters as reference. In terms of the typhoon track, the ECMWF forecasts still prevail, while the forecasts by the GRAPES_TYM and NCEP are comparable. In terms of the typhoon intensity, the GRAPES_ TYM exhibits relatively obvious advantages in average errors in both minimum central pressure and maximum wind speed forecasts.
It is known that the global forecast system is the core  Shen, X. S., J. J. Wang, Z. C. Li, et al. of any integrated NWP system. The time series of the anomaly correlation coefficient (ACC) of 500-hPa geopotential height from the GRAPES_GFS forecast on day 5 over the Northern Hemisphere is plotted in Fig. 4 to show the performance of the GRAPES_GFS in the past 10 years. Here, the ACC from the ECMWF, NCEP, and the T series are also included for comparison. The gradually better forecast skill of GRAPES_GFS can be obviously observed, demonstrating the effects of the improvement and optimization of the system as reviewed in the previous subsection, although there remain some gaps when comparing with the ECMWF and NCEP forecasts. This will be discussed further in Section 4. Although there exist some gaps in terms of ACC scores, a number of cases show good performance of GRAPES_GFS in rainfall forecast. In particular, during the 2019 flood season in China, the GRAPES_GFS heavy precipitation forecast has surpassed that of the ECMWF in several typical cases. The GRAPES_GFS notably outperformed the ECMWF in forecasting heavy precipitation induced by the large-scale warm low-level wind shear, which occurred over the middle and lower reaches of the Yangtze River. For example, the 36-h forecast of the GRAPES_GFS for the case that occurred on 26 May 2019 produced a TS score of 0.29 for the rainfall amount larger than 50 mm, while the TS of the ECMWF forecast was 0.16. Some other heavy rainfall cases induced by the cold low-level wind shear also exhibited better skill relative to the ECMWF, such as the cases that occurred on 6, 7, and 12 June 2019 over the middle and lower reaches of the Yangtze River as well as over South China. Figure 5 shows scatter plots of the rainfall larger than 25 mm for the 36-h forecasts and rain-gauge observations over the Chinese mainland, in which the GRAPES_GFS forecast was initiated at 2000 Beijing Time (BT) 24 May 2019, and verified at 0800 BT 26 May 2019. Clearly, the GRAPES_GFS outperformed the ECMWF in forecasting this heavy precipitation event.
The operational application of the GRAPES global and regional ensemble prediction systems has added a wealth of probabilistic forecasting products to the GRAPES "big family." These systems play a vital role in providing the probability of large-scale heavy precipitation a week ahead and in producing quantitative and fixed-point probability forecasts for strong weather conditions over local areas.

Redesign of the GRAPES model
During the initial design and development of the GRAPES dynamic core, the classic two-time-level SISL algorithm, which was the mainstream algorithm used in the then international operational NWP models, was adopted to address the stability issue in advection and fastwave calculations. While the two-time-level SISL algorithm is notably advantageous in saving computational memory and improving computational efficiency, it requires that several issues be paid special attention; namely, use of an isothermal reference atmosphere during linearization, selection of a weight coefficient in semi-implicit formulation, as well as temporal extrapolation of the nonlinear terms and the wind speed at midpoint of trajectory when calculating departure point. These issues are closely related to the computational accuracy and stability of the two-time-level SISL algorithm (Ritchie and Tanguay, 1996;Simmons and Temperton, 1997;Temperton et al., 2001), and were not carefully studied during the initial development stage of GRAPES. Recently, a predictor-corrector SISL algorithm has been developed for redesigning the GRAPES model, and a 3D reference atmosphere has been introduced into the linearization process (Côté et al., 1998). The use of the 3D reference reduces the order of magnitude of the nonlinear terms, compared with that using the isothermal state; and the iterative algorithm avoids the time extrapolation. These make the semi-implicit weight close to 0.5, and the time integration reaches nearly two-order accuracy. This algorithm recon-

Development and operational application of the GRAPES global 4D-Var
Since the ECMWF successfully realized operational application of 4D-Var in 1997, 4D-Var has been re-

Journal of Meteorological Research
Volume 34 garded as the advanced data assimilation technique in improving NWP skills (Rabier et al., 2000;Bannister, 2017;Bonavita et al., 2017;Kwon et al., 2018). Due to the enormous workload and difficulty involved in developing tangent linear and adjoint models required in 4D-Var, developing 4D-Var is a challenging task. The development of the GRAPES global 4D-Var system began in 2008, and its operational application was achieved 10 years later, in July 2018. Based on the global non-hydrostatic GRAPES model, the non-hydrostatic tangent linear and adjoint models were developed. To maintain a relatively high similarity of evolution trajectory between the tangent linear and nonlinear models within the assimilation time window, tangent linear physical processes were also developed, including cumulus convection, cloud microphysics, vertical diffusion, and subgrid orographic blocking drag. Such development work forms an important basis for the successful implementation of the GRAPES global 4D-Var. As mentioned previously, the GRAPES global 4D-Var uses an incremental analysis scheme, in which digital filtering is introduced into the cost function as a weak constraint to suppress high-frequency gravity waves. Multiple outer loops have been designed in the system (a single outer loop is currently used in operations), and the Lanczos-conjugate gradient (CG) algorithm is used as the minimization algorithm . The 4D-Var framework adopts a horizontally and vertically inseparable BC model. A secondorder auto-regressive approach is used as a horizontal correlation model, and its correlation scale varies with height. A vertical correlation model is directly obtained statistically by the ensemble samples. When it comes to equilibrium constraints, the equilibrium between rotational and divergent winds, between rotational winds and mass fields, and between non-equilibrium divergent winds and mass fields, is taken into account and achieved by a dynamic and statistical-combined scheme. The GRAPES 4D-Var system runs 4 times per day, and its initial time of the time window is 0300, 0900, 1500, and 2100 UTC, respectively. The 6-h forecast at each initial time is used as the background field of next 4D-Var cycle, whereas the 3-h forecast within the cycle is the initial condition of 10-day forecast at 0600, 1200, 1800, and 0000 UTC, respectively. Multiple types of observations can be assimilated in the GRAPES global 4D-Var system, including those from surface stations, aircraft, radiosondes, ships, retrieved products from satellites, as well as radiance data. For details, please refer to Table 3. The quantity of observational data used in the GRAPES global 4D-Var system is approximately 50% more than that used in the 3D-Var system.
If we take the analysis of ECMWF as a standard to see the quality of GRAPES assimilation, the GRAPES global 4D-Var system comprehensively surpasses the 3D-Var system, particularly in mass and wind fields. In terms of mean-square error, the background field (6-h forecast field) of the GRAPES global 4D-Var system is even better than the analysis of the global 3D-Var. There is a more pronounced improvement in the geopotential height in the Southern Hemisphere and high-level wind field in the tropics. Moreover, the GRAPES global 4D-Var system outperforms the 3D-Var system in terms of various forecast validity periods and variables. As demonstrated in Fig. 6, from the perspective of the ACC for 500-hPa geopotential height, the skillful forecast length with ACC larger than 0.6 has increased by 5 h in the Northern Hemisphere and by 7-8 h in the Southern Hemisphere on average. For precipitation forecast, the GRAPES global 4D-Var system has a significantly greater ability to forecast moderate and above-moderate levels of precipitation. In addition, both typhoon track and intensity forecasts are evidently improved, and the track errors can be reduced by about 15% for up to 5-day forecasts.

Development of satellite data assimilation technique and its improved application
Since the NCEP and the ECMWF successfully assimilated satellite radiance data of TOVS (Television and Infrared Observation Satellite Operational Vertical Sounder) in their operational 3D-Var in 1995 and 1996, respectively, satellite radiances have received special attention for being an important data source for global NWP, as well as for their significant role in improving NWP skills (Eyre, 2007). The ECMWF is thought to lead in this field, and its proportion of assimilated satellite data is more than 90% among the total assimilated observations. The investigations on satellite data assimilation techniques in GRAPES began in 2004. Direct assimilation of NOAA-16 AMSU (Advanced Microwave Sounding Unit)-A radiance was investigated in the early GRAPES isobaric-surface 3D-Var system (Zhang et al., 2004), but the operational capacity was not built up. Much of the capability for assimilating satellite radiances started from the GRAPES_GFS development (Shen et al., 2009). The NOAA series (15, 16, 17, and 18) advanced TOVS data, GPS occultation data, and geostationary and polar-orbiting satellite AMVs (atmospheric motion vectors) data were implemented into the quasioperational version of GRAPES_GFS in 2009. Since that time, data from more advanced instruments, including the AIRS (Atmospheric Infrared Sounder), IASI (In-AUGUST 2020 frared Atmospheric Sounding Interferometer), GOES sounder, AMSU-A on AQUA and AMSU-A, HIRS (High Resolution Infrared Radiation Sounder), and MHS (Microwave Humidity Sounder) on Metop satellites, as well as the Chinese FY series satellites, have been implemented into the GRAPES_GFS. The successful assimilation of multi-platform satellite data is an important factor that leads to continuous improvements in the GRAPES_ GFS forecast skill (Shen and Wang, 2015). Figure 7 shows the year-on-year changes in the types as well as the platform numbers of satellite data assimilated in the GRAPES_GFS. As demonstrated in Fig. 7, rapid progress has been made in assimilating the satellite data since 2010, particularly since the full operation of version 2.0 of GRAPES_GFS in 2016 (Shen et al., 2017). Currently, satellite data account for approximately 70% of all the assimilated observations in the GRAPES_GFS (Table 3).
The above-mentioned progress in satellite data assimilation relies on the details of the assimilation system. These details not only include the variational assimilation techniques and background error covariance, but also each part in the chain of assimilation, such as data quality control, bias correction, cloud detection, channel selection, and so on. It is worth pointing out that a constrained bias correction (CBC) scheme for satellite radiances has been developed, which is used to cope with the systematic errors in the background. As is known, the  bias correction is necessary when there are systematic differences between the simulations and observations, which cannot be explained by systematic errors in the background. Several practical issues may exist in operational application. For example, when there is a significant bias in the background, the observational information will be erroneously corrected; imperfect quality control will result in unreasonable observation bias estimation and correction. These two issues are particularly prominent in developing the GRAPES_GFS, and are primarily reflected by the interaction of forecast bias and observation bias correction during assimilation-forecast cycle, which may result in a drift of the observation bias correction toward the model bias. Taking the idea of the "minimum module solution" method in the regularization of inverse problems in mathematical physics, a CBC technique for satellite data was developed by introducing uncertainties existing in satellite radiance calibration and in the radiative transfer model into the cost function in the form of an inequality constraints via regularization parameters (Han, 2014;Han and Bormann, 2016). This technique reduces the impact of the background errors on radiance bias correction, removes the inherent systematic bias in data itself, and allows better use of observational information.
The CBC is an original technique proposed during the development of the GRAPES. In 2016, Han and Bormann (2016) further developed it into the constrained variational bias correction (CVarBC). The CVarBC effectively improved the quality of analysis in the stratosphere by mitigating the drift of the satellite bias correction into the model bias in the stratosphere. The operational application of the CVarBC began in the 2018 ECMWF operational version (Integrated Forecasting System CY46R1), and was planned to be used in ECMWF's reanalysis of the 6th generation.
With continuous advances in China's FY meteorological satellites, the GRAPES operational system has realized the full use of FY products, including those from both polar orbit and geostationary FY series. Particularly, through the quantitative assessment as well as assimilation in GRAPES, a good collaboration and feedback mechanism were built up between the FY satellite group and the GRAPES team, greatly promoting the concomitant progress. The FY-2 AMV products have been continually improving its accuracy, and have been playing a significant role in GRAPES global/regional forecasting system (Han et al., 2006;Wan et al., 2017Wan et al., , 2018. The Visible and Infrared Spin Scan Radiometer (VISSR) is a payload on board of FY-2 satellites. VISSR clear-sky water-vapor-channel radiance data have been assimilated in the GRAPES global 4D-Var system (Wang, 2017), taking fully advantage of its high temporal resolution observations.
To comprehensively achieve the full potential of the multiple loads on the same platform of the FY-3 satellites, an integrated cloud and precipitation detection algorithm was developed. The MWTS, MWHS, MWRI, and GNOS data have been effectively assimilated in both GRAPES_GFS and GRAPES_Meso (Li and Liu, 2016).
The FY-4A satellite, a new-generation Chinese geostationary meteorological satellite, was successfully launched in 2016 and equipped with a geosynchronous interferometric infrared sounder (GIIRS). This satellite, for the first time in the world, achieved geostationary-orbit hyperspectral sounding (Yang et al., 2017). Several core techniques, including a GIIRS observation operator (Di et al., 2018), a channel selection technique based on the BC over the observation region (Yin et al., 2019), and an online bias estimation and correction algorithm suitable for large array detectors (Yin et al., 2020), were developed in the GRAPES. In December 2018, GIIRS radiance data were assimilated for the first time in the world into the operational GRAPES global 4D-Var. In addition, by taking full advantage of the flexible observation capability of the detection instruments on board FY-4A as well as the cloud forecasting and the sensitive area detection techniques of the NWP system, an "intelligent rapid clear-sky observation mode" and a "high-impact weather target observation mode" were realized. In 2018, in realtime operational settings, parameters such as satellite observation region and observation frequency were optimally determined based on the sensitive regions identified by using the GRAPES SV-based initial perturbation approach for Typhoons Maria, Ampil, and Mangkhut. Through collaborative interaction of observations and forecasts, a forecast-object-oriented geostationary satellite detection mode was realized, and these intensified observations were assimilated into the GRAPES global 4D-Var immediately. These object-observation experiments obviously improved the typhoon track forecasts compared with the operational GRAPES_GFS for the three typhoon cases.

Development of 3D-Var for kilometer-scale NWP
Developing kilometer-scale data assimilation technique is of great importance for improving the forecasts of extreme weather events and near-surface meteorological parameters (Gustafsson et al., 2018). A GRAPES kilometer-scale data assimilation system was developed based on the GRAPES global and regional unified variational assimilation framework (Xue and Chen, 2008).

AUGUST 2020
Several modifications have been made to reconstruct the original regional 3D-Var system for developing the kilometer-scale NWP . First, the control variables in the minimum cost function were rewritten from the stream function ψ, the potential velocity χ, and π, to the zonal wind u, meridional wind v, temperature T, and surface pressure p s . In addition, the geostrophic equilibrium between the control variables was no longer taken into account. The scale of spatial correlations in the background errors was much small, corresponding to the new control variables; and these variables were direct measurements of the observation system. Figure 8 is plotted to show the analysis increment in the operational 10-km resolution 3D-Var and the kilometerscale 3D-Var with 3-km resolution. Here, only a singlepoint temperature observation was assimilated. Obviously, the analysis increment of the kilometer-scale 3D-Var system exhibits more local distribution. Second, a multi-scale analysis scheme was implemented for keeping the analysis quality for both the synoptic-scale and meso-to-small-scale weathers. The multi-scale analysis was realized by using a multi-scale horizontal correlation model (Wu et al., 2018), and by merging the synoptic-scale and meso-to-small-scale components of analysis increments from the global analysis and the kilometer-scale analysis, respectively. Also, a weak constraint was introduced into the cost function to maintain the synoptic-scale information (Yang et al., 2019). In addition, to capture the short life time and quick evolution of mesoscale systems, the kilometer-scale 3D-Var was run by a rapid assimilation and forecast cycle of 1-3 h with digital filter to suppress the noise in the rapid refresh cycle.
Assimilation of observations with high spatiotemporal density is critical for the kilometer-scale NWP. In the kilometer-scale 3D-Var system, the radar radial wind and wind profiler data are directly assimilated. The radar reflectivity is applied through a cloud analysis scheme, which is used to estimate the hydrometeors and the latent heating rate. Moreover, the assimilation technique for the FY-4A AGRI (Advanced Geosynchronous Radiation Imager) and GIIRS data is also developed. Furthermore, as significant enhancements to the kilometer-scale 3D-Var system, assimilation of the automatic weather station (AWS) observations, especially those over the complex terrain, is under development..
The above specific developments for the kilometerscale 3D-Var have achieved good results. Figure 9 shows its effect in precipitation forecast. Here, the ETS scores are plotted, which are the averages over East China from 1-month experiments in July 2018. In comparison, the experiments by using the downscaling initial values from the GRAPES_GFS analysis are also given. It should be pointed out that both of the two experiments utilized the cloud analysis. Obviously, the kilometer-scale 3D-Var can effectively improve the short-range precipitation forecast for 0-12 h, indicating the great potential in the kilometer-scale rapid refresh applications.

Development of ensemble techniques and their operational applications
The GRAPES global ensemble prediction system (GRAPES_GEPS) adopted the SV-based initial perturbation approach. The SVs are a set of orthogonal perturba-  tions that grow the fastest within a given time interval, which are obtained through minimizing the cost function defined by a weighted norm based on the GRAPES tangent linear and adjoint models. Here, the weighted norm is defined as the total energy norm of the horizontal wind component , the perturbation potential temperature , and the dimensionless perturbation pressure π'. Evaluation of the GRAPES SV total energy and its kinetic as well as potential energy components illustrates its rationality in representing the baroclinic characters of perturbation growth in the troposphere of mid-high latitudes (Liu et al., 2013). Based on the obtained SVs, an initial perturbation field has been constructed by the Gaussian sampling approach . The GRAPES_GEPS was put into operation in 2018 . Table 7 summarizes the performance of the GRAPES_GEPS for the period January-May 2019. As demonstrated in Table 7, the GRAPES_GEPS forecast lead time at which the ACC of the 500-hPa geopotential height reaches 60% for the extratropical Northern Hemisphere is about 8.6 days, which is 1 day longer than the GRAPES_GFS. In the Southern Hemisphere, half day longer forecast lead time can be obtained in GRAPES_GEPS.
The ETKF-based initial perturbations are applied in the GRAPES regional ensemble system (GRAPES_REPS; Zhang et al., 2014). Based on this, a multi-scale blending techniques are further developed, which merges the small-to-meso-scale component of ETKF with the largescale component from the GRAPES_GEPS (Zhang et al., 2015). A three-dimensional bias correction scheme has been proposed for post-processing the GRAPES_REPS products (Wang J. Z. et al., 2018). In 2019, a GRAPES_ GEPS-driven 10-km regional ensemble prediction sys-tem was implemented into operation, showing good capability in providing the useful probabilistic information on extreme events.
Several stochastic approaches have been developed in formulating the model uncertainties for GRAPES_GEPS and GRAPES_REPS, such as SPPT, SKEB, and SPP. In applying these stochastic schemes, the key issue is how to construct the stochastic field. A first-order Markov process with time correlation and Gaussian distribution is adopted in both GRAPES_GEPS and GRAPES_REPS (Yuan et al., 2016). Application of these stochastic approaches improves the spread/standard deviation relationship in the two systems, especially for the probability forecast skill of heavy rainfall in the GRAPES_REPS, and the skill over the tropics in the GRAPES_GEPS (Peng et al., 2019;Xu et al., 2019). Ma et al. (2018) implemented a prognostic cloud fraction scheme in the GRAPES_GFS based on the Tiedtke scheme (Tiedtke, 1993). This scheme takes into account the contributions from boundary layer, cumulus convection, and cloud microphysics. This prognostic cloud fraction scheme exhibits better simulation of different cloud formation processes in various regions across the globe, compared with the diagnostic scheme previously used in GRAPES_GFS (Xu and Randall, 1996). Moreover, diurnal variation of oceanic stratocumulus is well simulated. Implementation of the new cloud scheme significantly reduces the biases of the underestimated low and high clouds. Moreover, the improvement in clouds also results in the improvement in radiation, obviously reducing the positive bias of outgoing longwave radiation caused by less cloud in the original diagnostic scheme.

Outlook
An integrated NWP operational system in China has been successfully built up, based on the GRAPES achievements. However, several important issues should be paid specific attention in the future development, including the continuous improvement of model accuracy, advanced satellite data assimilation technique, etc. As a state-of-the-art NWP system, it is necessary to implement the assimilation system of land surface parameters as well as sea temperature and sea ice into the GRAPES system. Meanwhile, investigations on the future NWP system based on the earth system approach should be enhanced.
As a short-term strategy, the top priority over the next two to three years is to continually improve the GRAPES system to meet the increasing requirements from the forecasters. The following five items are the key steps.
(1) The ensemble information will be taken into account into the GRAPES global 4D-Var. This is expected to further improve the global analysis quality.
(2) The advanced satellite data assimilation techniques will be developed to further enhance the better use of satellite data. The focus should be placed on developing variational quality control and variational bias correction techniques for satellite data, as well as the assimilation techniques for cloud/rain-affected radiances.
(3) The land surface data assimilation and snow analysis system will be developed. The improved forecast skill of near surface weather parameters can be expected.
(4) It is necessary to continuously improve the GRAPES dynamic core in terms of computational accuracy and efficiency as well as scalability. In addition, it is also necessary to continually optimize the physical processes of the GRAPES model. The much higher resolution (about 10-15 km in horizontal) of GRAPES_GFS and the kilometer-scale rapid refresh GRAPES_Meso are planned in the coming three years.
(5) A prototype of next generation NWP model is being developed. This next generation model is designed to have global and local mass conservations with high-order accuracy and high-scalability. These new developments will be strengthened further through more extensive collaborations among the NWP teams in China.