Forecasting experiments of a dynamical-statistical model 1 of the sea surface temperature anomaly field based on the 2 improved self-memorization principle 3

Forecasting experiments of a dynamical-statistical model 1 of the sea surface temperature anomaly field based on the 2 improved self-memorization principle 3 Mei Hong, Xi Chen, Ren Zhang, Dong Wang, Shuanghe Shen and Vijay 4 P. Singh 5 1 Institute of Meteorology and Oceanography, National University of Defense Technology, Nanjing 211101, China 6 2 Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disaster, Nanjing University of 7 Information Science &Technology, Nanjing 210044, China 8


Introduction
The El Niño-Southern Oscillation (ENSO), the well-known coupled atmosphere-ocean phenomenon, was firstly proposed by Bjerknes (1969).The ENSO phenomenon can influence regional and global climates, so the prediction of ENSO has received considerable public interest (Rasmusson and Carpenter, 1982;Glantz et al., 1991).
Over the past two to three decades, one might reasonably expect the ability to predict warm and cold episodes of ENSO at short and intermediate lead times to have gradual improvement (Barnston et al., 2012).Many countries have been focusing on ENSO forecasts since the 1990s, and the ENSO forecast has become one of the important research topics in the International Climate Change and Predictability Research plan.The US International Research Institute for Climate and Society, US Climate Prediction Center, Japan Meteorological Agency (JMA) and European Centre for Medium-Range Weather Forecasts (ECMWF) have developed different coupled atmosphere-ocean models to Published by Copernicus Publications on behalf of the European Geosciences Union.
The forecast models can generally be divided into two types (Palmer et al., 2004;Zheng et al., 2017).The first type is typified by a dynamical model, which mathematically expresses physical laws that govern how the ocean and the atmosphere interact.The second type is typified by a statistical model, which requires a large amount of historical data and analyses the data to forecast (Chen et al., 1995;Moore et al., 2006).
Over the past three decades, ENSO predictions have made remarkable progress, reaching a stage where reasonable statistical and numerical forecasts (Jin et al., 2008) can be made 6-12 months in advance (B.Wang et al., 2009).However, there are three problems remaining to be resolved (Zhang et al., 2003a).(1) The current ENSO predictions are mainly limited to the short term, such as annual and seasonal predictions.(2) Although the representation of ENSO in coupled models has been advanced considerably during the last decade, several aspects of the simulated climatology and ENSO are not well reproduced by the current generation of coupled models.The systematic errors in sea surface temperature (SST) are often very large in the equatorial Pacific, and model representations of ENSO variability are often weak and/or incorrectly located (Neelinet al., 1992;Mechoso et al., 1995;Delecluse et al., 1998;Davey et al., 2002).
(3) Coupled models of ENSO predictions initialized from observed initial states tend to adjust towards their own climatological mean and variability, leading to forecast errors.The errors associated with such adjustments tend to be more pronounced during boreal spring, which is often called the "spring predictability barrier" (Webster et al., 1999).More efficient models are therefore desired (Belkin and Niyogi, 2003;Weinberger and Saul, 2006).Therefore, the idea of combining dynamical and statistical methods to improve weather and climate prediction has been developed in many studies (Huang et al., 1993;Yu et al., 2014a, b).By introducing genetic algorithms (GAs), Zhang et al. (2006) inverted and reconstructed a new dynamical-statistical forecast model of the tropical Pacific SST field using historic statistical data (Zhang et al., 2008).However, there is one flaw in the forecast model: the time-delayed SST field.It is because that ENSO is a complicated system with many influencing factors.To overcome information insufficiency in the forecast model, Hong et al. (2014) selected the tropical Pacific SST, sea surface wind (SSW) and sea level pressure (SLP) fields as three modeling factors and utilized the GA to optimize model parameters.
However, the above dynamical prediction equations, which were proposed by Hong et al. (2014), greatly depend on a single initial value, creating long-term forecasts over 8 months that diverged significantly.These unsatisfactory results indicate that this model needs to be improved.Cao (1993) first proposed the self-memorization principle, which transforms the dynamical equations with the self-memorization equations, wherein the observation data can determine the memory coefficients.This method has been widely used in forecast problems in environmental, hydrological and meteorological fields (Feng et al., 2001;Gu, 1998;Chen et al., 2009).The method can avoid the question of initial conditions for the differential equations, so it can be introduced here to improve the proposed dynamical forecast model.
Therefore, an improved dynamical-statistical forecast model of the SST field and its impact factors with a selfmemorization function was developed.The improved model can absorb the information from past observations.This paper is organized as follows: research data and forecast factors are introduced in Sect. 2. In Sect.3, the reconstruction of the dynamical model of the sea surface temperature anomaly (SSTA) field is described.To improve the reconstruction model, the self-memorization principle is introduced in Sect. 4. Model forecast experiments are described in Sect.5, and conclusions are given in Sect.6.

Data
The monthly average SST data were obtained from the UK Met Office Hadley Centre for the region of 30 • S-30 • N, 120 • E-90 • W. The gridded 1 • ×1 • Met Office Hadley Centre sea ice and SST data set (HadISST1; Rayner et al., 2003) includes both in situ and available satellite data.The sea areas provide important information on ocean-atmosphere coupling in the east and west Pacific Ocean and the El Niño/La Niña events.The reanalysis data, zonal winds and sea level pressures were obtained from the National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (Kalnay et al., 1996).The sea surface height (SSH) field was obtained from Simple Ocean Data Assimilation (SODA) data (Carton and Giese, 2008).Outgoing longwave radiation (OLR) was obtained from the National Oceanic and Atmospheric Administration (NOAA) satellites, at a resolution of 0.5 • ×0.5 • (Liebmann and Smith, 1996).The Southern Oscillation Index (SOI) data were obtained from the Climate Prediction Center (CPC).The time series of all data were from January 1951 to December 2010 (720 months in total).

