Spin-up characteristics with three types of initial fields and the restart effects on forecast accuracy in the GRAPES global forecast system

The spin-up refers to the dynamic and thermal adjustments made at the initial stage of numerical integration in order to reach a statistical equilibrium state. The analyses on the characteristics and effects of spin-ups are of great significance for optimizing the initial field of the model and improving its forecast skills. In this paper, three different initial fields are used in the experiments: the analysis field of four-dimensional variational (4D-VAR) assimilation, the 3 h prediction field in the operational forecasting system, and the Final (FNL) Operational Global Analysis data provided by National Centers for Environmental Prediction (NCEP). Following this, the characteristics of spin-ups in the version 2.3.1 of GRAPES (Global/Regional Assimilation and Prediction System) global forecast system (GRAPES_GFS2.3.1) under different initial fields are compared and analyzed. In addition, the influence of the lost cloud-field information on the spin-up and forecast results of the GRAPES model in the current operation is discussed as well. The results are as follows. With any initial field, the spin-up of GRAPES_GFS2.3.1 has to go through two stages – the dramatic adjustment in the first half-hour of integration and the slow dynamic and thermal adjustments afterwards. The spin-up in GRAPES_GFS2.3.1 lasts for at least 6 h, and the adjustment is gradually completed from lower to upper layers in the model. Therefore, in the evaluation of the GRAPES_GFS2.3.1, the forecast results in the first 6 h should be avoided, and the GRAPES_GFS2.3.1 with its own analysis field performs better than the one using FNL reanalysis data for the cold start in the spin-up because the variations in amplitude of the temperature and humidity tendency are smaller and the spin-up time is slightly shorter. Based on the 4D-VAR assimilation analysis field, the forecast in the operational model is artificially interrupted and restarted after 3 h of integration. In this process, as the cloud-field information is not retained, the spin-up should repeat in the model. The characteristics of spin-up are mostly consistent with those using the 4D-VAR assimilation analysis field as the initial field. However, as the cloud-field information is not retained in the current operation, the hydrometeor content in the atmosphere at the early stage of the forecast is underestimated, affecting the calculation accuracy of the radiation and causing a systematic positive bias of temperature and geopotential height fields at 500 hPa. In addition, the precipitation is also underestimated at the early stage of the simulation, affecting the forecast of typhoon tracks.

1 Introduction

Norwegian scholar Bjerknes (1904) first explicitly proposed the theory of numerical forecasting in the early 20th century. After more than a century of development, it has become an effective way for studying climate change and its causes, as well as forecasting climate and weather. In addition, higher requirements have been also raised for the improvement of numerical forecast accuracy (Bauer et al., 2015; IPCC, 2013).

The numerical forecast accuracy is determined by a variety of factors. The European Centre for Medium-Range Weather Forecasts (ECMWF) concluded that the steady improvement of the numerical forecast in the past 30 years can be mainly attributed to the improvement of the forecast model itself, the application of more observation data, and the development of data assimilation technology (Linus and Erland, 2013). Among them, the performance of the forecast model is mainly determined by the model resolution, the accuracy of finite difference methods, and the representativeness of the physical process parameterization schemes. Observation data mainly depends on the development of monitoring technology, especially the application of satellite data. Data assimilation integrates observation data from different sources with model forecast elements so that the observation data can be comprehensively used by the models. The main purpose of data assimilation is to create a simulated atmosphere state closer to the real atmosphere, reduce the bias of the initial atmosphere condition, and thereby improve the quality of the initial field. In data assimilation, observation data from many sources are used. The uncertainties in the observation data, the inconsistencies among observation elements, and the model flaws (caused by model dynamic assumptions, interactions between physical processes, static data initialization and the radiation balance adjustment, etc.) can lead to inconsistencies between the assimilated new observation input data and the original data in the model. Therefore, the model needs to readjust the dynamic and thermal processes at the initial stage of integration until a new statistical equilibrium state is reached. This process is called the spin-up in numerical modeling, and the time required to reach a new equilibrium state is called the spin-up time (Wolcott and Warner, 1981; Kasahara et al., 1992; Séférian et al., 2016; Sheng et al., 2006; Liu et al., 2008; Xue et al., 2017). During the dynamic and thermal adjustment in the spin-up, spurious gravity waves can be triggered, causing a rapid increase in the root-mean-square error of the forecast variables in the model and an underestimation of the forecast precipitation (Wehbe et al., 2019; Qian et al., 2003). It leads to unreliable forecast results during the spin-up. Therefore, many studies generally do not consider the forecast results during the spin-up when evaluating the model forecasts (Lo et al., 2008; Kleczke et al., 2014; Xie et al., 2013; Zhao et al., 2012). If the spin-up time is too long in the operational model, it would inevitably affect the forecast accuracy of the model. In addition, the overlong spin-up in the climate model or the ocean model can consume excessive computing resources (Duben et al., 2014). Therefore, studying the spin-up characteristics and reducing the spin-up time are of great significance for improving the model forecast and saving computing resources.

Due to different types and usages of numerical models, the spin-up time in different models is greatly different. For example, in global climate models, glacial models, and ocean circulation models, the spin-ups usually take decades to hundreds of years (Scher and Messori, 2019; Danek et al., 2019; Rimac et al., 2017). However, in a regional climate model or a land surface model, only several weeks to several months are needed (Zhong et al., 2008; Rimac et al., 2017; Senatore et al., 2015; Giorgi and Mearns, 1999; Chen et al., 1997). In addition, the spin-up time is also affected by factors such as the simulation domain, the simulation season, and the circulation intensity (Anthes et al., 1989; Errico et al., 1987). The spin-up time of short-term weather forecast models is relatively short, usually several hours to about a dozen hours (Weiss et al., 2008; Souto et al., 2003; Kasahara et al., 1988). To reduce the impact of overlong spin-up on the accuracy of numerical forecasts, many technical methods have been developed to shorten the spin-up time. For example, the “Distorted Physics”, “Matrix method”, “Jacobian-free Newton–Krylov” are used in marine models (Bryan, 1984; Khatiwala et al., 2005; Knoll and Keyes, 2004), and the cloud analysis method for assimilating unconventional observation data such as satellites and radars is used in the short-term weather forecast model to improve the initial humidity field and cloud field, shorten the spin-up time, and improve the short-term precipitation forecast (Li et al., 2011, 2018; Zhu et al., 2017; Xue et al., 2003, 2017; Zhi et al., 2010; Carlin et al., 2017).

