An Online Model Correction Method Based on an Inverse Problem: Part I——Model Error Estimation by Iteration

With the growing abilities of atmospheric models and the availability of computational resources, the capability of weather forecasting by numerical weather prediction (NWP) has improved dramatically in recent decades. However, the performance of NWP is generally hindered by its systematic errors, which are contributed by different sources. Besides the bias in the initial values, numerical models are defective in terms of their governing equations, numeric and physical parameterizations, and surface forcing. In recent decades, extensive effort has been made to understand and offset model biases. To handle model errors is still an enduring challenge from the model development point of view. On the other hand, from the practical point of view, researchers are making efforts to try to remove the forecast errors by model correction, regardless of possible sources of model bias.

Model correction is an efficient strategy to deal with model errors, which can be categorized as offline correction and online correction. Detailed reviews have been provided by (Danforth et al., 2007) and (Xue et al., 2013). Multiple offline corrections have been developed, such as model output statistics (MOS) (Glahn and Lowry, 1972; Carter et al., 1989), running-mean correction (Eckel and Mass, 2005; Hacker and Rife, 2007), and another correction method based on the Kalman filter (Monache et al., 2006; McCollor and Stull, 2008; Monache et al., 2011). MOS uses multi-year historical model outputs and observations to develop the offline model correction, while the others have been developed based on model output and observations for a period of less than one month. Offline bias correction does not act on the forecast dynamically, meaning the errors are allowed to interact nonlinearly throughout the model integration; whereas, online correction is conducted at the steps of the model integration to impede error growth. Moreover, it is possible to correct state-dependent errors by online correction, which has recently become a major focus in this field. (Leith, 1978) initiated a framework of an online correction method in which model bias was state-dependent. The scheme was modified and tested by (DelSole and Hou, 1999) with a nonlinear quasi-geostrophic model. The results showed that the modified Leith scheme largely improved the model forecast. Later, (Danforth et al., 2007) divided the model error into three components: model bias, periodic, and non-periodic. These components were corrected using methods including nudging, diurnal correction, and a simplified Leith scheme, respectively. The results of two models (a quasi-geostrophic model and a primitive equation model) were significantly improved.

Despite the different complexities of forecast error correction algorithms, the general idea is to estimate the forecast errors by considering the NWP as a direct problem. (Chou, 1974) proposed an alternative method by considering the NWP as an inverse problem, instead of the traditional view of the NWP. The model error (ME) due to the model deficiency in representing the real atmospheric motion is assumed as an unknown tendency term in the NWP model. If the formula of the ME can be solved inversely with observations or analyses, the model bias will be easily corrected. Unfortunately, it is impossible to know the exact mathematical law of the ME. It has been suggested that the model bias evolves linearly or quadratically in short-term forecasts (Vannitsem and Toth, 2002). Hence, the ME can be discretized into short intervals (e.g., 6 h) and considered as a constant or linear form in each interval. Given the past observations and NWP model, the discretized MEs in the past intervals can be solved iteratively as a constant or linear-increased tendency term in each interval. These MEs can be further used as the basis of the online corrections. Not only the systematic model errors, but also the diurnal model errors and state-dependent model errors, could be corrected by reasonably extrapolating the MEs in the past intervals to the forecast time. (Da, 2011) proposed the trapezoidal method to calculate the MEs in past intervals and the Lagrange polynomial to extrapolate them to the future time step. However, to apply Da's framework to a real case NWP model, there are two shortcomings. The first is that observations are supposed to be accurate, so that the order of the polynomial must be equal to the times of past data. In fact, errors are inevitably introduced by observations and analyses, which lead to Da's method becoming unpractical because a high order polynomial is sharply sensitive to the errors (Xue et al., 2013). The second is that the trapezoidal method recommended by Da makes sense in idealized cases but is not valid for real cases, because the model tendencies of the first step and last step are used. Tests have revealed that those tendencies cannot represent model errors well (not shown here).

In the present paper, we report a new method that has been developed for obtaining more accurate past MEs, which is the basis for suitable operational use of online model correction at relatively low computational expense. Part II of the series will report the development of a more steady correction formula compared to Da's polynomial. Part I, meanwhile, proceeds as follows: Section 2 describes the new iteration-approach method. Section 3 reports the results from idealized experiments and batches of real experiments to illustrate the validity of the approach and discuss the convergence of the method. Section 4 compares linearly estimated model correction based on two months of iteration results to the forecast errors of the Global and Regional Assimilation and Prediction System (GRAPES). Section 5 concludes the paper.