EOF deconstruction
The SSTA field can be calculated from the SST field and can be deconstructed into time (coefficients)-space (structure) using the empirical orthogonal function (EOF) method.Detailed information on the EOF method can be seen in the related references (Dommenget and Latif, 2002).We have used covariance matrix, because the covariance matrix was selected to diagnose the primary patterns of covariability in Ocean Sci., 14, 301-320, 2018 www.ocean-sci.net/14/301/2018/ the basin-wide SSTs, rather than the patterns of normalized covariance (or correlation matrix).We used the smoothing function with MATLAB to smooth the SSTA field before the EOF deconstruction, which is moving from five points, mainly filtering out some noise points and outliers.Then, an EOF analysis of smoothed anomalies was performed, and the first two SSTA EOFs are shown in Fig. 1a and c.The principal component (PC) time series corresponding to the first and second EOFs are shown in Fig. 1b  and d.The first EOF pattern, which accounted for 61.33 % of the total SSTA variance, represented the mature ENSO phase (El Niño or La Niña), and the corresponding PC time series was highly correlated (with a correlation coefficient of 0.85) with the cold tongue index (SST anomaly averaged over 4 • S-4 • N, 180-90 • W) over the whole period.The second EOF, accounting for 14.52 % of the total SSTA variance, indicated the ENSO signal beginning to enhance.Compared with the first mode, these were slightly attenuated in terms of the scope and intensity.The above analysis is similar to the EOF analysis of the SSTA field in the previous studies (Johnson et al., 2000;Timmermann et al., 2001).This indicates that the front two variance contribution modes can describe the main characteristics of the SSTA field and El Niño/La Niña.Therefore, we can choose the T 1 , T 2 time series EOF decomposition modes as the modeling objects.