The Global/Regional Assimilation Prediction System (GRAPES) is a numerical weather forecast model independently developed by the China Meteorological Administration (CMA). It has become the core of the national numerical forecast operational system in China. Numerical Weather Prediction Center of CMA has established a deterministic weather forecast model system with a global horizontal grid spacing of 25 km and a national horizontal grid spacing of 3 km (Shen et al., 2017; Zhang et al., 2019; Ma et al., 2018; Chen and Shen, 2006). Hao et al. (2013) used the three-dimensional variational (3D-VAR) system to perform the assimilation and analysis of initial fields in the GRAPES regional model, achieving a good forecast result. The research by Zhu et al. (2017) showed that the cloud analysis method in the GRAPES regional model can effectively shorten the spin-up time. After 1 h integration in the model, the precipitation forecast is very close to the observation, and this has a positive impact on the threat score of precipitation forecast within 12 h. Li et al. (2011) also showed similar findings. The assimilation module of GRAPES global forecast system (GRAPES_GFS) was upgraded from 3D-VAR to a four-dimensional variational (4D-VAR) assimilation system in July 2018. The analysis and forecast ability of a 4D-VAR assimilation system is significantly better than 3D-VAR (Zhang et al., 2019). However, there are still many unknowns to be answered. For example, what are the characteristics of the spin-up at the early stage of integration in GRAPES_GFS after the upgrade? In the research and development of the GRAPES-GFS, the widely-used FNL (Final Operational Global Analysis) reanalysis data provided by NCEP (National Centers for Environmental Prediction) (Kalnay et al., 1996) are usually adopted as the model's initial field to quickly evaluate the effects of modification in dynamic core and physical processes on the model forecast performance because the cold start simulation with FNL consumes less computing resources than that of a cycle assimilation simulation. Another question is what advantages the new 4D-VAR assimilation analysis fields have in spin-up process compared with the cold-start simulation with FNL. In addition, we should note that each forecast result of GRAPES_GFS is from the model integration forecast based on the 4D-VAR assimilation analysis field 3 h ago in the current operational forecast system. For example, the 12:00 UTC forecast result is based on the 4D-VAR assimilation analysis field at 09:00 UTC. Actually, for numerical weather prediction model's users (especially forecasters), they are usually accustomed to referring the forecast productions of model staring to integrate from 00:00 or 12:00 UTC (or more time, for example 18:00 UTC). Thus, considering the habits of users when using the forecast results, GRAPES_GFS integrates for 3 h (to 12:00 UTC) to retain the essential meteorological element fields (U, V, T, Q, H, TS Ps, etc.), and then the integration is terminated and restarts from 12:00 UTC by using the newly saved meteorological field data. The model forecast results thereafter are released, that is, the forecast results at 12:00 UTC are obtained by users. In this process, the cloud-field variables (the mass and concentration of hydrometeors and cloud cover) during the first 3 h of integration are not retained in the model, losing the cloud information formed after the 3 h spin-up. The reasons for the unretained cloud-field variables were mainly based on the following considerations: the hydrometeor contents are in very small amounts relative to water vapor, and they can be quickly created in the spin-up process when the model restarts. Moreover, this treatment can save storage space and input–output (IO) time. However, its impacts on the spin-up process and model forecast performance have not yet been carefully analyzed and evaluated. Therefore, we need to fully diagnose and analyze the necessity of the repetition of GRAPES_GFS spin-up during the reintegration, and the impact of the lost cloud-field information on the later forecast. In this regard, the characteristics of spin-ups in GRAPES_GFS using the 4D-VAR analysis data and the FNL data separately as the initial field are compared and analyzed, and the impacts of the cloud-field information loss in the current operation on the spin-up after the model restart and on later forecast results are discussed. This paper aims to provide the scientific basis for understanding the characteristics of GRAPES_GFS at the initial stage of integration and improving the assimilation system and operational procedure.

The paper is organized as follows. In Sect. 2, the GRAPES_GFS forecasting system and the experiment settings for one case study are introduced. In Sect. 3, the main research results are presented. Finally, in Sect. 4, the main conclusions are given, and some issues about spin-ups are discussed.

2 GRAPES_GFS2.3.1 and experiment setup

2.1 GRAPES_GFS2.3.1