Following (Xue et al., 2013), the NWP model can be considered as the following initial value problem: \begin{align} \label{eq1} &\dfrac{\partial{\psi}}{\partial t}={M}({\psi}) ,(1a)\\ \label{eq2} &{\psi}|_{t=0}={\psi}_0 , (1b)\end{align} where ψ is the vector form of prognostic variables including wind, temperature, moisture and pressure; and M is the model operator. Equation (3) is the initial condition. Supposing there exists an unknown tendency term E (ME) to make the model exact, the traditional NWP model, Eq. (2), can be rewritten as an inverse problem. This unknown tendency term can be considered as the overall effect of different sources of model imperfection. To solve the inverse NWP problem, the true states of the atmosphere in the past are necessary, which could be replaced by observations or analyses. This inverse problem is written as: \begin{align*} &\dfrac{\partial{\psi}}{\partial t}={M}({\psi})+{E} ,\tag{2a}\\ &{\psi}|_{t=0}={\psi}_0 ,\tag{2b0}\\ &{\psi}|_{t=-\delta}={\psi}_{-1} ,\tag{2b1}\\ &\cdots\cdots\\ &{\psi}|_{t=-n\delta}={\psi}_{-n} .\tag{2bn} \end{align*} where ψ₀,ψ_-1,…,ψ_-n, represent the true states at the past time iδ, and δ is the time interval. We further discretize E into series of datasets and denote E_i (i=-n,-(n-1),…,-1) as the ME in the past interval between iδ and (i+1)δ. (Vannitsem and Toth, 2002) examined the short-term model error dynamics using a low-order chaotic dynamical model, and their analysis indicated that, in light of the character of the model error sources (limited to white noise or an Ornstein-Uhlenbeck process), the mean-square error follows either a linear or a quadratic evolution in short-term simulation. Therefore, the ME is assumed as a constant or linearly dependent on time in this study. The problem now is how to determine E_i for an existent NWP model. Usually, an exact E_i is barely obtainable because the NWP model is a complex nonlinear system. Here, we have designed a technique to estimate E_i through an iterative method. We define an approximate value marked as \(\overline E_i\). By iteration, it is expected that \(\overline E_i\) converges to E_i. For brevity, the following discussion focuses on one prognostic variable of the model equation at the space-discretized point without loss of generality. With \(\overline E_i\), the NWP model in the interval between iδ and (i+1)δ can be written in the following form: \begin{equation} \label{eq7} \dfrac{\partial\psi}{\partial t}={M}(\psi)+\overline{E}_i .(3) \end{equation} With the first guess of ME labeled as \(\overline E_i,0\) and the remainder marked as \(\overline E'_i,0\), so that \(\overline E_i=\overline E_i,0 +\overline E'_i,0\), Eq. (8) can be rewritten as \begin{equation} \label{eq8} \dfrac{\partial\psi}{\partial t}={M}(\psi)+\overline{E}_{i,0}+\overline{E}'_{i,0} .(4) \end{equation} Integrating Eq. (9) from iδ to (i+1)δ, we find that \begin{equation} \label{eq9} \hat{\psi}_{i+1}-\hat{\psi}_i=\int_{i\delta}^{(i+1)\delta}{{M}_1(\hat{\psi})dt}+\int_{i\delta}^{(i+1)\delta}\overline{E}'_{i,0}dt . (5)\end{equation} Here, \(\hat\psi_i\) [i=-n,-(n-1),…,-1] and \(\hat\psi\) are the truths on discretize observation time and arbitrarily time respectively, and \(M(\hat\psi)+\overline E_i,0\) is marked as \(M_1(\hat\psi)\). However, it is impossible to know \(\hat\psi\) at every time step between iδ and (i+1)δ, so \(\overline E_i,0\) cannot be calculated directly. Instead, we hope the simulation ψ will converge to the truth \(\hat\psi\), so Eq. (10) can be approximately expressed as \begin{equation} \label{eq10} \hat{\psi}_{i+1}-\hat{\psi}_i\approx\int_{i\delta}^{(i+1)\delta}{M}_1(\psi)dt+\int_{i\delta}^{(i+1)\delta}\overline{E}'_{i,0}dt . (6)\end{equation} Integrating the revised model M₁(ψ), we find that \begin{equation} \label{eq11} \overline{E}'_{i,0}\approx\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\hat{\psi}_i)-\dfrac{1}{\delta}(\psi_{i+1,1}-\hat{\psi}_i) =\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\psi_{i+1,1}) . (7)\end{equation} Here, ψ_i+1,1 is the prediction of the revised model M₁(ψ). Therefore, we make Eq. (9) closed. If the forecast of Eq. (9) is exact enough, \(\overline E_i,0+\overline E'_i,0\) is the ME in the interval. Otherwise, denoting \(\overline E_i=\overline E_i,1+\overline E'_i,1=\overline E_i,0+\overline E'_i,0+\overline E'_i,1\) as the ME after the first iteration, a similar manipulation of Eqs. (10), (11) and (12) yields a recursive formula as \begin{eqnarray} \label{eq12} E_i&=&\overline{E}_{i,0}+\overline{E}'_{i,0}+\overline{E}'_{i,1}+\cdots+\overline{E}'_{i,{\rm n}}+\cdots\nonumber\\ &=&\overline{E}_{i,0}+\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\psi_{i+1,1})+\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\psi_{i+1,2})+\cdots\nonumber\\ &&+\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\psi_{i+1,k})+\cdots,(8) \end{eqnarray} where $$ \overline{E}_{i,k}=\overline{E}_{i,k-1}+\overline{E}'_{i,k-1}=\overline{E}_{i,k-1}+\dfrac{1}{\delta}(\overline{\psi}_{i+1}-\psi_{i+1,k}) $$ is the ME after the kth iteration, $$ \overline{E}'_{i,k-1}=\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\hat{\psi}_i)-\dfrac{1}{\delta}(\psi_{i+1,k}-\hat{\psi}_i) =\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\psi_{i+1,k}) $$ is the corresponding remainder, and ψ_i+1,k is the forecast of $$ {M}_{k}(\psi)={M}_{k-1}(\psi)+\overline{E}_{i,k-1} . $$ Therefore, the iterative equation is \begin{equation} \label{eq13} \overline{E}_{i,k}=\overline{E}_{i,k-1}+\dfrac{1}{\delta}(\hat{\psi}_{i+1}-\psi_{i+1,k}) . (9)\end{equation}