Selection of other prediction model factors
Considering the complexity of computation, the amount of variables in the equations of our model cannot be too large (usually three or four variable are best).This has been explained in our previous studies (Zhang et al., 2006(Zhang et al., , 2008)).If there are more than four variables in the modeling equation, it will cause the amount of parameters, such as a 1 , a 2 , . .., a n , b 1 , b 2 , . .., b n , . .., to be too large.The huge computation makes it difficult to be precisely modeled.Thus, the total number of parameters in the model of five variables was 102, which may cause an overfitting problem.Hence, when we selected the model of five or six variables, it entailed large amounts of computation that made precision difficult, and too many parameters might cause an overfitting phenomenon.If we choose only two or even fewer variables, the forecast performance is poor too.Too few variables cause reconstructed parameters to be too small, resulting in amounts of important information missing out in the model.Thus, four variables are best for modeling dynamically and accurately.Because we have chosen two time series in Sect.2.2 as the modeling objects, now we should select the other two ENSO intensity impact factors.
The ENSO intensity impact factor is an important issue in ENSO prediction.Previous studies have been completed in this area, which found that teleconnection patterns, temperature, precipitation, wind and SSH may affect ENSO strength.For example, Trenberth et al. (1998) noted that Pacific North American teleconnection (PNA), SOI and OLR in the Pa-cific Intertropical Convergence Zone (ITCZ) are all closely related to ENSO.Webster (1999) pointed out that, after 1970, the Indian Ocean dipole (IOD) was not only affected by ENSO but also affected the strength of ENSO (Ashok et al., 2001).Yoon and Yeh (2010) reported that the Pacific Decadal Oscillation (PDO) disrupts the linkage between El Niño and the following Northeast Asian summer monsoon (NEASM) by inducing the Eurasian pattern in the mid-to high latitudes.The vast majority of studies (Tomita and Yasunari, 1996;Zhou and Wu, 2010;Kim et al., 2017) have concentrated on the impacts of ENSO on the East Asian winter monsoon (EAWM).During the EAWM season, ENSO generally reaches its mature phase and has the most prominent impact on the climate.B. and C. Wang et al. (1999) suggested that the zonal wind factors in the eastern and western equatorial Pacific play a critical role in the phase of transition of the ENSO cycle, which could excite eastward propagating Kelvin waves and affect the SSTA in the equatorial Pacific.Zhao et al. (2012) analyzed the characteristics of the tropical Pacific SSH field and its impact on ENSO events.
Based on the above analysis, we have selected nine factors, which may be closely related with the ENSO index (Niño3.4).
1.The zonal wind in the eastern equatorial Pacific factor (u1) was calculated as the grid-point average of zonal wind in the area of 5 • S-5 • N, 150-90 • W.
2. The zonal wind in the western equatorial Pacific factor (u2) was calculated as the grid-point average of zonal wind in the area of 0-0 • N, 135-180 • E.
3. The PNA teleconnection factor was obtained from the CPC.
4. The dipole mode index factor (DMI) was obtained from SSTA for June-July-August (JJA) based on the Saji (1999) method.
5. The SOI factor was obtained from the CPC.
6.The PDOI factor was obtained from the Department of Atmospheric Sciences at the University of Washington.The website is http://tao.atmos.washinton.edu/pdo/RDO.latest.
8. The OLR in the ITCZ factor was calculated as the gridpoint average of OLR in the area of 10-20 • N, 120-150 • E. 9.The SSH factor was calculated as the grid-point average of the SSH data in the area of 10   A correlation analysis of the above factors was carried out and the results are shown in Table 1.
Table 1 shows that SOI and EAWMI have the stronger correlation with the front two time series (T 1 , T 2 ) than the other seven factors.The results are also consistent with previous research (Clarke and Van Gorder, 2003;Drosdowsky, 2006;Zhang et al., 1996;Wang et al., 2008;Yang and Lu, 2014).Therefore, the first time series (T 1 ), the second time series (T 2 ), SOI and EAWMI will be selected as prediction model factors.
3 Reconstruction of dynamical model based on GA Takens' delay embedding theorem (Takens, 1981) provides the conditions under which a smooth attractor can be constructed from observations made with a generic function.Later results replaced the smooth attractor with a set of arbitrary box-counting dimensions and the class of generic functions with other classes of functions.Takens had shown that if we measured any single variable with sufficient accuracy for a long period of time, it would be possible to construct the underlying dynamical structure of the entire system from the behavior of that single variable using delay coordinates and the embedding procedure.It was therefore possible to construct a dynamical model of system evolution from the observed time series.Introducing this idea here, four time series of the T 1 , T 2 , SOI and EAWMI factors were chosen to construct the dynamical model.
The basic idea of statistical-dynamical model construction is discussed in Appendix A and was introduced in our previous work (Zhang et al., 2006;Hong et al., 2014).
A simplified second-order nonlinear dynamical model can be used to depict the basic characteristics of atmosphere and ocean interactions (Fraedrich, 1987).Suppose that the following nonlinear second-order ordinary differential equations are taken as the dynamical model of reconstruction.In the equations, x 1 , x 2 , x 3 and x 4 were used to represent the time coefficient series of T 1 , T 2 , SOI and EAWMI.
Based on the parameter optimization search method of GA in Appendix A, the time coefficient series of T 1 , T 2 , SOI and EAWMI from January 1951 to April 2008 are chosen as the expected data to optimize and retrieve model parameters.In order to eliminate the dimensionless relationship between variables, data standardization x nor = (x − x min )/(x max −x min ) involves transforming data from different orders of magnitude to the same order of magnitude.Finally, we made forecast results revert back to the raw data magnitude by x = x nor (x max − x min ) + x min .
In order to quantitatively compare the relative contribution of each item of our model to the evolution of the system, we calculated the relative variance contribution.The formula is as follows: where n is the length of the data, and T i = a 1 x 1 , a 2 x 2 , . .., a 14 x 3 x 4 is the item in the equation.According to our previous research (Hong et al., 2007), the variance contribution of the real item reflecting the performance of the model has a large proportion, while the variance contribution of the false term is almost zero, so we delete the weak items of R i < 0.01.
After deleting the weak items, the nonlinear dynamical model of the first time series (T 1 ), the second time series (T 2 ), SOI and EAWMI can be reconstructed as follows: (2) The model required testing.Because the training period was from January 1951 to April 2008, we chose T 1 , T 2 , SOI and EAWMI of May 2008, which were not used as initial forecast data in the modeling.Next, the Runge-Kutta method was used to do the numerical integration of the above equations, and every step of the integration was regarded as 1 month's worth of forecasting results.As a result, forecast results of four time series over a period of 20 months were obtained.
Here, the focus was on the forecast results of T 1 and T 2 , as shown in Fig. 2. The Pearson correlation coefficient (CC) (W.C. Wang et al., 2009) and the mean absolute percentage error (MAPE) (Hu et al., 2001) are employed as objective functions to calibrate the model.The CC evaluates the linear relationship between the observed and predicting values and MAPE measures the difference between the observed and predicting values.
From Fig. 2, forecast performance of T 1 and T 2 within 5 months was better.Using T 1 as an example, the CC between model predictions and corresponding observations over the first 5 months of forecasts was 0.8966 and MAPE was 8.32 %.However, after 5 months, MAPE increased rapidly and was 31.29 % at 10 months.The model forecast then significantly diverged from observations, and the forecast became inaccurate.After 10 months, the forecast results became increasingly worse, which indicated that the forecast of the model after 5 months was unacceptable.The forecast results of T 2 were similar to those of T 1 .
The model's skill should be further assessed by crossvalidated retroactive hindcasts of the time series.As in the above example, omitting a portion of the time series (12 months, January 1951 to December 1951) from observations, we trained the model based on the data from January 1951 to December 2010 and then predicted the omitted segments (12 months, January 1951 to December 1951).Then, in the next prediction experiment, the omitted segment was from January 1952 to December 1952 and the training samples were from January 1951 to December 1951 and January 1953 to December 2010.So the forecast time series is from January 1952 to December 1952.We then repeated this procedure by moving the omitted segment along the entirety of the available time series.Each experiment has used a different training sample and established a different model equation (but the method is the same).The similar process of the cross-validated retroactive hindcasts has also been used in the previous literatures (Hu et al., 2017).
Finally, we obtained cross-validated retroactive hindcast results of T 1 and T 2 , as shown in Fig. 3.So the forecast results of 60 cross experiments (each experiment is the prediction of the 12 months as in Fig. 2) according to the time sequence can merge into a new time series (from January 1951 to December 2010), and then CC and MAPE can be calculated by the new prediction time series and the time series of the actual value.Figure 3 shows combined results of the 60 forecast experiments.As in Fig. 2, the forecast performance of T 1 and T 2 in Fig. 3 was not satisfactory.The model forecast significantly diverged from observations, and the forecast became inaccurate.The CCs of T 1 and T 2 between model predictions and corresponding observations were 0.3411 and 0.4176, respectively.Additionally, the MAPEs of T 1 and T 2 were 65.42 and 57.56 %, respectively.This indicates that the forecast of the model in the long term was inaccurate and unacceptable.
The forecast result may be inaccurate when the integral forecasting time is long.There will be a significant divergence which will cause an ineffective forecast.To improve the forecast accuracy, the forecast not only depends on the integral equation but also on a single initial value.Choosing the different initial value will cause different forecast accuracy.For example, in a total of 60 cross-validated retroactive hindcast examples, the minimum MAPE was 37.65 %, while the maximum MAPE was 89.88 %.A forecast, depending on a single initial value, will cause instability of the forecast re-sults.These two problems are addressed by introducing the self-memorization principle in the next section.