GRAPES is a global numerical weather prediction system that is composed of an atmospheric model and a variational data assimilation system (3D-VAR/4D-VAR). The framework of the atmospheric model is a fully compressible non-hydrostatic dynamical one with semi-implicit and semi-Lagrangian time difference scheme. In the horizontal direction, the equidistant latitude–longitude grid system with the Arakawa-C grid and central differencing of second-order accuracy for variable staggering is used, and in the vertical direction, the height-based terrain-following coordinate with the Charney–Phillips staggering is adopted. Forecast variables of GRAPES_GFS include the dimensionless air pressure (Exner function), potential temperature, three-dimensional wind field components, and specific humidity. It also introduces the Piecewise Rational Method (PRM) scalars (Su et al., 2013) into the model, which is a scheme of water vapor advection. The physical parameterization schemes used in the GRAPES_GFS operation mainly include the long-wave and short-wave radiation schemes (the rapid radiative transfer model, RRTMG) (Morcrette et al., 2008; Pincus et al., 2003), the land surface scheme (the Common Land Model, CoLM) (Dai et al., 2003), the planetary boundary layer scheme (Medium-Range Forecast, MRF) (Hong and Pan, 1996), the deep and shallow cumulus convection parameterization scheme (the New Simplified Arakawa–Schubert, NSAS) (Arakawa and Schubert, 1974; Liu et al., 2015; Pan and Wu, 1995). The cloud physics scheme includes the macro cloud scheme dealing with the condensation process under the unsaturated condition of grid-average water vapor, a double-moment cloud microphysical scheme, and a cloud cover prognostic scheme (Chen et al., 2007; Ma et al., 2018). On 1 July 2018, the GRAPES global 4D-Var data assimilation system came into operation (Zhang et al., 2019), which is called version 2.3.1 of GRAPES_GFS (abbreviated as GRAPES_GFS2.3.1). The GRAPES_GFS2.3.1 version is adopted in this research.

2.2 Experiment setup

In this paper, GRAPES_GFS2.3.1, with the operational forecast time of 00:00 UTC on 9 August 2019, is taken as an example, and three experiments are set up to analyze the similarities and differences in the spin-up characteristics of the model using different initial fields. The settings are shown in Table 1. In the first experiment, the analysis field provided by the 4D-VAR assimilation analysis system in the operational forecast at 21:00 UTC on 8 August 2019 is used as the initial field to directly perform model integration forecasts, and the initial time is 21:00 UTC on 8 August. This experiment is called G21. For the second experiment, called G00, its initial field adopts 3 h integration output of G21 without retaining cloud-field information. That is to say, at 00:00 UTC on 9 August, it retains the G21's 3 h forecast variables (u and v wind field components, potential temperature, water vapor, and dimensionless air pressure, etc.) required by the pre-processing system and stops the integration. During the process, the fields of all hydrometeor contents and cloud cover are lost considering the limitation of IO time and disk space. Then the model restarts at 00:00 UTC on 9 August with the reserved forecast-field information for forecasting in G00. Moreover, the model output of G00 is exactly the forecast results to be provided to users in the GRAPES_GFS2.3.1 operation. The third experiment uses the initial field from the NCEP FNL reanalysis data at 00:00 UTC on 9 August 2019 to perform the integration forecast. The purpose is to compare the spin-up characteristics of GRAPES_GFS2.3.1 model, respectively, using its own analysis field and FNL reanalysis field as the initial field. This experiment is called F00. To analyze the impacts of the initial field on the forecast, G00 and F00 produce a continuous 72 h forecast. As G21 starts the integration 3 h earlier than the other two, the forecast of G21 lasts for 75 h to ensure the same forecast and analysis period with G21 and G00.

Table 1 Model setup of three experiments used in this study.

All three experiments are based on the GRAPES_GFS2.3.1 operational model, with a horizontal grid spacing of 0.25^∘, 60 vertical layers, and a model integration time step of 300 s. The physical schemes used are from the operational setup (as described in Sect. 2.1), and the assimilation module is 4D-VAR assimilation system. To explicitly analyze the spin-up characteristics of the GRAPES_GFS2.3.1 at this early stage of integration, the results of each integration step are output, and the temperature tendency (TT) and water vapor tendency (WVT) fields at each model layer during the dynamic and physical processes are retained.

In addition, the cloud-field information has not been saved during the restart in the current operation. To examine its impact on the accuracy of the later forecast, this study investigates the super typhoon “Lekima” (no. 1909) that landed in China during the selected forecast period, and the forecast differences in cloud, precipitation field, and typhoon track during Lekima between G00 and G21 are analyzed.

3 Results

3.1 Characteristics of spin-ups

3.1.1 Characteristics of total WVT and total TT

To analyze the spin-up characteristics of GRAPES_GFS2.3.1, the initial fields in F00, G21, and G00 are used to perform the integration, and the temporal variations of the average total WVT and TT at different heights from 00:00 to 12:00 UTC are calculated, as shown in Fig. 1. Seen from the figure, both the WVT and TT show sharp fluctuations at the initial stage of the integration in the three experiments, especially during the first hour. After 3–6 h of spin-up adjustment, the variation magnitudes of WVT and TT gradually become gentle, but the variation characteristics vary with different initial fields. At the early stage of the integration, the WVT is adjusted in F00 and G21, with the amplitude of −4.5 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$. In G21, the water vapor adjustment occurs in the lower layers of the model (850 and 925 hPa), while the WVT is relatively gentle without an obvious adjustment in the upper and middle layers (500 and 300 hPa). In F00, the water vapor adjustment occurs at the upper levels of the model at the early stage of integration. The WVT at 300 hPa can reach −4.5 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$, but it weakens immediately afterwards, probably due to the supersaturated water vapor in the initial field from FNL data. In F00, the WVT in the lower layers of the model is also significantly larger than that in G21. For example, at 850 hPa, the WVT in F00 maintains about 1 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ for a relatively long time but in G21 mostly changes within 0.5 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$. The corresponding temperature adjustment processes in the two experiments present the same variation characteristics as the WVT adjustment. Therefore, the spin-up in the integration using the analysis field of GRAPES_GFS2.3.1 as the initial field is gentler than that using the FNL reanalysis data as the initial field.

Figure 1 Time evolution of global mean of the total water vapor tendency (WVT) and total temperature tendency (TT) at different vertical levels from 0 to 12 h simulated by F00 (a, b), G21 (c, d), and G00 (e, f) experiments. The unit of WVT and TT is $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ and K d⁻¹, respectively.