It is clear that the process described by Eq. (14) is the classical first-order stationary iterative method, except that ψ_i+1,k is computed by the NWP model integration. As indicated by the recursive formula, Eq. (13), the process of iteration actually means that the forecast error (the difference between the model forecast and the truth) plays a role as a nudging term to make the forecast converge to the truth, step by step. Theoretically, as an iterative method, it should be demonstrated that the iteration can be convergent in a unique value. However, it is difficult to demonstrate this theoretically, because ψ_i+1,k is calculated by a highly nonlinear NWP system. Instead, for the application value of this method, we have determined the convergence of the iteration by idealized and real tests directly.

2.2.1. Model description

In this study, the global forecast system of the Global/ Regional Assimilation and Prediction System (GRAPES-GFS) was employed. GRAPES was designed as a unified NWP system for both research and operational applications. The main features of GRAPES can be summarized as: (1) Fully compressible governing equations with options of hydrostatic or non-hydrostatic approximation; (2) A semi-Lagrangian and semi-implicit time-stepping scheme; (8) Height-based terrain following vertical coordinate; (9) Full physical package; (10) 3DVAR/4DVAR (three-dimensional variational/four-dimensional variational) for data assimilation; (11) Modularized and parallelized code architecture; (12) Latitude-longitude grid (Chen and Shen, 2006; Xue, 2006; Zhang and Shen, 2008).

A rapid radiative transfer model is used for the shortwave and longwave radiation from Atmospheric and Environmental Research (Mlawer and Clough, 1998; Iacono et al., 2000; Iacono et al., 2008). The convective parameterization uses a one cloud type Arakawa-Schubert convection scheme (Pan and Wu, 1995; Han and Pan, 2011), whereas the microphysics scheme adopts the WRF single-moment six-class method (Rutledge and Hobbs, 1983; Hong et al., 2004). The common land model (COLM) is used for the surface processes (Dai et al., 2003; Dai et al., 2004) and an implicit approach for vertical flux divergence for the planetary boundary layer processes (Troen and Mahrt, 1986; Hong and Pan, 1996).