Introduction of self-memorization dynamics to improve the reconstructed model
In the above discussion, it was shown that the accuracy of the forecast results of Eq. ( 2) were unsatisfactory.To improve long-term forecasting results, the principle of selfmemorization can be introduced into the mature model (Gu, 1998;Chen et al., 2009).The principle of self-memorization dynamics (Cao, 1993;Feng et al., 2001) can be seen in Appendix B.
Based on Eq. (B10) in Appendix B, the improved model can be expressed as follows: where y i is replaced by the mean of two values at adjoining times; i.e., y i ≡ 1 2 (x i+1 + x i ).F is the dynamical core of the self-memorization equation, which can be obtained from Eq. (2); and α and θ are the memory coefficients, the formula for which can be found in Appendix B.
If the values of α and θ can be obtained, Eq. ( 3) can be used to obtain the results of final prediction.The memory coefficients α and θ in Eq. (3) were calibrated using the least-squares method with the same data (January 1951 to April 2008) as those used in Sect.3. Equation (3) can be deconstructed as follows (M is the length of the time series): where Equation ( 4) can be written as (5) The memory coefficient vector W can be calibrated using the least-squares method: The memory coefficients a, θ can be obtained from Eq. ( 6).
We then made a prediction using the self-memorization Eq. ( 3), which used the p values before t 0 .
The coefficients in F and W were used with the same training data from January 1951 to April 2008.In the forecast examples, we trained both of the coefficients in F and W at the same time, but in the paper we describe them separately to facilitate better understanding for the reader.The training sample for the model was from January 1951 to April 2008.Here, from Eq. (3), the forecast results using T 1 , T 2 , SOI and EAWMI factors can be calculated; this is called a step-by-step forecast.
When the retrospective order p is confirmed, step-by-step forecasts can be carried out.For example, when the T 1 , T 2 , SOI and EAWMI values of May 2008 were forecast, y i was obtained from the previous p + 1 time of T 1 , T 2 , the SOI and the EAWMI data, and F i (x 1i , x 2i , x 3i , x 4i ) was obtained from the previous p times of T 1 , T 2 , the SOI and the EAWMI data.All four equations were integrated simultaneously.Taking these in Eq. ( 3), we can get the T 1 , T 2 , SOI and EAWMI values of May 2008, which these can be taken as the initial values for the next prediction step.Then, the T 1 , T 2 , SOI and EAWMI values from June 2008 and so on can be generated.

Determination of p
Based on the self-memorization principle, the selfmemorization of the system determines the retrospective order p (Cao, 1993).If the system forgets slowly, parameters a and θ will be small and the p value should be high.The SSTA field forecasts were on a monthly scale, the change of www.ocean-sci.net/14/301/2018/Ocean Sci., 14, 301-320, 2018 which was slow in contrast to large-scale atmospheric motion.So parameters a and θ were small, and generally, the p value was in the range of 5 to 15 (Cao, 1993).The retrospective order p was obtained by a trial calculation method.We selected the p values in the range of 4 to 16 to construct the model.The CCs and MAPEs of the longterm fitting test (from February 1951 to December 2010) are shown in Table 2, which can be used as the standard to determine the retrospective order p.
Table 2 indicates that when p = 6, the MAPE values of the long-term fitting test were the smallest and the CCs were the largest.Also, when p ranged from 5 to 9, the CCs were all more than 0.58, and the forecast results were all good, which is consistent with our interpretation of the physical mechanisms in Sect.6.2 below.SOI and the East Asian winter monsoon index (EMWMI) had 5-to 12-month lead relationships with SST (Xu et al., 1993;Chen et al., 2010;Wang et al., 2003).Using a cumulative period of SOI, EMWMI was 5-8 months ahead, as initial values can help improve the final forecast results.Our results in Table 2 are consistent with the actual physical ENSO process.Therefore, we selected the retrospective order as p = 6.
Then, the prediction experiments can be carried out, based on improved self-memorization (Eq.3).
The improved self-memorization equation of T 1 , T 2 , SOI and EAWMI can then be established.After the differential equation was discretely dealt with, the memory coefficients were solved by the least-squares method given in Sect. 4 (the training period is from January 1951 to April 2008).Finally, the improved prediction equation of T 1 , T 2 , SOI and EAWMI, based on the self-memorization principle, can be expressed as where (i = 0, 1, . .., 4; j = −6, −5, . .., 0).
The step-by-step forecast was performed.The retrospective order p = 6 means that the earlier data of seven observations (p + 1 = 7) should be used during the forecasting process.
The forecast results per month were saved for the next period predictions.

Long-term step-by-step forecasts of T 1 and T 2
To test the actual forecast performance of the aboveimproved model, long-term step-by-step forecasts of T 1 and T 2 from May 2008 to December 2010 for 20 months were carried out, as shown in Fig. 4. The forecast results of T 1 and T 2 were good.Within 8 months, the CCs of T 1 and T 2 were 0.9163 and 0.9187.MAPEs of T 1 and T 2 were small (only 5.86 and 6.78 %).The forecast time series from 8 months to 14 months gradually diverged, but the trend was acceptable.The CCs of T 1 and T 2 reached 0.8375 and 0.8251, and MAPEs of T 1 and T 2 were 8.32 and 9.11 %.After 14 months, the forecast began to diverge and the error started to increase, but the CCs of T 1 and T 2 remained about 0.6899 and 0.6782, and MAPEs reached 18.31 and 19.44 %, which can be acceptable.