In G21 and G00, the variations of both WVT and TT are very consistent, indicating that G00 has inherited the temperature and humidity structure of G21 well. However, G00 still needs to go through the spin-up during which a gradually stable adjustment process follows a sharp fluctuation at the early stage of integration; i.e., the dynamic and thermal adjustments are required to reach a statistical equilibrium state in the model. At the initial stage of integration in G00, the variation amplitudes of WVT and TT are smaller than those in G21, but greater than those in G21 after the 3 h integration. It shows that although G00 can retain the temperature and humidity structure of G21, the loss of cloud-field information in the operation still has a destructive effect on the model equilibrium state after 3 h adjustments. Based on the variation of TT, the spin-up time required for G00 is generally less than that for G21. It takes about 6 to 8 h to reach a TT equilibrium state in G21, but it is less than 6 h in G00.

3.1.2 Tendency characteristics of the model dynamical and physical processes

In GRAPES_GFS2.3.1, the total temperature tendency of the model (ALL) is determined by dynamic core (DYN), radiation process (RAD), turbulent mixing in planetary boundary layer process (PBL), cumulus convection process (CONV) and cloud physical process (CLOUD). Among them, the total temperature tendency of all physical processes (PHY) is defined as the sum of the last four items (PHY = RAD + PBL + CONV + CLOUD). Likewise, the total water vapor tendency for ALL and PHY are same to those of temperature tendency except for the radiation process (RAD). Figure 2 shows the temporal variation of mean WVT due to dynamic and physical processes at different heights in F00, G21 and G00. In the middle and upper layers of the model (Fig. 2a and d), there is a drastic adjustment in the atmosphere at the early stage of the integration in F00. It may be due to the supersaturated water vapor in the initial field from FNL data, which causes the cloud to condense very quickly, and thus a relatively stable state is reached after three integration steps. At these levels in G21 (Fig. 2b and e), the total WVTs at the first few integration steps are slightly larger than those at the subsequent integration steps. The variations of the WVTs from dynamic core and turbulent mixing process in the planetary boundary layer are much less than those from the cumulus convection process and cloud physical process, and the latter two processes jointly determined the variation of WVTs at 300 and 500 hPa. There is not much difference in the dynamic field tendencies between G21 and F00. The magnitudes of the WVTs in the dynamic processes of the two experiments are also very close: around 0.5 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ at 500 hPa and 0.25 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ at 300 hPa. Therefore, the differences of the upper middle-level water vapor adjustments in the spin-ups between G21 and F00 are mainly caused by physical processes, and there is a good consistency in the dynamic process between the two experiments. At 925 hPa (the lower layer of the model), the total WVT stays around 1 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$ in F00 after three integration steps, humidifying the atmosphere. In G21, it reaches a relatively stable state after six integration steps, and water vapor decreases overall. As the WVTs of the dynamical processes in F00 and G21 have the same magnitude around 0.25 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$, the difference of the total WVT between G21 and F00 is mainly caused by physical processes. The effect of the boundary layer on the WVT is similar in both experiments, and the WVT is almost 3 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$. The greatest difference between the two experiments is mainly caused by the convection scheme. The convection in F00 is relatively gentle, and the WVT from convection is around −1 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$. In contrast, due to the strong dehumidification ability of convection in G21, the WVT is between −5 and −2.5 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$, which is significantly stronger than that in F00. At 925 hPa, the water vapor mainly decreases due to the strong convection process in G21. Such a significant difference in the convection processes between F00 and G21 may be related to the low-level temperature and humidity structures and the triggering conditions for convection. Meanwhile, it can be seen that the difference in the initial field of the model can significantly affect the physical processes.

Figure 2 Time evolution of mean water vapor tendency (WVT) of the dynamical core and each physical process at 300, 500, and 925 hPa heights from 0 to 1 h simulated by the F00 (a, d, g), G21 (b, e, h), and G00 (c, f, i) experiments (values given in $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$).

In summary, in the middle and upper atmosphere, the fluctuation of WVT in G21 is weaker than that in F00, indicating the advantage of using the data assimilation cycling as the initial field. Both experiments quickly reach a quasi-equilibrium state after dramatic adjustments over several integration steps. The water vapor adjustment in spin-ups mainly occurs in the lower atmosphere of the model. The difference is mainly caused by different convection schemes. At the same time, different initial fields of the temperature and humidity structure may lead to a great difference in the dehumidification ability of convection. For G00 and G21, the WVTs of the dynamic and physical processes have roughly the same characteristics. At all of the three levels, the WVTs in G00 are slightly lower than those in G21.

Figure 3 The same as Fig. 2 but for the results of temperature tendency. (values given in K d⁻¹).

In the middle and upper layers of the model, the dramatic change of the TT in F00 mainly occurs within the first half-hour of the integration (Fig. 3a and d). Among all the TTs at the first integration step, the cloud physical process leads to the largest one, followed by convection process, and they are related to the water vapor condensation process (Fig. 2a and d). For example, at 500 hPa, the global average heating produced by the cloud microphysical condensation process at the initial time can exceed 5 K d⁻¹, and it takes four integration steps to reach a relatively stable state. However, at this level, the TT caused by the convection process is 3 K d⁻¹, and it only needs one integration step with the drastic adjustment to get relatively stable. In addition, the TT caused by the dynamic process fluctuates greatly at the first half-hour of the integration. For example, at 300 hPa, the TT fluctuates between 1.1 and 1.5 K d⁻¹, and it requires extra 3 or 4 integration steps to reach a relatively stable state compared to the physical processes. Nevertheless, after half an hour of severe fluctuations, the TT caused by dynamic and physical processes tends to be relatively stable. Overall, the temperature increases by 0.25 to 0.5 K d⁻¹ in the middle and upper atmosphere in F00. Compared with that in middle and upper layers, the TT variation caused by the dynamic and physical processes in the lower layer of the model (Fig. 3g) shows a relatively small and rapid adjustment at the first integration step. However, no drastic adjustment is shown afterwards, and its variation is relatively stable. The TT of the convection process at 925 hPa in F00 varies between 1.5 and 2 K d⁻¹, which is mainly caused by condensing and dehumidifying of the atmosphere (Fig. 2g). Except for the cloud physical process, which has a relatively small positive tendency in the first four time steps, the TTs of dynamic core and other physical processes are all negative. Overall, in F00 the total atmospheric temperature is reduced with an amplitude of about −1.2 K d⁻¹ in the first hour of the integration at 925 hPa.