The prognostic variables of GRAPES-GFS include horizontal velocity (u and v), vertical velocity in the hight-based terrain following vertical coordinate system, perturbation from reference potential temperature θ’, Exner pressure perturbation \(\Pi’\), and water substances.

Figure 1.  Flow chart of the iteration method.

2.2.2. Application to estimate model error of GRAPES-GFS

As mentioned previously, iteration needs the true states of the atmosphere, which are impossible to obtain. Hence, a surrogate has to be sought for the truth. The observation is what most resembles the truth, but it first has to be interpolated into the model grids and converted to the model variables for the application of iteration. This will inevitably involve interpolative errors and bring problems for use compared with analyses, which are assimilated from observations and model simulations by model dynamics. Although analytical errors are inevitably embodied in analyses, it is most economical and practical to employ an analysis dataset. Following (Danforth et al., 2007), the initial data and verification data used in this study were NCEP (National Centers for Environmental Prediction) FNL (final) data. These NCEP FNL Operational Global Analysis data are on 1^°× 1^° grids prepared operationally every six hours. This product is from the Global Data Assimilation System, which continuously collects observational data from the Global Telecommunications System, and other sources, for many analyses. The analyses are available at the surface and at 26 levels from 1000 hPa to 10 hPa. Variables include surface pressure, sea level pressure, geopotential height, temperature, sea surface temperature, soil state variables, ice cover, relative humidity, horizontal winds, vertical motion, vorticity and ozone. The NCEP RTG_SST (real-time global surface sea temperature) analysis dataset is used to drive GRAPES-GFS with a 1^°× 1^° horizontal resolution.

The model initialization module was employed to interpolate and convert the analyses into the variables on the GRAPES-GFS grids. Due to the credibility problem of vertical velocity in observations or analyses, the iteration will take no account of the vertical velocity. Additionally, there are several water substances in GRAPES-GFS (such as cloud water, rain water and ice crystal) and only one prognostic variable (relative humidity) in FNL analyses. Therefore, the iteration described above applied in GRAPES-GFS involved forecast variables of u, v, θ’ and \(\Pi’\). Given the model forecast and the verification, the revised ME could be calculated by running a specific module. The model resolution was set to 1^°× 1^° and the model time step to 10 minutes. Figure 1 is a flow chart of the iteration process, which can be summarized as involving the following steps:

(1) Integrate the model and determine the exit condition. The model forecast error (difference between model forecast and analysis) would be a good index to judge the condition and, sequentially, a critical value is necessary. Instead, the iterative step can be another index to determine the exit condition. For simplicity, 20 iterative steps were set as the exit condition in this study.

(2) Obtain the first guess of ME \(\overline E_i\). If the exit condition is not satisfied, the ME should be calculated. As the first step of the iteration, the first guess of ME needs to be provided. One approach is the so-called trapezoidal method (Xue et al., 2013), which uses the equation $$ \overline{E}_{i,0}\!=\!(\hat{\psi}_{i+1}\!-\!\hat{\psi}_i)/\delta t\!-\![(\psi_{i,0}\!-\!\hat{\psi}_i)/dt+(\psi_{i+1,0}-\hat{\psi}_{i+1})/dt]/2 $$ to calculate \(\overline E_i,0\). Here, \((\psi_i,0-\hat\psi_i)/dt\) and \((\psi_i+1,0-\hat\psi_i+1)/dt\) are the model tendencies of the first model step at the corresponding time iδ and (i+1)δ. The more direct approach considers the time weighted forecast error at time (i+1)δ as \(\overline E_i,0\), i.e., $$ \overline{E}_{i,0}=(\psi_{i+1,0}-\hat{\psi}_{i+1})/\delta t . $$

(3) Apply the ME in the model. \(\overline E_i,0\) is actually a database including the model error of u, v, θ’ and \(\Pi’\) at every model grid, which need to be input before model integration. Then, for an arbitrary step, \(\overline E_i,0\) is taken as a tendency term and multiplied by the model step interval to update the corresponding prognostic variables after a full integration step (i.e., dynamics, physics and filter etc.).

(4) Continued iteration. With the new model output at (i+1)δ, if the exit condition is not satisfied, the new ME can be calculated based on Eq. (14) and then repeat step (8).

(5) Exit. If the exit condition is satisfied, the iterative loop is exited.

Before application, the convergence of the iteration should first be confirmed. To fulfill this objective requires the discussion of two uncertainties. The first is confirming the approaches to calculate the first guess of the ME in step (1), and the second is how to choose the form of the ME in step (2).