Cross-validated retroactive hindcasts of time series T 1 and T 2
As in Sect.3, the model's skill should be further assessed by cross-validated retroactive hindcasts of the time series.Because our step-by-step forecasts need the earlier data of seven observations (p+1 = 7), we can obtain cross-validated retroactive hindcast results of T 1 and T 2 from August 1951 to December 2010, as shown in Fig. 5. From Fig. 5, the forecast performance of T 1 and T 2 was good.The CCs of T 1 and T 2 were 0.7124 and 0.7036, respectively.The MAPEs of T 1 and T 2 were small at only 19.57 and 19.79 %, respectively.The peaks and valleys of T 1 and T 2 were also forecasted accurately.The forecast results indi- Many researchers (Zhang et al., 2003b;Smith, 2004) have used the Oceanic Niño Index (ONI) which is used by the US NOAA Climate Prediction Center to determine the El Niño and La Niña years.It was defined that when the ONIs of 5 consecutive months in winter are all more than 0.5 (less than −0.5), it is an El Niño (La Niña) year.Based on the above criterion, we can divide the total 60 years (1951-2010) into three categories.It includes the 18 examples of El Niño years (such as 1958, 1964, 1966, etc.), 22 examples of La Niña years (such as 1951, 1955, 1956, etc.) 3.
From Table 3, the average CC of both T 1 and T 2 of 60 experiments within 6 months was more than 0.84, and MAPE was less than 8 %.The average CC within 12 months was more than 0.74, and MAPE was less than 12 %.According to the literature (Tofallis, 2015), when MAPE was less than 15 %, it meant the error was not great and the forecast results were good.Obviously, the forecast results of the El Niño/La Niña experiments were a little worse than those of neutral examples, which means the forecast ability of our model for the abnormal situation was a little worse than that for the normal situation.However, even for El Niño/La Niña experiments, the average CC was still more than 0.7 and MAPE was less www.ocean-sci.net/14/301/2018/Ocean Sci., 14, 301-320, 2018 than 15 %, which means the error was not too large and was still within an acceptable range.

Forecast of the SSTA field
When we obtained the forecast results of the time coefficient series T 1 and T 2 , we submitted them into the following equation to reconstruct the forecast SSTA field: where E n , T nt are the EOF space fields and forecast time coefficients, respectively, and xij is the forecast SSTA field reconstructed by EOF.
After reconstruction of the space mode (treated as constant) and time coefficient series (model prediction), the forecast of the SSTA fields was obtained based on the forecast results of T 1 and T 2 in Sect.5.2.For economy of space, we cannot draw all of the forecasted SSTA fields, so we selected a strong El Niño event (December 1997), a strong La Niña event (December 1999) and a neutral event (November 2002) as examples.
Figure 6 shows the forecast SSTA field during a strong El Niño event.From the actual SSTA field in December 1997 (Fig. 6a), an obvious warm tongue structure occurred in the area of 10 • S-5 • N, 90-150 • W in the eastern equatorial Pacific, and a warm anomalous distribution arose in the west Pacific, which indicated a weak El Niño event.The forecasted SSTA field of December 1997 is shown in Fig. 6b.Although the range of warm tongue was a litter bigger than the actual situation, the forecast shape was similar to the actual field and also the contour lines were similar.The average MAPE between the forecast field and the actual field is 8.56 %, which was controlled within 10 %.The forecast results of the improved model event were quite good for the El Niño event.
Figure 7 shows the forecasted SSTA field of a strong La Niña event.From the actual SSTA field in December 1999 (Fig. 7a), an obvious cold pool occurred in the area of 10 • S ∼ 10 • N, 120 • W ∼ 180 • W in the equatorial Pacific, which covered the Niño3.4area.This SSTA field presented a strong strength La Niña event.The forecast SSTA field from December 1999 is shown as Fig. 7b.Although the strength of the cold pool was weaker than the actual situation, the forecast shape was similar to that of the actual field.The average MAPE between the forecast field and the actual field was 9.69 %.The errors were larger than those of the El Niño event, but they can be controlled within 10 %, which is acceptable.
Figure 8 shows the forecasted SSTA field of a neutral event.From the actual SSTA field in November 2002 (Fig. 8a), a warm pool occurred in the area of 10 • S-10 • N, 120-180 • W in the equatorial Pacific, which covered the Niño3.4area.However, the warm pool was small and weak, which represented a neutral event.The forecasted SSTA field from November 2002 is shown in Fig. 8b.Comparing Figs.6-8, we can see that the forecasted SSTA field of a neutral event was a little worse than that of the El Niño and La Niña events.The forecasted shape of the SSTA field basically described the actual situation, but the warm pool in the Niño3.4area was stronger and bigger than that of the actual situation, which indicated a borderline El Niño event.The average MAPE between the forecasted field and the actual field was 14.50 %, which was big but can be accepted.
We obtained the average values of MAPE of 18 El Niño events, 22 La Niña events and 20 neutral events, which were 9.52, 9.88 and 14.67 %, respectively, representing a good SSTA field forecasting ability of our model.

Forecast of ENSO index
The ENSO index can be represented as the SSTA in the Niño3.4region (5 • N-5 • S, 120-170 • W) and the ENSO index forecast was the 3-month forecast (Barnston et al., 2012).So we also can pick up the ENSO index from the aboveforecasted SSTA field.The forecast results of the ENSO index within 20 months can also be obtained.The definition of lead time can be seen in the reference (Barnston et al., 2012).Therefore, similar to the forecast experiment in Sect.5.1,  a succession of running 3-month mean SST anomalies with respect to the climatological means for the respective prediction periods, averaged over the Niño3.4region, can be obtained, as demonstrated in Fig. 9.
The evaluation criterion of the ENSO index is the temporal correlation (TC); its definition and specific calculation steps can be seen in the literature (Kathrin et al., 2016;Nicosia et al., 2013).The TC is often used to measure the prediction  effect of the ENSO index.For example, Barnston et al. in 2012 also used the TC to compare the forecast skill of 21 real-time seasonal ENSO models.
The forecast results within lead times of 18 months are shown in Fig. 9, which demonstrate that the forecast results of the ENSO index are good.Within the lead time of 12 months, the TC was 0.8985 and the MAPE value was small at only 8.91 %.In addition, the borderline La Niña event in 2008-2009 was predicted well.After lead times of 12 months, forecasts began to diverge and the errors started to increase.Although the TC remained approximately 0.61, MAPE reached 18.58 %.Therefore, a moderatestrength El Niño event that occurred in 2009/2010 was not predicted.4.
From Table 4, the average TC of 60 experiments was 0.712, and the average MAPE was 7.62 % within 12 months for all seasons of lead time, which indicates that the overall ENSO forecast ability of our model was good.The forecast results of the El Niño examples were significantly worse than those of La Niña examples, while the forecast results of La Niña examples were significantly worse than those of neutral examples, which show the model forecast ability of the abnormal state was worse than the normal state of the ENSO index.Even for the forecast results of El Niño examples, the average TC was still above 0.6 and the average MAPE can be controlled below 10 %, which means the forecast results were still in the acceptable range.Our model not only accurately predicted the stronger El Niño and La Niña phases but also the neutral states.
The ENSO forecast often had a spring predictability barrier (Webster, 1999), which was most prominent during decades of relatively poor predictability (Balmaseda et al., 1995).To test our model, the skill should be computed over the entire time series and separately for seasonal subsets of the time series.From Table 4, we can see that although the forecast results of the present model in the spring were worse than in the autumn, the margin was not high, which means  the model can overcome the "spring predictability barrier" to some extent.