In G21, the TT in the middle and upper layers also experiences a dramatic adjustment in the first half-hour of the integration (Fig. 3b and e), and the main reason for the fluctuation is the dehumidification and heating in the convection process, which is different from that in F00 caused by the cloud physical process. The temperature increase caused by the convection process in G21 is 1 to 2.5 K d⁻¹, which is about twice that in F00. The TT caused by the cloud physical process in G21 varies relatively gently. Similar to F00, the TT caused by the dynamic process in G21 also shows obvious fluctuations, which may be caused by the drastic variations of physical processes. In the lower layer of 925 hPa (Fig. 3h), the positive TT in G21 is also caused by convective dehumidification and heating, while other processes lead to cooling. In terms of the total TT (dynamic core plus all physical processes), F00 has a cooling effect with a value of −1 K d⁻¹, while G21 has a warming effect with a value within 1 K d⁻¹. The temperature increase rate of G21 gradually decreases with the integration step.

The characteristics of the TT variation in G00 are consistent with those in G21 (Fig. 3c, f and i). In the first few time steps, G00 also has an adjustment process, with the adjustment amplitudes of TT close to half those in G21 at all levels. After half an hour, the temperature tends to be relatively stable. The TT variation in G00 indicates that although G21 has undergone a 3 h spin-up, G00 needs to undergo it again due to the loss of cloud-field information during the restart, and its fluctuation amplitude is not substantially smaller than that of G21.

3.1.3 Evolution characteristics of the cloud field

The comprehensive adjustment effect of the dynamic and the physical processes on the water vapor and temperature in the numerical model can be presented by the cloud state. To reveal the dynamic and thermal adjustment processes in GRAPES_GFS2.3.1 system at the beginning of the integration and the time required for the model to reach the statistical equilibrium state (spin-up time), this section uses the total grid number of cloud (TGNC) in the model as the index for analyses. Although the cloud is changing locally, the total area covered by cloud can be regarded as a constant globally on average. Therefore, TGNC is used as the analysis index, and the model is considered to have completed the spin-up when the TGNC gets relatively stable. The total hydrometeors content (THC, THC = cloud water + raindrop + cloud ice + snow + graupel) greater than 1.0 × 10⁻⁴ g kg⁻¹ in GRAPES_GFS2.3.1 is defined as the grid with cloud, and the TGNC at a global scale or a certain height is the sum of all the grids in the corresponding cloud area.

Figure 4 Vertical distribution of total number of cloud points at different forecast time simulated by F00 (a), G21 (b), and G00 (c) experiment, respectively (values given in number ⋅ 10 000).

Figure 5 Distributions of all hydrometeor content at 400 hPa at different forecast times (5 min, 1, 3, 6 h) simulated by F00 (a–d), G21 (e–h), and G00 (i–l) experiments, respectively (values given in g kg⁻¹).

Figure 4 shows the vertical distributions of TGNC at different lead times in three experiments. It can be seen that regardless of whether the GRAPES_GFS2.3.1 model is cold-started with reanalysis data (F00, Fig. 4a) or warm-started with the 4D-VAR analysis field as the initial field (G21, Fig. 4b), the TGNC experiences rapid generation and growth during the 3 h after the beginning of integration in the two experiments, especially in the middle- and low-cloud regions below 300 hPa. After 3 h of integration, the TGNC grows relatively slowly, while after 6 h of integration, the TGNC becomes basically stable. However, the time required for the TGNC to reach the equilibrium state is slightly different at different heights. In F00, the integration time required for the TGNC to gradually reach the statistical equilibrium state below 850 hPa is 6 h. Note that the statistical equilibrium state is defined when the difference of TGNC with respect to the 24 h integration is insignificant (the difference is less than 20% of TGNC at 24 h). However, it takes 6–12 h for the TGNC to get stable and it completes the spin-up above 850 hPa. For G21, the TGNC of the middle and low cloud below 300 hPa needs 6 h to reach the statistical equilibrium state, while the TGNC of the high cloud above 300 hPa needs 6–12 h. It can be seen that the GRAPES_GFS2.3.1 using the analysis field from its own data assimilation cycling enables the cloud field in middle and upper layers to reach the equilibrium state earlier than that using FNL data for the cold start. In addition, GRAPES_GFS2.3.1 is gradually adjusted from the lower to the upper layers of the model to reach the equilibrium state, which is consistent with the evolution characteristics of the thermodynamic process in the troposphere. For the cloud above 500 hPa, the TGNC in F00 is significantly more than that in G21, which is related to a higher relative humidity of the initial field. Compared with G21, F00 has a wetter water vapor environment at the upper levels (Fig.6d), which tends to quickly condense the water vapor into more hydrometeors through the cloud scheme to eliminate supersaturated water vapor at the beginning of the integration (Fig. 2a). Thus, F00 has a higher hydrometeor content value and a wider distribution of cloud region (Fig. 5a and e), and its TGNCs are also larger than those of G21 at the upper layers.

In G00 (Fig. 4c), the growth of TGNC is found to be much slower than that in G21, especially for the TGNC of the middle and upper cloud. For example, at 3 h after the beginning of G00, the TGNC of the middle cloud is mostly between 15 and 20, while the TGNC in G21 can reach 25–30. The reason may be that the humidity and temperature fields of the model in G21 are already in a relative equilibrium state after 3 h spin-up. Meanwhile, as the restart of GRAPES_GFS2.3.1 has lost the cloud-field information (dotted light blue line) from the 3 h integration, the TGNC cannot reach the previous magnitude in the middle and upper layers even if it has been integrated for 24 h in G00 (Fig. 4c, solid purple line).