For the second issue, the form of the ME should be based on the real model forecast error evolution in a short interval (such as δ=6 h). As mentioned earlier, the mean-square error initially follows either a linear or a quadratic evolution in the short-term dynamics of ME using a low-order chaotic dynamical model (Vannitsem and Toth, 2002). Therefore, the ME was assumed as constant or a linear form here on. For the constant case, the value of the first guess could be calculated by the direct method and trapezoidal method discussed previously, and the updated values were based on Eq. (14). For the linear case, the form was taken as A× t, where t is time and A is a coefficient. The integrated effect of the linear term should be equivalent to that of the constant term. Therefore, the coefficient A could be obtained through $$ \int_{i\delta}^{(i+1)\delta}A_ktdt=\int_{i\delta}^{(i+1)\delta}\overline{E}_{i,k}dt , $$ where k is an arbitrary iterative step and \(\overline E_i,k\) can be calculated in the same way as in the above constant case. For the purpose of convergence discussion, it was necessary to perform some idealized experiments, the details of which are provided in Table 1. These tests were all idealized in order to obtain the proper first guess of the iteration and form of the ME applied in the model. All the tests were initiated using NCEP 1^°× 1^° FNL data at UTC 0000 6 June 2009, and the iteration intervals were 6 h.

Test1 was designed to determine the validity of the trapezoidal method, in which the first guess of the ME was produced by the trapezoidal method. In addition, the model forecast was considered as the truth, meaning the ME should converge to zero if the trapezoidal method was valid. In Test2 and Test3, artificial model error sources were introduced in GRAPES-GFS, and both constant and linear forms were tested. In Test2, the constant form of the ME was employed, while the linear form was used in Test3. These two tests were designed in order to choose a proper form of the ME, countering either constant or linear error source. In Test2 and Test3, the first guesses of the MEs all employed the direct method described in the previous section.

The forecast errors were calculated based on the difference between the model forecast and the truth (or analysis) to examine the performance of the iteration. Figure 2 shows the mean absolute forecast errors of u, v, θ’ and 1000 \(\times\Pi’\) at hour 6 with iteration steps for the designed tests. Note that the mean absolute forecast errors of the four variables increased with the iterative steps, as shown in Fig. 2a. Actually, some other tests have also shown the same results for the trapezoidal method, whatever form of the ME is chosen (not shown). These results indicate that the iteration is divergent if the trapezoidal method is employed to obtain the first guess. Similarly, applying the ME as a constant term into the model also made the iteration divergent, both for the constant form of error source (shown in Fig. 2b) or the linear form (not shown). Conversely, applying the ME as a linear form resulted in the iteration being convergent against the constant form of error source (not shown) or the linear form (Fig. 2c). The reason might be that the model needs to gradually adapt the ME from zero to the actual magnitude.

Figure 2.  Mean absolute forecast errors of the field variables u (m s^-2), v (m s^-2), θ’ (^°C) and 1000 \(\times\Pi’\) at 6 hour: (a) Test1; (b) Test2; (c) Test3.

Figure 3.  Mean absolute forecast errors of the field variables u (m s^-2, black), v (m s^-2, black), θ’ (^°C, green) and 10000 \(\times \Pi’\) (blue) at hour 6 involved in the iterative steps for (a) JA 2009 and (b) JF 2010.

Figure 4.  Evolution of mean absolute forecast errors of u (m s^-2, red), v (m s^-2, blue), T (^°C, green) and Z (gpm, black) at hour 6 at the 500 hPa pressure surface before iteration (upper solid lines) and after iteration (lower dotted lines) for (a) JA 2009 and (b) JF 2010. The abscissa shows the date.

Based on the theory of limit, if only one maximum (or minimum) point exists, the numerical solution will be convergent to the solution no matter what the first guess is; while there are several maximum (or minimum) points, the first guess should influence the results. For choosing the first guess in this study, it should represent the model error as accurately as possible. Theoretically, the trapezoidal method recommended by (Da, 2011) is more reasonable, but it is based on the two model tendencies at the beginning and end of the model integration according to the formula, which may introduce too many errors. We have examined the first guess using the trapezoidal method, and the value of it was found to be unreasonable (not shown). Another technical problem is how to apply the ME in the model. We assumed the form of the ME as a constant in a specific interval, as described in the previous section, but it caused the iteration to be divergent. Applying the ME as a constant, the NWP model was impacted at the very beginning of the model integration, which was unable to make the complicated and nonlinear model adjust to the outside influence immediately. Meanwhile, using the linear form of the ME to obtain the first guess of the ME, the forecast errors were linearly distributed in the integrated interval, which caused the model to be gradually influenced by the ME. Therefore, the iteration became converged. In summary, using the trapezoidal approximation for the first guess or applying the ME as a constant term makes the iteration divergent, while using the linear form of the ME together with the direct approach to obtain the first guess of the ME allows the iteration to become converged.