Compared with six mature models
Barnston et al. ( 2012) compared many ENSO forecast models.Based on his research, we selected four high-quality dynamical models, including ECMWF, JMA, the National Aeronautics and Space Administration Global Modeling and Assimilation Office (NASA GMAO) and the National Centers for Environmental Prediction Climate Forecast System (NCEP CFS; version 1).Two high-quality statistical models also were selected, including the University of California, Los Angeles Theoretical Climate Dynamics (UCLA-TCD) multilevel regression model and the NOAA/NCEP/CPC constructed analogue (CA) model.The details of the above models can be found in these references (Reynolds et al., 2002;Luo et al., 2005;Barnston et al., 2012).
We then compared the forecast ability of the above six models with that of our model.All of the experiments of our model and six other models were conducted under the same conditions using the same historical data for modeling and the same initial values to forecast.On the CPC website, there are detailed explanations of the six models' training samples and the initial values.So we do not need to install all these models on their own machines and run them for forecasting.We just made training samples, and initial values of our model were the same as those of the six selected models.At an 8-month lead time, the TC of our model for all seasons combined was 0.613 (Fig. 10).In brief, the forecast ability of the ECMWF model was slightly better than that of our model but the ability of the other five models was worse than that of our model.However, in regard to the forecast length, the TC within 12 months of our model is greater than 0.6, which was superior to the ECMWF model.In addition, the forecast results of the UCLA-TCD model and the CPC CA model reduced quickly after 5-month lead times, so the forecast ability of our model was more stable than theirs.
The root mean square error (RMSE) was also examined to assess the performance of discrimination and calibration.Barnston et al. (2012) believed that all seasonal RMSE values contributed equally to a seasonally combined RMSE.So we drew Fig. 11 to show seasonally combined RMSE.
www.ocean-sci.net/14/301/2018/Ocean Sci., 14, 301-320, 2018 From Figs. 10 and 11, we can see the highest correlation tends to have lower RMSE.So the RMSE of our model was slightly higher than that of the ECMWF model, but it was much lower than those of the other five models.Figures 10  and 11 show the average TC and RMSE of the 240 experiments compared with six mature models and cover a variety of different types of ENSO and different lead times.So those samples should be really representative.

Conclusions
A new forecasting model of the SSTA field was proposed based on a dynamical system reconstruction idea and the principle of self-memorization.The approach of the present paper consisted of the following steps.
The SST field can be time (coefficients)-space (structure) deconstructed using the EOF method.Take T 1 , T 2 , SOI and EAWMI and consider them as trajectories of a set of four coupled quadratic differential equations based on the dynamical system reconstruction idea.The parameters of this dynamical model were estimated using a GA.
The forecast results of the dynamical model can be improved by the self-memorization principle.The memory coefficients in the improved self-memorization model were obtained using the GA method.
The long-term step-by-step forecast results and crossvalidated retroactive hindcast results of time series T 1 and T 2 are all found to be good, with a CC of approximately 0.80 and a MAPE of less than 15 %.
The improved model was used to forecast the SSTA field.The forecasted SSTA fields of three types of events are accurate.Not only is the forecast shape similar to the actual field but also the contour lines are similar.
The improved model was also used to forecast the ENSO index.The average TC of 60 examples within 12 months is 0.712, and the MAPE value is small at only 7.62 %, which proves that the improved model has better forecasting results of the ENSO index.Although the forecast results of the model in the summer were worse than in the winter, the margin was not high, which means that the model can overcome the spring predictability barrier to some extent.Finally, compared with the six mature models, the new dynamicalstatistical forecasting model has a scientific significance and practical value for the SST in the eastern equatorial Pacific and El Niño/La Niña event predictions.