Figure 5 shows the distributions of THC at 400 hPa at different forecast time in the three experiments. It can be seen that the temporal variation characteristics of THC and its horizontal distribution at 400 hPa have consistent results with those shown in Fig. 4. In F00 and G21, as supersaturated water vapor is removed from the initial field, the cloud is quickly generated at the first integration step of the model. The THC rapidly increases within 1 h, and the cloud area with high hydrometeor content is constantly expanding. For example, at 1 h into the integration in F00, the THC in most areas of the Pacific Warm Pool is 0.2 g kg⁻¹. With the further adjustment of the spin-up, the THC in this area gradually decreases and maintains a relatively equilibrium state after 6 h of integration. The variation characteristics of the THC in the storm track area (60–30^∘ S) in the Southern Hemisphere are similar to those in the warm pool area but are less significant.

Figure 6 Distribution of water vapor content (WVC) (a, b) simulated by F00 and G21, and the differences of WVC (c) and relative humidity (RH) (d) between F00 and G21 (F00-G21) at 400 hPa in their initial fields. WVC and RH are given in units of g kg⁻¹ and %, respectively.

Experiments using the 4D-VAR analysis field to provide the initial field (Fig. 5e–h) show that the variation characteristics of THC at 400 hPa are generally consistent with those in F00. After the first integration step of GRAPES_GFS2.3.1, cloud areas are quickly generated in tropical and midlatitude areas. Due to the rapid development of convection processes in tropical areas, more cloud with THC of 0.0–0.05 g kg⁻¹ appears. After 3 h of integration, the development of the cloud area gradually weakens. After 6 h of integration, the variations of the range and shape of the cloud area are no longer obvious, and it can be considered that a relatively equilibrium state is reached. From the view of absolute value of THC in the cloud area, although the difference in the distribution range of the cloud is insignificant, the THC in G21 is significantly less than that in F00 due to the different temperature and humidity conditions in their initial fields (Fig. 6).

Since G00 does not retain the cloud-field information after 3 h of integration in G21 (the THC in Fig. 5g), the model needs to undergo a new cloud-generation process when restarting the integration. However, as the dynamic and thermal fields are obtained after 3 h of adjustments in G21, the relative humidity has undergone a condensation process, making the atmosphere of G00 have a much weaker supersaturation at the initial time than that in G21. Therefore, unlike F00 (Fig. 5a) or G21 (Fig. 5e), in which large-scale cloud appears instantaneously, the cloud field in G00 can only be gradually generated by the dynamic and physical processes of the model. It can be seen from Fig. 5i–k that this process is relatively slow, and a relatively stable cloud distribution does not appear until 3 h after the integration. The cloud range in G00 at that time is smaller than that in G21, and it generally reaches the equilibrium state after 6 h of integration. The influence of slower generation and smaller range of the cloud in G00 on the model forecast results will be analyzed and explained in Sect. 3.2.

To reveal the reason why the TGNC (Fig. 4) and the THC (Fig. 5) in the upper layers of the model in F00 are significantly higher than those in G21, the difference of water vapor content and relative humidity at 400 hPa is analyzed, and the results are shown in Fig. 6. Figure 6c shows that the specific humidity in the initial field of F00 is generally higher than that of G21 in the tropical areas and the midlatitude and high-latitude areas of the Northern Hemisphere, especially in the tropical warm pool area where the difference is mostly over 0.2 g kg⁻¹. The relative humidity reflects the degree of water vapor saturation. Figure 6d shows that the humidity of the initial field from the FNL reanalysis data is high relative to that from the 4D-VAR analysis field in the tropical warm pool, Intertropical Convergence Zone (ITCZ), and midlatitude and high-latitude areas at 400 hPa. This means that the water vapor is more likely to get saturated using the FNL reanalysis data as initial field. Thus, the cloud area is larger and the THC is higher at the beginning of the integration. It is not difficult to conclude that there are differences in the structure of atmospheric temperature and humidity among different initial field data, which significantly impacts the spin-up characteristics of the model and the cloud formation and development. It also suggests that we need to pay more attention to the analysis quality of water vapor in data assimilation (DA). It has been also confirmed by previous studies (Weygandt et al., 2002; Ge et al., 2013) that having an accurate moisture initial field by DA is an effective way to improve the forecast performance of supercell storms in numerical weather prediction models.

3.2 Impacts on later forecast results

It can be seen from Sect. 3.1 that the cloud-field information formed in the first 3 h of integration has not been saved operationally; thus, the model must restart the spin-up, and THC appears to be significantly less in the new spin-up. In order to discuss the impact of the restarted spin-up and the decreased THC on the later forecasts by GRAPES_GFS2.3.1, the global radiation field and synoptic field (temperature and geopotential height) are analyzed in this section. The cloud and precipitation fields and the track of the super typhoon Lekima that made landfall in China during the simulation period will be analyzed as well.

3.2.1 Impacts on global radiation

