
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 nonhydrostatic gridpoint atmospheric models. It should be noted that the route to select the gridpoint model for designing the GRAPES was scientific and reasonable at that time. In addition, the strategy adopted to develop assimilation systems, which involved the use of 3DVar first and a gradual transition to 4DVar, also conformed to the pattern of scientific development.

By the end of the first 5yr project (2001–2005), an SISL nonhydrostatic dynamic core, physical processes suitable for mesoscale NWP, regional 3DVar 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 highaccurate 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 3DVar. 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 3DVar 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 4DVar. Based on this new 3DVar 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 3DVar to the kilometerscale NWP by optimizing the background error covariance (BC), and introducing new balance constraints and horizontal correlation models as well as largescale constraints. Meanwhile, a cloud analysis system was developed based on the scheme in the ARPS (Advanced Regional Prediction System), which is capable of using minutelevel sounding data, ground observations, Fengyun (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 (Huang et al., 2017), and the vertical resolution was also increased from 30 to 50 levels (especially in the boundary layer). In 2016, a GRAPES_Meso quasioperational system with 3km 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 nonhydrostatic dynamic core and the isobaricsurface regional 3DVar, 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 quasioperational 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° and there are 36 levels in the vertical. Main achievement during this period is not only the knowledge accumulation of global NWP, but also the establishment of an NWPoriented operational receiving and processing procedure for global multisatellite platform data in CMA. This is an excellent example of how the independent and innovative development of NWP drove fullchain progress in meteorological operations.
It is noted that the classic twotimelevel 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 semiimplicit 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 righthand side of the equations. This led to computational noise, and further resulting in computational instability in some cases. Moreover, due to the use of nonconservative equation set and the semiLagrangian 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 3DVar on the standard isobaric surface, with geopotential height, stream function, velocity potential, and relative humidity as analysis variables on the Arakawa Agrid. 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 problemtackling stages followed the quasioperation 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 3DVar 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 4DVar. The details of the assimilation framework were further addressed, including the computational efficiency, time slot of observations, multiouterloop configurations, wind–pressure balance relationship, BC, etc. In particular, linearized physics, including vertical diffusion, deep cumulus convection, subgrid orographic blocking drag, and largescale 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 4DVar system (Zhang et al., 2019). In July 2018, the GRAPES_GFS 2.4 with the 4DVar 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 stochastic physics schemes, including the Stochastically Perturbed Physics Tendency (SPPT) Scheme and Stochastic Kinetic Energy Backscatter (SKEB). The GRAPES_REPS on the 15km horizontal resolution was brought into operation in 2015. In 2019, the GRAPES_REPS was further upgraded to 10km 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 fullchain support system was also built up.

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 longterm 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 twotimelevel 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 doublemoment 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.
Dynamic core Physical process GRAPES 1.0 GRAPES operational system GRAPES 1.0 GRAPES operational system 1. Nonhydrostatic, fully compressible equations (shallow atmospheric approximation)
2. Isothermal reference atmosphere
3. Conventional twotimelevel SISL algorithm
4. 3D vector discretization of the momentum equations
5. Solving the potential temperature prediction equation using the conventional SISL algorithm
6. Quasimonotone positivedefinite scalar advection [quasimonotone semiLagrangian (QMSL)]1. Nonhydrostatic, fully compressible equations (shallow atmospheric approximation)
2. 3D reference atmosphere
3. Predictor–corrector SISL algorithm
4. 3D vector discretization of the momentum equations + W damping
5. Solving the potential temperature prediction equation using the vertical noninterpolation SISL algorithm
6. PRMbased highprecision conservative scalar advection1. NCEP microphysical scheme with 3–5 water substances
2. Betts–Miller cumulus parameterization scheme
3. Dudhia shortwave and RRTM longwave radiation schemes
4. Mediumrange forecast (MRF) PBL scheme
5. SLAB landsurface processes
6. Dynamics–physics feedback: linear interpolation1. Doublemoment microphysics or WRF single moment6 (WSM6) (for mesoscale model)
2. Cloud macrophysics and prognostic cloud scheme (for global model)
3. New simplified Arakawa–Schubert (NSAS) cumulus parameterization scheme or mesoSAS (for mesoscale model)
4. RRTMG longwave and shortwave radiation scheme
5. CoLM (for global model) or Noah (for mesoscale model) landsurface process scheme
6. New MRF (NMRF) or Yonsei University (YSU) (for mesoscale model) PBL scheme
7. Subgrid orographic gravity wave drag
8. Smallscale orographic turbulence form drag
9. Dynamics–physics feedback: cubic spline and quasimonotonic cubic spline interpolationTable 1. Comparison of the current operational GRAPES system and its early developmentstage version
Table 2 is given to show the advances in data assimilation system. Compared with the isobaricsurface 3DVar in the earlystage GRAPES Var system, the GRAPES global 4DVar and a regional 3DVar suitable for kilometerscale 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 4DVar is a big step for further advancing the GRAPES NWP system.
GRAPES 1.0 regional isobaricsurface 3DVar Operational system regional GRAPES 3DVar Operational system global
GRAPES 4DVarFormulation Incremental analysis Incremental analysis Incremental analysis Grid A grid in the horizontal direction and standard isobaric surfaces in the vertical direction C grid in the horizontal and terrainfollowing height coordinate in the vertical direction C grid in the horizontal and terrainfollowing height coordinate in the vertical direction Analysis variable Φ, ψ, χ, and RH/q T, p_{s}, u, v, and q/RH/RH^{*} π, ψ, χ, and q/RH/RH^{*} BEC NMC method NMC method NMC method Equilibrium relationship Linear balance equation Linear balance equation + statistics Linear balance equation + statistics Optimization algorithm Limitedmemory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm LBFGS algorithm Lanczosconjugate gradient (CG) algorithm Initialization Digital filter Incremental digital filter Incremental digital filter with weak constraints 1. Nonhydrostatic tangent and adjoint models
2. Linear physics (vertical diffusion, subgrid orographic blocking drag, deep cumulus convection, and cloud physics)Table 2. Comparison of the current and early developmentstage versions of the GRAPES in terms of the Var 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 nonconventional data, including nationwide radar reflectivity and radial wind data, groundbased Global Positioning System (GPS) precipitable water data, domestic and foreign geostationary and polarorbiting satellite multiplatform multichannel 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.
Observation type Instrument Platform Variable Proportion (%) Total proportion (%) Conventional data TEMP u, v, p, q 3.4 30.1 SYNOP p 3.1 SHIP p 0.2 BUOY u, v 0.1 AIREP u, v, T 23.3 Satellite data AMSUA NOAA15, 18, 19; MetOpA, B Radiance 9.6 69.9 MWHS2 FY3C, 3D Radiance 0.2 GIIRS FY4A Radiance 0.7 ATMS NPP Radiance 4.1 IASI MetOpA, B Radiance 39.8 AIRS EOS2 Radiance 2.3 GNSSRO COSMIC1; MetOpA, B; FY3C, 3D Refractivity 2.1 AMVsIR/
VIS/WVsFY2E, GOES13, GOES15, MTSAT2,
METEOSAT_10u, v 10.4 ASCAT MetOpA, B u, v 0.7 Table 3. Statistics of observational data used in the GRAPES_GFS during 0300 UTC 18–0900 UTC 27 July 2019
Observation type Instrument Platform Variable Proportion (%) Total proportion (%) Conventional data TEMP u, v, T, RH 4.22 6.24 SYNOP p, RH 0.74 SHIP p, RH 0.04 AIREP u, v, T 1.24 Unconventional data Doppler radar SA/SB/CB/SC/CD/CC VAD wind, 0.09 93.76 radial wind, 53.35 refractivity 20.85 WPR Wind 0.57 GPSPW Precipitable
water0.14 GNSSRO COSMIC1; MetOpA, B; FY3C, 3D Refractivity 0.11 AMVsIR/
VIS/WVsFY2E, GOES13, GOES15, MTSAT2, METEOSAT_10 u, v 0.36 Rain Rain (nudging) 18.29 Table 4. Statistics of observational data used in the GRAPES_Meso during 0600 UTC 1–0000 UTC 7 July 2018

By the end of 2018, based on the GRAPES achievements, an integrated NWP operational system with appropriate configurations of high and lowresolutions 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 not only extensive operationoriented transformation of research and development achievements both in models and assimilation, but also reconstruction and development of various subsystems, from the data retrieval, preprocessing, and monitoring of observational data to the postprocessing, 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 highresolution (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 GRAPESdriven 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.
GRAPES_GFS Mesoscale prediction system for the Asia–Pacific region (GRAPES_TYM) The kmscale highresolution prediction system (GRAPES_Meso3 km) Resolution 0.25° × 0.25°, L60 0.09° × 0.09°, L68 0.03° × 0.03°, L50 Model area Global; model top at 3 hPa 15°S−60°N, 40°E−180°; 10 hPa 10°−65°N, 70°−145°E; 10 hPa Forecast time 0000, 0600, 1200, 1800 UTC 0000, 0600, 1200, 1800 UTC 0000, 0600, 1200, 1800 UTC Forecast length 240 h 120 h 36 h Assimilation system 4DVar Downscaling of global analysis + vortex initialization Downscaling of global analysis + cloud analysis Physical process 1. Doublemoment cloud physics
scheme
2. Cloud macrophysics and prognostic
cloud scheme
3. NSAS cumulus parameterization
scheme
4. RRTMG longwave and shortwave
radiation scheme
5. CoLM landsurface processes scheme
6. NMRF PBL scheme
7. Subgrid orographic gravity wave
drag scheme
8. Smallscale orographic turbulent
drag scheme1. MesoSAS cumulus parameterization
scheme
2. WSM6 cloud physics scheme
3. RRTM longwave/Goddard shortwave
radiation schemes
4. Noah landsurface processes scheme
5. YSU PBL scheme1. WSM6 cloud physics scheme
2. RRTMG longwave and shortwave
radiation scheme
3. Noah landsurface processes scheme
4. NMRF PBL scheme and vertical
diffusionTable 5. Configuration of the GRAPES deterministic operational prediction system
GRAPES_GEPS Mesoscale GRAPES_REPS Model area Global; model top at 3 hPa ${10^ \circ } \!\! \!\! {65^ \circ }{\rm N}, {70^ \circ }\!\! \!\! {145^ \circ }{\rm E}$; 10 hPa Forecast time 0000, 0600, 1200, 1800 UTC 0000, 0600, 1200, 1800 UTC Forecast length 360 h 84 h Assimilation system Upscaling of global 4DVar Downscaling of global analysis + cloud analysis Initial perturbation Singular vector ETKF Model perturbation SPPT and SKEB SPPT Lateral boundary condition GRAPES_GEPS Number of ensemble members 30 15 Table 6. Configuration of the GRAPES ensemble prediction system
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 rainyarea satellite data assimilation, convectionpermitting scale ensemble system, kilometerscale data assimilation techniques, land surface data assimilation, etc. These important issues need to be addressed in the future.
Figure 1 shows the yearonyear 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 3km 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 3h rainstorm forecasts produced by the Shanghai Meteorological Service (SMSWARMS 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_Meso3 km was significantly advantageous in the 36h rainstorm forecasts, especially in the first 20h forecast range.
Figure 1. Monthlyaveraged threat score (TS) and its linear trend of 24h rainfall forecast by the GRAPES_MESO since its operational application in July 2006. TS is calculated by using 2510 raingauge observations.
Figure 2. ETS of 3h heavy rainfall forecast averaged over 10 June–18 August 2019 by GRAPES_Meso 3km version (red), in comparison with that produced by the Shanghai Meteorological Service WRF model (SMSWARMS; blue) and that by the ECMWF model (green). This figure is drawn based on the data from the CMA internal network (http://idata.cma/areaHighResolution). The ETS is calculated based on 11,000 Automatic Weather Station (AWS) observations over China. The forecast of each model starts from 1200 UTC.
Figure 3 shows temporal evolution of the 5day 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.
Figure 3. Temporal evolution of 5day hindcast mean average error (MAE) in (a) track (TRK, km), (b) minimum central pressure (p_{min}, hPa), and (c) maximum wind (v_{max}, m s^{−1}) of all the typhoons over the western Pacific during 2016–2018 by the 2019 operational version (V3.0) of the GRAPES_TYM (blue line), the NCEP GFS (green line), and the ECMWF IFS (red line).
It is known that the global forecast system is the core of any integrated NWP system. The time series of the anomaly correlation coefficient (ACC) of 500hPa 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.
Figure 4. Time series of the anomaly correlation coefficient (ACC) for the Northern Hemisphere 500hPa geopotential height forecast on day 5 by GRAPES_GFS (red line), ECMWF (black line), NCEP (green line), and T series (blue line) since 2010.
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 largescale warm lowlevel wind shear, which occurred over the middle and lower reaches of the Yangtze River. For example, the 36h 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 lowlevel 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 36h forecasts and raingauge 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.
Figure 5. Scatter plots of the 12–36h accumulated precipitation (larger than 25 mm) forecasts and the rainfall observations at 2400 raingauge stations over the Chinese mainland from (a) GRAPES_GFS and (b) ECMWF. The model initial time is 2000 Beijing Time 24 May 2019. Unit: mm.
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 largescale heavy precipitation a week ahead and in producing quantitative and fixedpoint probability forecasts for strong weather conditions over local areas.

During the initial design and development of the GRAPES dynamic core, the classic twotimelevel 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 twotimelevel 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 semiimplicit 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 twotimelevel 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 semiimplicit weight close to 0.5, and the time integration reaches nearly twoorder accuracy. This algorithm reconstruction improves the accuracy of the GRAPES dynamic core, and enhances the computational efficiency and stability significantly (Su et al., 2018).
Another issue facing all the semiLagrangian models is the positivedefinite and conservation issue in dealing with the scalar advection (e.g., water vapor). The semiLagrangian algorithm produces a forecast for the subsequent time step by searching for departure points and interpolating the physical quantities at those points. The interpolation process will inevitably result in nonconservation and dissipation, to a certain extent. Bermejo and Staniforth (1992) proposed a quasimonotone semiLangrangian algorithm (QMSL). This algorithm is simple and easy to be implemented, but has relatively low computational accuracy for regions with strong gradients. In addition, this algorithm is also nonconservative. Later, Bermejo and Conde (2002) proposed a corrected conservative algorithm based on the QMSL. Due to the simplicity and high computational efficiency, this algorithm is extensively used in semiLagrangian models, and the ECMWF SISL model with this algorithm remains in use today (Diamantakis and Flemming, 2014). A number of researchers have also developed highly accurate and conservative “volumeremapping”type semiLagrangian algorithms, e.g., the cellintegrated semiLagrangian advection scheme (CSLAM; Lauritzen et al., 2010) and the semiLagrangian inherently conserving and efficient scheme (SLICE; Zerroukat et al., 2004, 2009). However, due to their relatively high computational cost, these algorithms are relatively infrequently used in operational NWP models. In view of this situation, a scalar advection scheme was developed based on the PRM. By rewriting the Lagrangianform prognostic equations of scalar variables to a flux form, a conservative advection scheme is easily constructed, while the fluxes at the cell interfaces are updated by the semiLagrangian algorithm. The PRM allows highly accurate, conservative, and highly efficient water vapor advection calculations (Su et al., 2013).

Since the ECMWF successfully realized operational application of 4DVar in 1997, 4DVar has been regarded 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 4DVar, developing 4DVar is a challenging task. The development of the GRAPES global 4DVar system began in 2008, and its operational application was achieved 10 years later, in July 2018. Based on the global nonhydrostatic GRAPES model, the nonhydrostatic 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 4DVar. As mentioned previously, the GRAPES global 4DVar uses an incremental analysis scheme, in which digital filtering is introduced into the cost function as a weak constraint to suppress highfrequency gravity waves. Multiple outer loops have been designed in the system (a single outer loop is currently used in operations), and the Lanczosconjugate gradient (CG) algorithm is used as the minimization algorithm (Liu et al., 2018). The 4DVar framework adopts a horizontally and vertically inseparable BC model. A secondorder autoregressive 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 nonequilibrium divergent winds and mass fields, is taken into account and achieved by a dynamic and statisticalcombined scheme.
The GRAPES 4DVar system runs 4 times per day, and its initial time of the time window is 0300, 0900, 1500, and 2100 UTC, respectively. The 6h forecast at each initial time is used as the background field of next 4DVar cycle, whereas the 3h forecast within the cycle is the initial condition of 10day forecast at 0600, 1200, 1800, and 0000 UTC, respectively. Multiple types of observations can be assimilated in the GRAPES global 4DVar 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 4DVar system is approximately 50% more than that used in the 3DVar system.
If we take the analysis of ECMWF as a standard to see the quality of GRAPES assimilation, the GRAPES global 4DVar system comprehensively surpasses the 3DVar system, particularly in mass and wind fields. In terms of meansquare error, the background field (6h forecast field) of the GRAPES global 4DVar system is even better than the analysis of the global 3DVar. There is a more pronounced improvement in the geopotential height in the Southern Hemisphere and highlevel wind field in the tropics. Moreover, the GRAPES global 4DVar system outperforms the 3DVar system in terms of various forecast validity periods and variables. As demonstrated in Fig. 6, from the perspective of the ACC for 500hPa 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 4DVar system has a significantly greater ability to forecast moderate and abovemoderate 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 5day forecasts.
Figure 6. The 10day evolution of the anomaly correlation coefficient (ACC) of the 500hPa geopotential height forecast for the (a) Northern Hemisphere and (b) Southern Hemisphere. The dashed and solid curves show the GRAPES 3DVar and 4DVar results, respectively, which are the averages for the period from June 2016 to May 2017. The ACC difference for 4DVar and 3DVar and the 95% confidence threshold are also plotted in the bottom panels. The horizontal axis is the forecast time in hour.

Since the NCEP and the ECMWF successfully assimilated satellite radiance data of TOVS (Television and Infrared Observation Satellite Operational Vertical Sounder) in their operational 3DVar 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 NOAA16 AMSU (Advanced Microwave Sounding Unit)A radiance was investigated in the early GRAPES isobaricsurface 3DVar 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 polarorbiting 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 (Infrared Atmospheric Sounding Interferometer), GOES sounder, AMSUA on AQUA and AMSUA, 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 multiplatform 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 yearonyear 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).
Figure 7. Changes of data types and platform numbers of the satellite data assimilated by GRAPES_GFS since 2010. The right most color bar shows the status at the advanced NWP center of ECMWF. GPS: global positioning system; GEO: geostationary (satellite); LEO: low earth orbit (satellite); and AMVs: atmospheric motion vectors.
The abovementioned 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 FY2 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., 2017, 2018). The Visible and Infrared Spin Scan Radiometer (VISSR) is a payload on board of FY2 satellites. VISSR clearsky watervaporchannel radiance data have been assimilated in the GRAPES global 4DVar 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 FY3 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 FY4A satellite, a newgeneration 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 geostationaryorbit 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 4DVar. In addition, by taking full advantage of the flexible observation capability of the detection instruments on board FY4A as well as the cloud forecasting and the sensitive area detection techniques of the NWP system, an “intelligent rapid clearsky observation mode” and a “highimpact 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 SVbased initial perturbation approach for Typhoons Maria, Ampil, and Mangkhut. Through collaborative interaction of observations and forecasts, a forecastobjectoriented geostationary satellite detection mode was realized, and these intensified observations were assimilated into the GRAPES global 4DVar immediately. These objectobservation experiments obviously improved the typhoon track forecasts compared with the operational GRAPES_GFS for the three typhoon cases.

Developing kilometerscale data assimilation technique is of great importance for improving the forecasts of extreme weather events and nearsurface meteorological parameters (Gustafsson et al., 2018). A GRAPES kilometerscale data assimilation system was developed based on the GRAPES global and regional unified variational assimilation framework (Xue and Chen, 2008).
Several modifications have been made to reconstruct the original regional 3DVar system for developing the kilometerscale NWP (Wang R. C. et al., 2018). 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 10km resolution 3DVar and the kilometerscale 3DVar with 3km resolution. Here, only a singlepoint temperature observation was assimilated. Obviously, the analysis increment of the kilometerscale 3DVar system exhibits more local distribution. Second, a multiscale analysis scheme was implemented for keeping the analysis quality for both the synopticscale and mesotosmallscale weathers. The multiscale analysis was realized by using a multiscale horizontal correlation model (Wu et al., 2018), and by merging the synopticscale and mesotosmallscale components of analysis increments from the global analysis and the kilometerscale analysis, respectively. Also, a weak constraint was introduced into the cost function to maintain the synopticscale information (Yang et al., 2019). In addition, to capture the short life time and quick evolution of mesoscale systems, the kilometerscale 3DVar 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.
Figure 8. Vertical crosssections of temperature analysis increment forced by singlepoint sounding temperature observations at 850 hPa at location (31.93°N, 118.9°E) by (a) 10km and (b) 3km variational assimilation systems of GRAPES_Meso.
Assimilation of observations with high spatiotemporal density is critical for the kilometerscale NWP. In the kilometerscale 3DVar 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 FY4A AGRI (Advanced Geosynchronous Radiation Imager) and GIIRS data is also developed. Furthermore, as significant enhancements to the kilometerscale 3DVar 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 3DVar 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 1month 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 kilometerscale 3DVar can effectively improve the shortrange precipitation forecast for 0–12 h, indicating the great potential in the kilometerscale rapid refresh applications.
Figure 9. ETS scores of precipitation with different intensity for (a) 0–6h and (b) 6–12h forecasts. The black and red bars show ETS scores using initial values from downscaling of global analysis and from 3 kmscale 3DVar combined with cloud analysis, respectively. The ETS scores are averages based on 2400 raingauge observations over China. The red line box represents the 95% confidence interval.

The GRAPES global ensemble prediction system (GRAPES_GEPS) adopted the SVbased initial perturbation approach. The SVs are a set of orthogonal perturbations 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
$\left({u,\;v} \right)$ , the perturbation potential temperature$\left({\theta '} \right)$ , 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 (Li and Liu, 2019; Li et al., 2019). The GRAPES_GEPS was put into operation in 2018 (Li and Liu, 2019). 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 500hPa 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.Northern Hemisphere Southern Hemisphere Ensemble mean Control forecast Improvement Ensemble mean Control forecast Improvement 8.620 7.557 1.062 7.459 6.936 0.523 Table 7. Period (day) of usable forecasts given by the GRAPES_GEPS operational system for the period January–May 2019
The ETKFbased initial perturbations are applied in the GRAPES regional ensemble system (GRAPES_REPS; Zhang et al., 2014). Based on this, a multiscale blending techniques are further developed, which merges the smalltomesoscale component of ETKF with the largescale component from the GRAPES_GEPS (Zhang et al., 2015). A threedimensional bias correction scheme has been proposed for postprocessing the GRAPES_REPS products (Wang J. Z. et al., 2018). In 2019, a GRAPES_GEPSdriven 10km regional ensemble prediction system 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 firstorder 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.
Dynamic core  Physical process  
GRAPES 1.0  GRAPES operational system  GRAPES 1.0  GRAPES operational system  
1. Nonhydrostatic, fully compressible equations (shallow atmospheric approximation)
2. Isothermal reference atmosphere 3. Conventional twotimelevel SISL algorithm 4. 3D vector discretization of the momentum equations 5. Solving the potential temperature prediction equation using the conventional SISL algorithm 6. Quasimonotone positivedefinite scalar advection [quasimonotone semiLagrangian (QMSL)] 
1. Nonhydrostatic, fully compressible equations (shallow atmospheric approximation)
2. 3D reference atmosphere 3. Predictor–corrector SISL algorithm 4. 3D vector discretization of the momentum equations + W damping 5. Solving the potential temperature prediction equation using the vertical noninterpolation SISL algorithm 6. PRMbased highprecision conservative scalar advection 
1. NCEP microphysical scheme with 3–5 water substances
2. Betts–Miller cumulus parameterization scheme 3. Dudhia shortwave and RRTM longwave radiation schemes 4. Mediumrange forecast (MRF) PBL scheme 5. SLAB landsurface processes 6. Dynamics–physics feedback: linear interpolation 
1. Doublemoment microphysics or WRF single moment6 (WSM6) (for mesoscale model)
2. Cloud macrophysics and prognostic cloud scheme (for global model) 3. New simplified Arakawa–Schubert (NSAS) cumulus parameterization scheme or mesoSAS (for mesoscale model) 4. RRTMG longwave and shortwave radiation scheme 5. CoLM (for global model) or Noah (for mesoscale model) landsurface process scheme 6. New MRF (NMRF) or Yonsei University (YSU) (for mesoscale model) PBL scheme 7. Subgrid orographic gravity wave drag 8. Smallscale orographic turbulence form drag 9. Dynamics–physics feedback: cubic spline and quasimonotonic cubic spline interpolation 