Discussion
L' Heureux et al. (2013) reported that using different data sets and time periods, the second EOF is not stable, entirely due to the strong trend.So we need to do more experiments to prove that we choose the second mode of EOF to be appropriate and whether different time periods will make our forecast unstable or not.Our original data are the monthly average SST data from January 1951 to December 2010 (60 years).We will increase the length of the data for 20 years (January 1931-December 2010) and for 10 years (January 1941-December 2010), and decrease the length of the data for 10 years (January 1961-December 2010) and for 20 years (January 1971-December 2010).Then, we will use the same method to reconstruct a model and forecast the ENSO index as in Sect.5.4.The results show that, in the 60 experiments, the difference among forecast results of both TC and MAPE of five different sample data is lower, and no abnormal changes that are suddenly worse or better appear.All this indicates that using different data sets and time periods may have a certain impact on the pattern of the second EOF, but the impact on our forecast is not great and it will not make our forecast unstable.
Actually, the amount of variables and which variables are used in our model become key issues to be resolved.We have developed a complex coupled model of four-factor differential equations, so we are more concerned with the correlations between each of them.The correlation must be considered as an important criterion to select the factors, but in order to further verify the correctness of the selection criterion, we have carried out the prediction experiments (the 60 crossvalidated retroactive hindcast experiments of the ENSO index for all seasons combined at lead times of 8 months) of different variables.
We can see that for all the forecast results of the models of different variables, the prediction results of T 1 , T 2 and SOI are the best among the three factors, and the prediction results of T 1 , T 2 , SOI and EAWMI are the best among the four factors.However, the prediction results of T 1 , T 2 , SOI and EAWMI are the best among all the factors, which proves that our selection factors are correct.In our previous study (Hong et al., 2015), the model of the western Pacific subtropical high was established by using the correlations as a criterion to select factors, and their forecast results are also good.Now, we use the correlations as a criterion to select factors also in line with our previous research.
Based on the definition of overfitting and the previous studies (Golbraikh et al., 2003;Everitt and Skrondal, 2010), there is no evidence that more parameters will result in overfitting.We can judge whether a model is overfitting or not by the accuracy of prediction results of independent samples (Golbraikh and Tropsha, 2002;Qin and Li, 2006).
In the sample training, our model does not purposely pursue the high degree of the training sample fitting and improve the effectiveness of the independent generalization.In fact, in our paper, the forecast results of the cross-validated retroactive hindcasts (Sect.5.2) and the independent sample validation (Tables 3 and 4) are both good.Especially in the independent sample validation of the ENSO index (Table 4), we have carried out the 240 independent sample validation predictions of four seasons of different ENSO events, and the coverage of independent samples test is very wide.More-over, compared with six mature prediction models, the forecast results of our model are also good, which proves that the overfitting problem does not exist in our model.So according to the definition of overfitting, we can say the overfitting phenomenon does not exist in our model.
Compared with the original model, the reasons why the improved model has good forecast results and can overcome the spring predictability barrier to some extent are as follows.
Recently, many studies have pointed out that spring is the most unstable season of the air-sea interaction and the error is likely to develop or grow in the spring, resulting in the spring predictability barrier (Zhang et al., 2012;Philander et al., 1992).When the original model uses the indexes in summer as the initial values to predict, the SOI factor representing the air-sea interaction is most unstable in the spring and the EMWMI factor does not have much influence on ENSO in summer, so the forecast results using the indexes in summer as the initial values are certainly much worse than those using the indexes in the winter as the initial values.That is why our original model does not overcome the spring predictability barrier.
However, the introduction of the self-memorization dynamics principle can help our model overcome the spring predictability barrier to some extent.Although the lead time is still summer (such as JJA), the information of the initial value actually contains the previous p + 1 month (in this case, p = 6), which contains the information of the previous 7 months, including the information of the T 1 , T 2 , SOI and EMWMI factors in winter (January, February), spring (March, April, May) and summer (June and July).From the dynamical analysis, in this situation, the information and interaction relationship of four factors has been accumulating for a long period (from winter to summer), containing many air-sea interaction processes, and the winter monsoon contains abnormal information, so the forecast results of our improved model will be much better than the original model which simply uses only one initial value.That is why the improved model overcomes the spring predictability barrier to some extent.
The forecast results of our model are good, but it still has some problems: 1.The inclusion of these terms and the physical processes these terms represent in Eq. ( 2) are important, especially for the discussion of dynamical characteristics of the dynamical model.However, now it is difficult to give a clear meaning.Now, the main work of our paper is the prediction experiments of the model.For the reasons of time and length, this paper mainly discusses the prediction results of the model.The physical processes do these terms represent and the discussion of the dynamical characteristics of the model will be the focus of our next work.Before this, we have also used Takens' delay embedding theorem to reconstruct the dynamical model of the western Pacific subtropical high (WPSH).
Based on the reconstructed dynamical model, dynamical characteristics of WPSH are analyzed and an aberrance mechanism is developed, in which the external forcings resulting in the WPSH anomalies are explored, which have been published (Hong et al., 2016).We also study the bifurcation and catastrophe of the west Pacific subtropical high ridge index of a nonlinear model (Hong et al., 2017).Based on our previous method and work, our next work is to analyze the physical processes and the dynamical characteristics of the SST field.
2. The experiments in the present study have proven that the forecasting results of the improved model are good for large-scale systems, such as ENSO events, and the forecasting period has been extended.However, for small-scale systems, such as hurricanes, whether the forecast results could be improved using the present improved model needs to be further verified.
3. Our paper focuses primarily on these defined indices with T 1 , T 2 to reconstruct a prediction model.Maybe we can select variables (predictors) based on EOF analysis and our model may be a more physically oriented model.Maybe we can learn from Yim et al. (2013Yim et al. ( , 2015) ) to draw correlation maps between these fields and the SSTA field, and select the predictors from physical considerations.All the above questions require a lot of experiments to be carried out.
These items will be our future work.
The others data of our paper all can be obtained from the website and the persistent URLs is https://www.esrl.noaa.gov/psd/data/(Kalnay et al., 1996).
Appendix A: The principle of dynamical model reconstruction Suppose that the physical law of a nonlinear system going over time can be expressed as the following difference form: where f i is the generalized nonlinear function of q 1 , q 2 , . .., q i , . .., q N , N is the number of variables, and M is the length of observed data.f i (q j t 1 , q j t 2 , . .., q j t i , . .., q j t N ) can be assumed to contain two parts: G j k representing the expanding items which contain variable q i , and P ik just representing the corresponding parameters which are real numbers (i = 1, 2, . .., N, j = 1, 2, . .., M, k = 1, 2, . .., K).
It can be supposed as follows: Parameters of the above equation can be determined through inverting the observed data.Vector P , which satisfies the above equation, can be solved based on a given vector D.
Assuming q is unknown, it is a nonlinear system.However, assuming P is unknown, it is a linear system.With the restriction S = (D − GP ) T (D − GP ) as a minimum, GA is introduced as an optimization solution search in the model parameter space.
Assuming that the parameter matrix P is the population (solutions), the S = (D − GP ) T (D − GP ) is an objective function, l i = 1 S is the value of individual fitness, and L = n i=1 l i is the value of total fitness.The operating steps of GA include creation and coding of initial population (solutions), fitness calculation, the choice of male parents, crossover and variation, etc.A detailed theoretical explanation can be obtained from Wang (2001).The step length is 1 month during the calculation.After optimization searches and genetic operations, the target value can be rapidly converged and each optimal parameter of the dynamical equations can be obtained.
Through the above approach, we can obtain parameters of a nonlinear dynamical system and reconstruct the nonlinear dynamical equations from observed data.