Figure 7 shows the zonal mean distributions of averaged column cloud water content (CCWC), the outgoing longwave (OLR) at the atmosphere top and the downward longwave at ground (GDLW) level simulated by G21 and G00 from 00:00 to 03:00 UTC on 9 August 2019, as well as the distributions of difference between them. It can be seen from Fig. 7a that the total zonal-averaged CCWC forecasted in G00 is systematically smaller than that forecasted in G21. The areas with smaller CCWC are mainly located in the Southern Hemisphere storm track, tropical low-latitude areas, and midlatitude and high-latitude areas in the Northern Hemisphere with active cloud. Among them, the area with the smallest CCWC is the active area of Southern Hemisphere storm track, with the CCWC difference reaching 240 g m⁻², and there are also some areas with the CCWC difference over 200 g m⁻² in the Northern Hemisphere. From the OLR and GDLW predicted in the two experiments, it can be seen that the OLR predicted in G00 is systematically larger than that in G21, with the maximum bias (20 W m² s⁻¹) appearing in the Southern Hemisphere storm track. This is due to the interaction between clouds and radiation, as well as the underestimation of the CCWC. In terms of GDLW, the reduced CCWC weakens the atmospheric warming effect, resulting in systematically smaller GDLW in G00 than in G21. In most areas, the GDLW is smaller than the observation by over 10 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and the regions with the largest bias are the midlatitude and high-latitude areas of the Southern Hemisphere and high-latitude areas of the Northern Hemisphere.

Figure 7 Zonal means and their differences of (a) 3 h-averaged cloud water path (CWP), and (b) the outgoing longwave (OLR) at the top of atmosphere and the downward longwave at ground (GDLW) simulated by G21 and G00 experiments for 00:00–03:00 UTC, 9 August 2019. The units of CWP and OLR/GDLW are g m⁻² and W m⁻², respectively.

3.2.2 Impacts on the global temperature and geopotential height fields

The change in the calculation of the radiation flux induced by cloud would seriously affect the atmospheric temperature field and geopotential height field. Figure 8 shows the difference distributions of the 500 hPa temperature field and the geopotential height field at four lead time between G00 and G21. It can be found that as there is less hydrometeor in the cloud in G00 than in G21, the temperature field in G00 at different forecast times shows a systematic warming of more than 0.1 K in the tropical low-latitude and midlatitude and high-latitude areas with active cloud. With the increase of the lead time, the warming area is expanding and the degree of warming gradually increases. For example, after 72 h of integration, the warming in many areas is larger than 0.2 K, and it can reach 0.5 K in some areas. Systematic biases also appear in the corresponding geopotential height field. Compared with those in G21, the geopotential height fields in G00 have also systematic positive biases. For example, in the first 24 h of integration, the systematic biases in the geopotential height field are above 0.5 gpm, and the positive bias can exceed 1 gpm in areas with active cloud. After 72 h of integration, the geopotential height field in the tropical area still shows a systematic positive bias, while in the midlatitude and high-latitude areas, the bias of the geopotential height field shows the structure with an alternation of positive biases and negative biases due to the biases of the weather system location predicted in the two experiments, but in most areas the forecast fields are still higher than the observation.

Figure 8 Distribution of the differences (G00 minus G21) of temperature field (a–d) and geopotential height field (e–h) at 500 hPa simulated by G00 and G21 experiments. The units of temperature and geopotential height are K and gpm, respectively.

3.2.3 Impacts on typhoon forecasts

This section analyzes the biases of the cloud field, precipitation field, and the track of the super typhoon Lekima (no. 1909) and typhoon “Krosa” (No. 1910) in 2019 during the forecast period to evaluate the impact of the lost hydrometeor information on typhoon forecast operation in GRAPES_GFS2.3.1. During the forecast, Lekima and Krosa appear as double typhoons in the western Pacific. Lekima made landfall in northern China, while Krosa remained offshore. Since the conclusions for both Lekima and Krosa are the same, only Lekima will be presented in this study. Here, we show the impact on the cloud and precipitation of Lekima by the lost hydrometeor information on typhoon forecast operation of GRAPES_GFS2.3.1. In the last part, the path-forecast biases for the two typhoons are both given.

Figure 9 Time evolution of the sum of averaged column cloud water content (CCWC) and column cloud ice content (CCIC) at the typhoon Lekima region (22–34^∘ N, 117–130^∘ E) simulated by G00 and G21 experiment (values are given in g m⁻²).

Figure 9 shows the evolutions of the averaged CCWC and column cloud ice content (CCIC) within the main cloud area of Lekima (22–34^∘ N, 117–130^∘ E) simulated in G00 and G21 from 00:00 UTC on 9 August to 00:00 UTC on 10 August 2019. It can be seen that the CCIC predicted in G00 at the early stage of integration is obviously underestimated. The averaged CCIC values in G21 are maintained within 850–1000 g m⁻² from 00:00 to 09:00 UTC on 9 August, while the CCIC is only 480 g m⁻² at the initial time of G00. G00 needs to restart the spin-up. During the spin-up, the CCIC predicted in G00 increases rapidly, with the greatest increase during 00:00 to 06:00 UTC. After 3 h of the integration, the CCIC increases rapidly from 480 to 820 g m⁻². After 6 h of integration, the CCIC is close to 900 g m⁻². In G00, the CCIC is not as large as that in G21 until 9 h after the beginning of integration.

Figure 10 Distribution of the differences (G00 minus G21) of 3-hourly and 24 h accumulated precipitation (since 00:00 UTC 8 August 2019) of the typhoon Lekima simulated by G00 and G21 experiments (values are given in mm).

Figure 10 shows the difference distributions of both 3 and 24 h accumulated precipitation (since 00:00 UTC 8 August 2019) of Lekima between forecasts of G00 and G21.The most significant difference of the 3 h cumulated precipitation appears within the first 3 h of integration in G00. The 00:00–03:00 UTC precipitation forecasted in G00 presents a systematic underestimation when compared with G21, and the biases are all above 1 mm. The precipitation bias in the center of Lekima can even exceed 5 mm (Fig. 10a). As shown in Fig. 9, after 3 h of adjustments, the total CWP and CCIC in the typhoon system in G00 grows rapidly and gets close to the magnitudes in G21. Therefore, the difference of the 3 h precipitation between forecasts of G00 and G21 is no longer significant during 03:00–06:00 and 06:00–09:00 UTC, and there is no more systematic bias (Fig. 10b and c). The phase differences of the weather system lead to the structure with alternating positive biases and negative biases for the precipitation difference.

It can be found from Fig. 10d that the lack of cloud-field information has a significant impact on the simulation of the accumulated precipitation in the first 24 h of Lekima. The negative biases dominate the central area of the typhoon; i.e., there is an underestimation of precipitation with the maximum bias of 5–10 mm. In contrast, in the spiral cloud zone around the typhoon, there is a structure with an alternation of positive and negative biases, which is related to the location bias of the weather system simulated in the two experiments in this area.

Figure 11 Time evolution of the forecasted track errors of G00 and G21 experiments for the typhoons Lekima and Krosa during the forecast period of 72 h (values are given in km).

Figure 11 shows the forecast track evolution of Lekima and Krosa in G00 and G21 within the lead time of 72 h. Overall, G21 performs better than G00 in predicting the tracks of these two typhoons, and there are different characteristics for the track forecast biases of the two different typhoons. Lekima landed on the coast of Chengnan Town, Wenling City, Zhejiang Province, at 15:45 UTC on 9 August 2019. There is not much difference in the biases of the track forecast between G00 and G21 before the Lekima landing. In contrast, the biases appear to be different after the landfall (16:00 UTC), and the track forecast in G21 is slightly better than that in G00 around the landfall. After the landfall, the track biases change continuously during the 27th to 42th hour and 54th to 60th hour of the forecast, the track bias in G21 is smaller than that in G20. The maximum difference between the two track forecasts can reach 32 km. From the 65th to 72th hour, the forecast track bias in G21 is slightly larger. For Krosa, during the first 42 h, the biases of the tracks forecasted in G00 and G21 are not much different. But the forecast tracks of the two become different after the 42th hour, with the track bias in G00 becoming larger. In most forecasts after the 42th hour, the track biases in G00 are over 20 km and larger than those in G21. The abrupt track difference after 42th hour is most likely caused by the continuous accumulation of the direct cloud-radiation process and the systematic temperature bias in the typhoon peripheral cloud area during the re-undergone spin-up of G00 experiment, along with their impacts on the typhoon eye (track) through dynamic processes with the model integration.

Overall, G21 performs better than G00 in the track forecasts of Lekima and Krosa within the lead time of 72 h, especially in the forecast of Krosa. For Krosa, the forecast track on the ocean is less affected by other factors, so the forecast track biases at the later stage of the forecast are significantly smaller. It shows that GRAPES_GFS2.3.1 performs better in continuous-integration forecasts, and the interruption in the operation is destructive to the typhoon track forecast.

4 Conclusions and discussion

To analyze the characteristics of the spin-up at the early stage of integration in GRAPES_GFS2.3.1, this study adopted three different initial fields, namely the 4D-VAR analysis field (G21), the field obtained by interrupting and restarting the 4D-VAR analysis field after 3 h of integration (G00), and the field based on FNL reanalysis data for a cold start (F00). Moreover, the differences between G00 and G21 on the later model forecast results were analyzed to evaluate the impact of current operational procedure on GRAPES_GFS2.3.1 forecasts. The main conclusions are as follows.

All three experiments using different initial fields show that the spin-up of GRAPES_GFS2.3.1 has to go through two stages: the dramatic adjustment in the initial half-hour of integration and the slow dynamic and thermal adjustment afterwards. In the middle and lower layers of the model, the spin-up takes 6 h to reach the equilibrium state and takes longer in the upper layers. The dynamic and thermal adjustment is gradually completed from the lower to the upper layer of the model.

The GRAPES_GFS2.3.1 using its own analysis field as the initial field (G21) is gentler in the water vapor and temperature adjustment in the spin-up than the GRAPES_GFS2.3.1 using FNL reanalysis data for cold start (F00), and the time required is slightly shorter. Due to the different structures of temperature and humidity in the two initial fields, the differences of physical processes in the model spin-up adjustment are obvious, especially regarding the convections and cloud physical processes. However, the differences in dynamic processes are not obvious. G00 needs to repeat the spin-up. Its dynamic and thermal adjustments are similar to that in G21. The temperature and humidity adjustment in G00 is slightly weaker than that in G21, and its spin-up is slightly shorter.

In G00, the cloud-field information is not retained during the current operation of GRAPES_GFS2.3.1. It shows that G00 significantly underestimates the atmospheric CCWC and CCIC at the early stage of forecast, which would affect the calculation accuracy of radiation and result in systematic positive biases in temperature and geopotential height fields at 500 hPa. Due to the lack of cloud-field information, the accumulated precipitation in the first 3 h of integration in G00 is significantly underestimated. The 24 h accumulated precipitation in the typhoon center is also less than that in G21, and a destructive effect is made on the typhoon track forecast.

Regarding the influence of the lost cloud-field information in the GRAPES_GFS2.3.1 operation on the forecast results, this paper mainly analyzes the differences of simulation results between G21 and G00, and evaluates the possible changes brought to the GRAPES_GFS2.3.1. But an in-depth analysis of how the simulation results can improve the forecast performance is absent in this paper. The reason is that the forecast biases of the numerical model result from a combination of various factors, and it is difficult to explain the improvement of the GRAPES_GFS2.3.1 forecast system just with a single case. Therefore, a batch of experiments are needed later in our future study. Since the absence of cloud-field information at a single time can bring systematic biases to the simulated temperature field and geopotential height field, in the cycling numerical forecasting operational system, the cloud-field information that has formed should be retained as much as possible. Moreover, the temperature and humidity structure in the initial field, especially the water vapor, can significantly affect the dynamic and physical processes in the numerical model. Thus, in addition to the improvement of dynamic and physical processes, more attention should be paid to the assimilation of water vapor data, to improve the data quality of water vapor in the initial field of GRAPES_GFS2.3.1.