In this part, results are reported from carrying out iteration experiments with the implementation of the ME tendency term into GRAPES-GFS for July-August (JA) 2009 and January-February (JF) 2010. The iteration of each case was carried out at every six hours. Again, the initial data and verification data were from the NCEP 1^°× 1^° FNL dataset. The ME was set in its linear form and the first guess used in the direct way. Instead of setting a critical iterative value, all cases were simply iterated for 20 steps. That meant there were four tests per day, and each test iterated 20 times. Every iteration step needed to integrate the model for 6 h. To verify the performance of the iteration, the forecast errors at hour 6 were used, which were the differences between the forecasts at hour 6 and NCEP analyses.

Figure 5.  Mean wind differences (u and v, m s^-2) between model forecasts of GRAPES-GFS at hour 6 and NCEP FNL before iteration (a, c) and after iteration (b, d) at the 500 hPa pressure surface for JA 2009. The errors after iteration have been amplified by 10.

Figure 6.  Zonally averaged latitude-height cross sections of errors of u (m s^-2), v (m s^-2) and T (^°C) between model forecasts of GRAPES-GFS at hour 6 before iteration (a, c, e) and after iteration (b, d, f) for JA 2009. The errors after iteration have been amplified by 10. The vertical axis shows pressure.

Figure 7.  Mean forecast errors and estimated mean error corrections of 10000 \(\times\Pi’\) at hour 6 at the 13th model level for JF 2010. Mean forecast errors are the mean differences between model forecasts at hour 6 and analyses, while the estimated mean error corrections are calculated by linearly multiplying time on the mean MEs.

Figure 8.  Zonally averaged latitude-height cross sections of forecasts errors (left column) and estimated error corrections (right column) of u (m s^-2), v (m s^-2), θ’ (^°C) and 2000 \(\times\Pi’\) at hour 6 for JF 2010. The vertical axis shows the model level.

Figure 9.  Forecast errors and estimated error correction evolution of (a, b) 10000 \(\times\Pi’\) and (c, d) θ’ (^°C) for JF 2010.

The global mean absolute forecast errors of the field variables u, v, θ’ and \(10000\times \Pi’\) at hour 6 for JA 2009 and JF 2010 are shown in Fig. 3. The bunches of lines represent horizontal velocity, perturbation from reference potential temperature and Exner pressure perturbation. Each line, which is related to a single test's single variable, tends towards zero as the iteration step increased. Clearly, the iterations of all of the cases were convergent for all the variables in the experiments of JA 2009 and JF 2010. The evolution of the global mean absolute forecast errors of u, v, T and Z at hour 6 at 500 hPa before and after iteration are also shown (Fig. 4). The errors reduced to about 10% of their original value after 20 steps of iteration for both JA 2009 and JF 2010. It is interesting that the forecast errors before iteration remained flat with an obvious diurnal cycle evolution for both JA 2009 and JF 2010. This implies that, besides climatological systematic errors, the component of the MEs related to the diurnal cycle should not be ignored. Besides systematic errors, the diurnal errors were reduced substantially after iteration. The amplitude decreased greatly compared to the original mean error.

To support the conclusion of iterative convergence, Fig. 5 shows the mean difference of wind at 500 hPa between the forecasts of GRAPES-GFS at hour 6 and NCEP FNL for JA 2009. Note that the magnitude of the errors after iteration has been amplified 10 times, but the errors are the same size as those of the original errors. The results of other variables were similar (not shown). This means that the errors of each variable at the isobaric surface were convergent by iteration. The zonally averaged latitude-height cross sections of errors of u, v and T of GRAPES-GFS for JA 2009 are shown in Fig. 6. Similarly, the magnitude of the errors after iteration has again been amplified 10 times. Apparently, the zonally averaged errors were all convergent by iteration for vertical sections. The results for JF 2010 were also similar to the result for JA 2009 (not shown).

Based on the results reported in section 3, the convergence of iteration relating to all variables of the cases was confirmed for the isobaric surface and vertical sections, and the forecast errors at hour 6 were reduced to 10% of their original value after 20 steps of iteration. To understand whether the MEs can represent the model forecast errors, it was necessary to compare the mean MEs to the mean forecast errors. As we know that the MEs are a kind of model tendency, the mean MEs therefore need to multiply time to estimate the off-line model forecast error corrections. As pointed out by (Jung and Tompkins, 2003) and (Jung, 2005), systematic errors approximately evolve linearly in short-term forecasts; therefore, we estimated the error correction by linearly multiplying time to the 2-month mean value of MEs. In this section, the results from performing offline comparisons between the mean error corrections estimated by MEs and mean forecast errors are reported.

Figure 7 shows the mean forecast errors and estimated error corrections of 13-model-level \(10000\times \Pi’\) for JF 2010. Note that GRAPES-GFS underestimates the pressure over high-latitude regions and the Arctic, but pressure over low-latitude regions is over-estimated. That means the pressure gradient along the equator to pole was underestimated. Moreover, with this pressure error pattern, the errors tended to increase linearly. It should be noted that the patterns of the mean forecast errors and estimated error corrections were similar but the signs were opposite. For example, the positive forecast errors corresponded to negative estimated correction in the Arctic region, and the negative forecast errors were related with the positive estimated correction in the midlatitudes. Figure 8 presents the zonally averaged latitude-height cross sections of forecast errors and estimated error corrections of u, v, θ’ and \(10000\times \Pi’\) at hour 6 for JF 2010. The patterns of the forecast errors and estimated errors matched well but the signs were also opposite, except for the errors of v and \(\Pi’\) at the boundaries of some regions. These results indicate the forecast errors can be expected to be offset by introducing online error correction into the model integration. However, the estimated error corrections increased quicker than the forecast errors (Figs. 7 and 9). Figure 9 shows the forecast and estimated root-mean-square error evolution of \(\Pi’\) and θ’ at the 13th model level for JF 2010. The root-mean-square errors evolved linearly for \(\Pi’\). It can be seen that the estimated root-mean-square errors of \(\Pi’\) increased quicker than the systematic errors by about fivefold in both hemispheres, and about tenfold in the tropics. As for θ’, the error evolution rates were linear during the first two days, and then decreased slightly after the two-day forecast. Similar to \(\Pi’\), the estimated error rates of θ’ were also steeper than the forecast error's. These results suggest that the patterns of the estimated model error corrections can match well to the mean forecast errors but the linear estimation should be weighted in real model correction.

To improve NWP model performance, an unknown tendency term (ME) was introduced into GRAPES-GFS, which was presumed to represent the imperfection of the NWP model. This part (Part I of this series) of the study presents an iterative method for obtaining MEs in the past to provide the basis for further uses. These MEs were compared with model forecast errors and were expected to offset the systematic MEs. The following key conclusions can be drawn:

(1) Having developed this new iteration method to obtain the MEs in past intervals, and having tested its convergence in idealized experiments, the indication was that using the linear form of the ME together with the direct method for the first guess caused the iteration to be convergent.

(2) Batches of iteration test results indicated that the forecast errors at hour 6 were reduced to 10% of their original value after 20 steps of iteration.

(8) By comparing the error corrections estimated by MEs to the mean forecast errors, the patterns of estimated errors were considered to agree well with those of the forecasts, but the linear growth rates of the estimated errors were sharper than those of the forecasts. Therefore, it can be expected that a larger proportion of the forecast errors will be offset by properly introducing the estimated ME correction into GRAPES-GFS.

The purpose of this study was to provide a method and basis for estimating the MEs in past intervals, to be used to correct model forecasts. Since forecast errors are composed of different components of different frequencies, they should be decomposed into different components and dealt with separately. In this study, the systematic errors were the focus, which could be canceled out by introducing the corresponding error correction into the model for another study. However, it should also be noted that the forecast errors at hour 6 evolved flatly with clear diurnal fluctuation (Fig. 4). This implies that, except for long existent climatological systematic errors, the errors relating to the diurnal cycle and initial state cannot be ignored. Therefore, spectral analysis is required for the forecast errors and different strategies will be introduced associated with different components of errors.