Appendix B: The mathematical principle of self-memorization dynamics of systems
The dynamical equations of a system can be expressed as ∂x i ∂t = F i (x, λ, t) i = 1, 2, . .., J, (B1) where J is an integer, x i is the ith variable of the system state, and λ is the parameter.Equation (B1) represents the relationship between a source function F and a local change of x.Obviously, x is a scalar function with time t and space r 0 .A set of time T = [t −p . ..t 0 . ..t q ] can be considered, where t 0 is an initial time.A set of space R = [r a . ..r i . ..r β ] can be considered, where r i is a spatial point.An inner product in space L 2 :T × R is defined by Accordingly, a norm can be defined as For a completion L 2 , it can become a Hilbert space H .A generalized one in H can be regarded as a solution of the multitime model.By introducing a memorization function β(r, t), we can obtain where r in β(r, t) can be dropped through fixing on the spatial point r 0 .Supposing that function β(r, t) and variable x, etc. are all continuous, differentiable and integrable, an integration by the left parts of Eq. (B3) can be made as As a matter of convenience, we set β t ≡ β(t), β 0 ≡ β(t 0 ), x t ≡ x(t), x 0 ≡ x(t 0 ); the following text uses similar notations.Then, Eq. (B7) can be expressed as S 1 is called a self-memory term and S 2 is called an exogenous effect term.
For the convenience of calculations, the above selfmemorization equation can be discretized.The differential by difference and the summation can replace the integration in Eq. (B9) and the mean of two values which are at adjoining times; i.e., x m i ≈ 1 2 (x i+1 + x i ) ≡ y i can simply replace x m i .Taking an equal time interval t i = t i+1 − t i = 1 and incorporating β i and β t , we can obtain a discretized selfmemorization equation as follows: where F is the dynamical kernel of the self-memorization equation, α i = (β i+1 −β i ) β t ; θ i = β i β t .Based on Eq. (B10), the above technique performed computations and the forecast can be called a self-memorization principle.

Figure 1 .
Figure 1.(a, c) First and second modes of the EOF deconstruction of the SSTA field and (b, d) the corresponding PC time series.

Figure 2 .
Figure 2. Forecast results of the first time coefficient series (a) and the second time coefficient series (b) of the SSTA field by the original model.

Figure 3 .
Figure 3.The cross-validated retroactive hindcast results of the first time coefficient series (a) and the second time coefficient series (b) of the SSTA field by the original model.

5
Model prediction experiments 5.1 Forecast of time series T 1 and T 2

Figure 4 .
Figure 4. Long-term step-by-step forecast results of the first time coefficient series (a) and the second time coefficient series (b) of the SSTA field by the improved model.

Figure 5 .
Figure 5.The cross-validated retroactive hindcast results of the first time coefficient series (a) and the second time coefficient series (b) of the SSTA field by the improved model.
and the remaining 20 experiments of neutral years.Since the details in Fig. 5 are not clear, we list the forecast results of 60 experiments (including 18 El Niño examples, 22 La Niña examples and 20 neutral examples) in Table

Figure 6 .
Figure 6.The forecast SSTA field (a) and the actual SSTA field (b) of an El Niño event (December 1997).

Figure 7 .
Figure 7.The forecast SSTA field (a) and the actual SSTA field (b) of a La Niña event (December 1999).

Figure 8 .
Figure 8.The forecast SSTA field (a) and the actual SSTA field (b) of a neutral event (November 2002).

Figure 9 .
Figure 9.The improved dynamical-statistical model prediction of the ENSO index.
We should give more examples to test the ENSO prediction ability of our model.As in Sect.5.3, we can divide 60 examples into three types, which are examples of the El Niño year, La Niña year and neutral year.Finally, we can obtain the forecast results of different types of examples in different lead times, as shown in Table

Figure 10 .
Figure 10.Temporal correlation between model forecasts and observations for all seasons combined, as a function of lead time.Each line highlights one model.
Fig .11. RMSE in standardized units, as a function of lead time for all seasons combined.Each line Figure 11.RMSE in standardized units, as a function of lead time for all seasons combined.Each line highlights one model.

Table 1 .
The correlation analysis between the front two time series (T 1 , T 2 ) and nine impact factors.

Table 2 .
The CCs and MAPEs of the long-term fitting test when the retrospective order p is different.

Table 3 .
The forecast results of T 1 and T 2 in different examples within 6 and 12 months.

Table 4 .
The TC and the MAPE between model forecasts and observations within 12 months for Nov-Jan, Dec-Feb and Jan-Mar as lead times for winter, for Feb-Apr, Mar-May and Apr-Jun as lead times for spring, for May-Jul, Jun-Aug and Jul-Sep as lead times for summer and for Aug-Oct, Sep-Nov and Oct-Dec as lead times for autumn.
= ∂β(t)/∂t.The mean value theorem can be introduced into the third term in Eq. (B4), and the following equation can be obtained: Because the x value, which is at initial time t 0 and middle time t m , only on the fixed point r 0 itself, relates to the first term and the second term in Eq. (B6), they are called selfmemory terms.Also, we can call the third term an exogenous effect, i.e., which is contributed by other spatial points.Similarly to Eq. (B4), for multi-time t i , i = −p, −p + 1. .., t 0 , t, it gives After the same term β(t i )x(t i ), i = −p + 1, −p + 2, ..., 0 is eliminated, we haveβ(t)x(t) − β(t −p )x(t −p ) − m (t 0 ) ≡ x(t m ), t 0 < t m < t.Substituting Eqs.(B4) and (B5) in Eq. (B3) and carrying out an algebraic operation, the following equation can be obtained: