Improved ENSO Prediction Skill Resulting From Reduced Climate Drift in IAP‐DecPreS: A Comparison of Full‐Field and Anomaly Initializations

When initiated from observational conditions, coupled climate models used in seasonal predictions generally experience climate drifts. How climate drift affects El Niño–Southern Oscillation (ENSO) prediction is important but not clearly understood. Here, we investigate this issue by comparing seasonal hindcasts using two distinct initialization approaches, namely, anomaly and full‐field initializations, based on a climate prediction system named IAP‐DecPreS. The differences between the two approaches are mainly evident in the drift behavior. We find that the hindcasts based on anomaly initialization (Hindcast‐A) have higher ENSO prediction skill compared to those based on full‐field initialization (Hindcast‐F). The climate drifts are largely reduced in the Hindcast‐A as expected. In contrast, the Hindcast‐F features a growing warming of the equatorial central eastern Pacific with increasing lead times. To investigate the impact of drift on the prediction, the 1997/1998 and 2015/2016 El Niño cases are analyzed. At a 7‐month lead, the Hindcast‐A reasonably predicts the two events, while the Hindcast‐F shows large errors in both evolution and amplitude. Budget analyses show that the underestimation of warming tendency in the Hindcast‐F is caused by cooling effects of excessive anomalous surface shortwave radiative flux and anomalous temperature advection by mean horizontal currents, both of which are associated with the climate drift. Our results imply that the use of the AI scheme can improve ENSO predictions through a reduction in climate drift in IAP‐DecPreS. The drifts can dynamically influence ENSO predictions, and their impact cannot be thoroughly removed via the empirical bias correction. Thus, reducing drift impacts is necessary.


Introduction
El Niño-Southern Oscillation (ENSO) is the strongest interannual climate fluctuation on Earth and has a striking impact on global climate through atmospheric and oceanic teleconnections (Philander, 1983;Aceituno, 1992;Alexander et al., 2002;McPhaden et al., 2006;Wu et al., 2017aWu et al., , 2017bYeh et al., 2018). Because of its far-reaching influences, ENSO has been treated as an important source of seasonal climate predictability (Kim et al., 2012;Latif et al., 1998;Peng et al., 2009). Reliable ENSO prediction can aid in reducing its socioeconomic impacts .
Because of the improvements in understanding ENSO mechanisms, modeling techniques, observations, and data assimilation techniques over the past decades, significant progresses have been made in ENSO predictions Cane et al., 1986;Jin et al., 2008;Latif et al., 1998;Tippett et al., 2017). The ENSO predictive skill by coupled general circulation models can be obtained at a lead time of 2-3 seasons (Barnston et al., 2012Jin et al., 2008;Wang et al., 2009) and even greater than 1 year (Luo et al., 2008(Luo et al., , 2017. Although it was suggested that today's best ENSO forecast systems might already be close to the predictability limit in the cold tongue region (Newman & Sardeshmukh, 2017). The satisfactory predictions of ENSO's amplitude at long lead time are still challenging (Barnston et al., 1999McPhaden et al., 2006;Tang et al., 2018).
Generally, there are four factors limiting the ENSO predictive skill: inherent predictability limits, gaps in observational systems, model flaws, and suboptimal available observational data use (Chen, 2010). An important model flaw is model systematic biases, which represent the difference between the model simulated and the observed climatology (Magnusson, Alonso-Balmaseda, Corti, et al., 2012;Wang et al., 2014). Model systematic biases, particularly the "cold-tongue" bias in the tropical Pacific, greatly influence ENSO simulations (Guilyardi, 2006;An et al., 2010;Chen et al., 2016;Kim, Cai, Jin, & Yu, 2014;. Specifically, the underestimated negative shortwave cloud radiative feedback in response to ENSO-related SSTAs increases the simulated ENSO amplitude, which is a common bias among the models from the CMIP3 to the CMIP5, and related to the coldtongue bias (Chen et al., 2013;Chen et al., 2018;Chen, Sun, et al., 2019). To reduce the influence of systematic biases on prediction, the full-field initialization (FFI) approach is a common approach in operational prediction systems that replaces the initial model state with the best possible available estimate of the observed state. FFI can efficiently improve the climatology at the beginning of prediction in the coupled models. Nevertheless, because of model deficiencies, there is still a distance between the initial conditions and model inherent dynamics. As the prediction continues, its trajectory tends to drift away from the observations and back to its preferred state, no matter how small the initial error. This type of tendency is known as "model drift" or "initial shock" Gupta et al., 2013;Smith et al., 2012;Sanchez-Gomez et al., 2016). The drifts aggravate the model mean state errors and may contaminate the forecast signal, making it difficult to extract the predicted climate variability (Smith et al., 2007).
Several solutions have been presented to overcome the problem of model drift, such as empirical bias correction, flux correction, and AI. Empirical bias correction uses retrospective predictions to remove the lead time-dependent biases caused by the model drift from the raw predictions (Stockdale, 1997;ICPO 2011). Flux correction, such as heat flux and momentum flux correction, is applied in the coupling interface between the atmosphere and the oceanic components by adding an additional correction term during model integrations (Manganello & Huang, 2009;Magnusson, Alonso-Balmaseda, Corti, et al., 2012). AI assimilates the observed anomalies added to the model climatology (Pierce et al., 2004;Schneider et al., 1999;Smith et al., 2007). In theory, the forecasts do not drift (Smith et al., 2013). Although the aforementioned solutions can partly reduce the impacts of climate drift impacts, weaknesses remain evident. The empirical bias correction strategy neglects nonlinear interactions between the model drift and predictive anomalies. Hence, residual biases still emerge in the predictions (Kirtman et al., 2013). Utilizing the flux correction may result in unreal physical processes in the balanced model physics (Neelin & Dijkstra, 1995). For the AI, because the model systematic biases are fully developed, and the observational anomalies are often not consistent with those of the model, it still requires bias correction (Magnusson, Alonso-Balmaseda, Corti, et al., 2012). Therefore, model drift remains an unsolved problem in climate predictions (Gupta et al., 2013).
How model drifts interact with the predicted variability and the prediction bias generation mechanisms has been a focus of the World Climate Research Group Working Group on Subseasonal to Interdecadal Prediction (WGSIP; https://www.wcrp-climate.org/wgsip-projects/lrftip). Previous studies have found that there is a close relationship between model drifts and ENSO predictions. For example, based on the CMIP3/Clipas and ENSEMBLES, the bias of the simulated tropical annual cycle, which intensifies with increasing leads, is related to the predictive ENSO variability (Jin et al., 2008). By comparing the impacts of three different initialization schemes on the prediction skill in ECMWF System3, the most skillful scheme for the prediction of the Niño 4 region is the scheme that has the minimal model drifts in this region (Balmaseda & Anderson, 2009). Based on the DEMETER multimodel hindcasts, there is a robust inverse relation between skill and model drift, particularly for the tropical Pacific (DelSole & Shukla, 2010). The influence of three different initialization strategies on ENSO predictions are compared in the ECMWF IFS model. The highest skill at seasonal time scale is obtained using momentum flux correction because it avoids the positive feedback responsible for a strong cold bias and thus reduces the model drift in the tropical Pacific (Magnusson, Alonso-Balmaseda, Corti, et al., 2012). However, the detailed physical processes involved in the interactions between the model drifts and the predicted variability remains unknown.
Recently, a process-based study found that the growing mean bias of the equatorial Pacific upper oceanic temperature in an operational seasonal forecast model is a possible cause of the growing forecast errors in ENSO amplitude through modulation of the thermocline feedback (Kim et al., 2017). The growing forecast errors can result from either the model drift or the intrinsic model physics. The relative contribution of the model drift to the forecast errors is unclear. Hence, understanding the physical processes involved in the interactions between drift and ENSO prediction is of urgent need. Here, we address this issue by comparing the seasonal hindcasts using two distinct initialization approaches, AI and FFI. The hindcasts based on the FFIs cannot avoid model drifts, while the hindcasts based on AIs can reduce the drifts under the condition of balanced model physics unaffected by the artificial additional flux term. We address the following three questions: (1) How well does the IAP-DecPreS perform in ENSO prediction based on the two distinct initialization approaches? (2) What are the differences in the drift behaviors between the two hindcast types? (3) What are the fundamental physical processes involved in the interactions between the model drifts and predicted variability?
The remainder of the paper is organized as follows. Section 2 introduces the prediction system and analysis methods. In section 3, we first show the model performance in the simulation of the mean state of the tropical Pacific and ENSO in sections 3.1 and 3.2. Section 3.3 presents the ENSO predictive skill results in seasonal hindcasts based on full-field and AIs. The causes of these differences are discussed in section 3.4. Finally, discussion and summary of the results are provided in sections 4 and 5.

Climatic Prediction System and Experiments
IAP-DecPreS is a near-term climate prediction system based on a state-of-the-art fully coupled global climate model, FGOALS-s2, which was developed by the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG) at the Institute of Atmospheric Physics (IAP), Chinese Academy of Sciences (CAS) (Bao et al., 2013, Zhou, Yu, et al., 2014. In FGOALS-s2, the atmospheric component is the Spectral Atmosphere Model in IAP LASG version 2 (SAMIL2), with a horizontal resolution of~2.81°(longitude) × 1.66°(latitude) and 26 vertical levels extending from the surface to 2.19 hPa (Bao et al., 2013). The oceanic component is the LASG IAP Climate System Ocean Model version 2 (LICOM2), with 30 levels in the vertical direction. The zonal resolution of LICOM2 is 1°. To better resolve the equatorial waves, the meridional resolution is 1°outside 20°S/°N, refining from 1°to 0.5°between 10°a nd 20°S/°N, and turning into 0.5°between 10°S and 10°N (Lin et al., 2016;Liu et al., 2012). The land component is the Community Land Model version 3 (CLM3, Oleson et al., 2004). The sea ice component is the Community Sea Ice Model version 5 (CSIM5, Briegleb et al., 2004). FGOALS-s2 has been used in the Coupled Model Intercomparison Project Phase 5 (CMIP5) experiments (Zhou, Chen, et al., 2014;Zhou et al., 2018).
FGOALS-s2 is initialized through a weakly coupled assimilation, in which observational ocean temperatures are assimilated into the ocean component of the coupled model using a newly designed assimilation scheme, ensemble optimal interpolation (EnOI), and incremental analysis update, referred to as EnOI-IAU (Wu et al., 2018;Hu et al., 2019). The assimilated observational data include observational sea surface temperatures (SSTs) from HadISST version 1.1 (Rayner et al., 2003) and upper 700 m oceanic temperature profiles from the EN4 data set version 1.1 (Good et al., 2013). Please see Wu et al. (2018) for detailed information regarding EnOI-IAU and IAP-DecPreS.
We used two different initialization approaches, FFI and AI. The FFI assimilates raw observational ocean temperature, while the AI uses the constructed ocean temperature, which is the sum of the observational ocean temperature anomalies and the model climatology. The anomalies of the ocean temperature profiles are calculated with respect to the climatology of the corresponding objective analyses from the EN4 data set for the period of 1950-2005. We have three runs for both the AI and the FFI, which are started from three historical runs.
An ensemble of 14-month hindcasts was conducted based on the initial conditions offered by the FFI (AI) runs that started from February, May, August, and November in each year over the period 1979-2016. These hindcasts are referred to as Hindcast-F and Hindcast-A, corresponding to the FFI and AI, respectively. For each hindcast, both Hindcast-F and Hindcast-A have nine members, which started from three initialization runs with initial perturbations. More details of the experiment's designs are shown in Table 1.

Observational Data
The following observational and reanalysis data sets were used to assess the predictive skill: (1) observational precipitation data from the Global Precipitation Climatology Project (GPCP, Adler et al., 2003), 2) observational SST from the Met Office Hadley Centre Sea Ice and SST Data set (HadISST version 1.1, Rayner et al., 2003), (3) circulation and tropospheric temperature data from the ERA-interim reanalysis (Dee et al., 2011), (4) the monthly ocean temperature and current were provided by the NCEP Global Ocean Data Assimilation System (Saha et al., 2006) produced at the National Centers for Environmental Prediction (NCEP), and (5) the surface heat fluxes are from TropFlux.v1 (Praveen Kumar et al., 2012). All the data sets are monthly mean and cover the period 1979-2016.

Analytical Methods
The Niño 3.4 index is defined as the area-averaged SSTAs in the equatorial central and eastern Pacific (5°S to 5°N, 170-120°W). We used the anomaly correlation coefficient (ACC) to evaluate the deterministic prediction skill, that is, the ensemble mean prediction. The significance level of the ACC is tested by the Students t test.
The predicted anomalies are calculated with respect to the predicted climatology during the period 1979-2016. By using the anomalies, mean state biases of the prediction for each start month and lead time are approximately removed, equivalent to empirical bias correction. The lead time of n months of the prediction represents the time interval between the predictive target and initial date.
The model drift is defined as the change in the predictive climatological state deviated from the model intrinsic climatological state with the lead time during a prediction (Sanchez-Gomez et al., 2016), which is quantified as follows: where X p (i,j) is the climatological state of prediction initiated from date i (4 times for this study) for target time j, estimated as the average of the raw output data by using all hindcast experiments. X m (j) is the corresponding model intrinsic climatological state, which can be estimated based on historical simulations. X d (i,j) represents the model drift, which is the difference between X p (i,j) and X m (j). As in the historical (RCP4.5) run, before (after) 2005 Hindcast-F As in Hindcast-A Initial states produced by FFI 9 As in the historical (RCP4.5) run, before (after) 2005 To understand the physical processes responsible for the ENSO prediction in the hindcast runs, the mixed layer heat budget is diagnosed for both the reanalysis and Hindcast runs. The mixed layer ocean temperature equation (Chen et al., 2015;Li et al., 2002;Ray et al., 2018aRay et al., , 2018b is written as follows: where Q ′ and Res are the net surface heat flux and the residual term, respectively, and D ′ O is the oceanic temperature transport effect because of 3-D advections. D ′ O can be further decomposed into three components as follows: where T is the mixed layer temperature and (u,v,ω) are the zonal, meridional and vertical velocity, respectively. The primes (bars) represent the anomaly (climatology). The temperature gradient is calculated using centered finite difference operation.
The net surface heat flux converted to temperature units is equal to the following: The net heat flux consists of the surface latent heat (Q ′ LH ), sensible heat (Q ′ SH ), longwave radiation (Q ′ LW ), and shortwave radiation (Q ′ SW ) fluxes. The parameters ρ and C P are the seawater density and specific heat at a constant pressure; h is the mixed layer depth, which is defined as a depth at which ocean temperature is 0.5°C lower than that at the surface level. We also tried to use 0.75°C and 1°C as the threshold. The results here are not sensitive to the definition of the mixed layer depth.
The entrainment term in (2) arises from the changes in the mixed layer depth (Ray et al., 2018a(Ray et al., , 2018b, which is Where the T −h is the temperature at the layer base, and T D is the vertical-averaged mixed layer temperature. The subgrid-scale oceanic mixing, and other additional processes such as nonlocal vertical mixing, submesoscale mixing, and sensible heating due to rainfall and rivers are involved in the residual term (Ray et al., 2018a(Ray et al., , 2018b).

Mean State Simulation
To evaluate the model performance in simulating the climate mean states of the tropical Pacific, we compare the annual mean SST and precipitation from the FGOALS-s2 historical runs to the observations (Figure 1). For the observations, there are remarkable north-south and west-east asymmetries for both the precipitation and SST (Figures 1a and 1d). For FGOALS-s2, although the asymmetries are reasonably represented in the simulation, there are evident biases in the spatial pattern and intensity (Figures 1b and 1e). A cold bias is present over the extratropical western Pacific, while the southeastern and northeastern Pacific show a warm bias. Over the equatorial Pacific, the cold tongue is narrow, equatorially confined, and associated with a cold SST bias greater than -0.8°C (Figure 1c).
Along with the SST bias, the precipitation bias features a "double Intertropical Convergence Zone" (Figure 1f). These features of the model systematic errors are consistent with the multimodel mean of the CMIP5 models Wang & Li, 2017), suggesting that the results based on FGOALS-s2 have a certain degree of representativeness.

ENSO Simulated by the Historical Runs
The ENSO-related SSTAs in the tropical Pacific are reproduced by the FGOALS-s2 historical runs (Figures 2a and 2b). The major deficiencies are that simulated warm SSTAs in the equatorial central eastern Pacific are stronger than that in the observation, and its center is shifted westward relative to the latter. In the observation, the warm SSTAs enhance convection over the equatorial central Pacific, which in turn decreases downward shortwave radiative flux and thus dampens the underlying warm SSTAs (Figure 2c). The FGOLAS-s2 reproduces the positive precipitation anomalies but greatly underestimate the damping effect of the deep-convection-induced negative downward shortwave radiative flux anomalies (Figure 2d). This is closely associated with the cold tongue bias in the model, which is a common and long-lasting problem from the CMIP3 to CMIP5 (Chen et al., 2013;Chen et al., 2018).
ENSO is a combination of a variety of positive and negative feedbacks (McPhaden et al., 2006). To evaluate the model performances in simulating these processes, we compare the mixed layer budgets in the equatorial central eastern Pacific during ENSO developing phase between the reanalysis and FGOALS-s2 (Figure 3). Though the model simulates the three major positive feedbacks, Bjerknes SST-wind-thermocline feedback (−ω∂T ′ =∂z), zonal advective feedback (−u ′ ∂T=∂x), and the Ekman pumping-induced anomalous upwelling feedback (−ω ′ ∂T=∂z), their magnitudes are smaller than the counterparts in the reanalysis, which is a common problem in many CMIP5 models (Bellenger et al., 2014;Kim, Cai, Jin, Santoso, et al., 2014). For the two major damping effects during the ENSO developing phases, the model reproduces the negative latent heat flux anomalies associated with Newtonian cooling effect, but fails to simulate negative shortwave radiative flux anomalies induced by the enhanced deep convection. The net effects of weaker positive feedbacks and weaker negative feedbacks cause the growth rate of ENSO in the FGOALS-s2 is close to that in the observation (Figures 2a and 2b).
These results indicate that the model can reproduce basic features of the ENSO, consistent with Bellenger et al. (2014), which suggested that the FGOALS-s2 is one of the best CMIP5 models in the ENSO simulations

Journal of Advances in Modeling Earth Systems
in terms of amplitude and seasonality. However, it is worth noting that the detailed feedback strengths for the ENSO development in the model are not consistent with those in the reanalysis.

Predictive Skill for ENSO
To compare the predictive skill for the wintertime tropical Pacific SST variability, we calculate the anomaly correlation coefficient (ACC) skill scores of the Hindcast-A and the Hindcast-F (Figure 4). At a lead time of 1 month, the highest predictive skills are in the tropical central eastern Pacific for both hindcasts, with ACC values greater than 0.7 (Figures 4a and 4e), suggesting that both hindcasts can skillfully predict ENSOrelated wintertime SSTAs when initialized from November. With increasing lead times, the ACC in Hindcast-F more rapidly decreases over the tropical central eastern Pacific than that in Hindcast-A (Figures 4b-4d and 4f-4h). Specifically, over the tropical central Pacific, the ACC values are higher than 0.5 at a 10-month lead in Hindcast-A, which is significant at the 1% level according to Student's t test. In contrast, the skill values of Hindcast-F are smaller or even negative. The results suggest that Hindcast-A has higher skills in SST variability prediction at long lead times.
To examine the ENSO prediction skill, the time series of the Niño 3.4 index derived from the two hindcasts are compared to the observations. As shown in Figure 5, the observed Niño 3.4 indices are well predicted at a 1-month lead time. The correlation coefficients between the hindcasts and observations are both 0.86 for Hindcast-A and Hindcast-F ( Figure 5a). As the lead time increases, the prediction biases tend to increase (Figures 5b and 5c). Correspondingly, the correlation coefficients are 0.71, 0.60, and 0.45 (0.61, 0.36, and 0.24) in Hindcast-A (Hindcast-F) at 4-, 7-and 10-month lead times, respectively. Specifically, the 1982Specifically, the /1983Specifically, the , 1997Specifically, the /1998Specifically, the , 2002Specifically, the /2003Specifically, the , 2006Specifically, the /2007Specifically, the , and 2015Specifically, the /2016 El Niño events are predicted at a 7-month lead, and the 1988/1989, 1999/2000, and 2007/2008 La Niña are predicted even at 10month lead times by Hindcast-A.
Among the samples, three events including 1982/1983, 1997/1998, and 2015/2016 are termed extreme El Niño events, which are characterized by an exceptional warming in the eastern equatorial Pacific, leading to an equatorward shift of the Intertropical Convergence Zone, and hence intense rainfall in the equatorial eastern Pacific where cold and dry conditions normally prevail . The skillful predictions of these events are of great significance, but reasonable prediction of the ENSO amplitude at long lead times remains a challenge (Barnston et al., 1999Tang et al., 2018). Notably, Hindcast-A and Hindcast-F markedly differ in their predictions for the amplitude of the three extreme El Niño events at a 7-month lead time (Figure 5c). Compared to Hindcast-F, Hindcast-A can successfully predict the onset and development of these three events, and hence their amplitudes are reasonably reproduced ( Figure 6). In contrast, the growth rate of the ENSO intensity in Hindcast-F is much weaker than that in Hindcast-A and the observation. Hence, the ENSO intensity is too weak in Hindcast-F during its entire life cycle. Notably, the individual members show a larger spread in Hindcast-A compared to Hindcast-F for the 1997/1998 and 2015/2016 events after integrating for 5 months, exceeding 95% confidence level according to the F test, which implies that the larger uncertainty emerges in the prediction in Hindcast-A. For the 1982/1983 episode, the spread of the individual members increases in Hindcast-F. Meanwhile, the ensemble mean prediction shows a weak La Niña condition, indicating failure (Figure 6c).
In summary, the aforementioned results demonstrate that the ENSO prediction skill in Hindcast-A is higher than that in Hindcast-F, particularly at long lead times.

Factors Causing the Differences in Predictive Skill
The predictive skill can be influenced by (1) initial conditions, (2) model drifts, and (3) ENSO dynamic physical processes. To investigate the difference in predictive skill between Hindcast-A and Hindcast-F, we examine the impact of these three factors via case analyses of predictions for 1997/1998 and 2015/2016 events at a 7month lead time; that is, the hindcasts are initiated in May for the target season of December-January-February (DJF), as describe in this section.

Initial Conditions
Seasonal prediction can be seen as an initial value problem (e.g., Doblas-Reyes et al., 2003). To investigate the differences in the initial conditions of Hindcast-A and Hindcast-F, we compare the conditions produced by the AI and FFI (Figure 7). The observed tropical Pacific SSTAs during April 1997 and 2015 are positive over the equatorial and extratropical eastern Pacific, indicating that the warming signals emerge during boreal spring for these strong El Niño events. Both the AI and FFI reasonably reproduce the observed features, with comparable intensity and similar spatial patterns (Figure 7). The pattern correlation coefficients between the FFI (AI) and the observational reference are 0.77 (0.78) and 0.57 (0.59) for 2015 and 1997, respectively. Correspondingly, the root-mean-square errors (RMSEs) are 0.32 (0.31) and 0.34 (0.29). The results suggest the differences between the initial conditions derived from the FFI and AI are quite small Figure 6. Temporal evolutions of the observed and predictive Niño 3.4 index for El Niño events of (a) 2015/2016, (b) 1997/1998, and (c) 1982/ 1983. The curves are the Niño 3.4 index for the observations (black), Hindcast-A (red), and Hindcast-F (blue), respectively. Thin dashed and thick solid lines represent individual ensemble members and ensemble mean. and can be neglected. Thus, the accuracy of initial conditions is not the major cause of the differences in the predictive skill.

Model Drift in the Tropical Pacific
To examine model drift behavior, we calculate the drifts in the tropical Pacific based on equation (1) for the hindcasts initialized from May ( Figure 8). The model drifts of the seasonal mean SST and surface winds in Hindcast-A are largely reduced as expected (Figures 8a-8c). Conversely, Hindcast-F suffers from larger drifts from June-July-August to DJF with increasing integrations (Figures 8d-8f). The most prominent feature of drifts in Hindcast-F is an increasing warming of the equatorial central eastern Pacific, which develops from June-July-August and finally covers the whole tropical central eastern Pacific during DJF (Figures 8d-8f). A similar drift behavior for SST was found in the tropical Pacific when assimilating observational subsurface ocean temperature (Sanchez-Gomez et al., 2016), while some studies detect an opposite behavior in drifts to the model intrinsic state (Kim et al., 2017;Magnusson, Alonso-Balmaseda, Corti, et al., 2012); thus, the drift behavior may depend on the model and initialization scheme. As a response to the SST drifts over the central eastern Pacific, low-level winds converge to the equator from the subtropical regions, which are proportional to the intensity of the SST drifts and strengthen with increasing lead times (Figures 8d-8f). The results demonstrate that there are large differences between Hindcast-A and Hindcast-F in terms of model drift behavior. How does the drift affect the predicted signals? We investigate this issue in the next section.

ENSO Physical Processes
To understand the physical processes that cause the prediction differences, we compare El Niño evolution between the two hindcast types. During the developing phase (May to September), the El Niño growth rate in the Hindcast-A is about 2 times as strong as that in the Hindcast-F, so that the former El Niño intensity is closer to the observation in El Niño mature winter, while the latter is much weaker (Figures 6 and 9).
To further show the detailed physical processes, we performed the ocean mixed layer heat budget analyses in the equatorial central eastern Pacific (5°S to 5°N, 170-85°W) during the El Nino developing phase (May to

Journal of Advances in Modeling Earth Systems
September) for both 1997/1998 and 2015/2016 events ( Figure 10). The results indicate that the three typical positive feedbacks, Bjerknes SST-wind-thermocline feedback (−ω∂T ′ =∂z ), zonal advective feedback ( −u ′ ∂T=∂x ), and the Ekman pumping-induced anomalous upwelling feedback (−ω ′ ∂T=∂z), have no significant contribution to the underestimated warming tendency in the Hindcast-F. The differences in the warming tendency are mainly contributed by damping effects of the anomalous shortwave radiation (Q ′ SW =ρC p h, about 48%), and the advection of anomalous temperature by mean zonal currents ( −u∂T ′ =∂x , about 32%).
The Q ′ SW is closely linked to the overlying total cloud cover anomalies, with their spatial patterns highly consistent in both types of the hindcasts (Figures 11a-11d). For the Hindcast-A, the Q ′ SW exhibits a zonal dipole pattern, with increased clouds in the tropical central Pacific and reduced clouds in the tropical eastern Pacific (Figures 11a and 11b). In contrast, for the Hindcast-F, there are uniform positive total cloud anomalies and negative Q ′ SW over the tropical central eastern Pacific (Figures 11c and  11d). The difference in the damping effect of the Q ′ SW between the Hindcast-F and the Hindcast-A is caused by their difference in the total cloud fraction anomalies over the tropical central eastern Pacific (Figures 11e and 11f). Moreover, the anomalous total cloud fractions in the tropical central eastern Pacific are mainly contributed by the low clouds (Figures 11g and 11h).
The reduced low cloud over the tropical eastern Pacific in the Hindcast-A is closely associated with the equatorial cold tongue bias in FGOALS-s2 (Figure 12a) through following two processes. First, the excessive cold tongue prevents the total SST from exceeding the convection threshold in the tropical eastern Pacific even during El Niño year (Chen, Hua, et al., 2019;Graham & Barnett, 1987;Zhang, 1993). Second, the stronger climatological cold tongue leads to more stable marine boundary layer and thus more amount of marine stratiform ( Figure 12d). As a result, the cloud-SST relationship in the tropical central eastern Pacific is dominated by a positive feedback, instead of the negative feedback in the observation. The warm ENSO-related SSTAs reduce stability of the marine boundary layer and thus decrease marine stratiform cloud and increase downward shortwave radiative flux Philander et al., 1996).
In contrast, the equatorial central eastern Pacific in the Hindcast-F is dominated by negative feedback as in the observation. The warm SSTAs enhanced the deep convection and thus increase the total cloud and decrease the downward shortwave radiative flux (Figures 12c and 12d). This is because the El Niño-like SST drift in the Hindcast-F convert the air-sea interaction regime in the tropical central eastern Pacific from the stratiform-related positive feedback to deep convection-related negative feedback (Figures 12c and 12f).
The −u∂T ′ =∂x represents the restoring effect of the mean state caused by the thermal advection via zonal mean currents. It works as a negative feedback in the observation (Jin et al., 2006, Figure 3). The −u∂T ′ =∂ x term is very weak in the Hindcast-A, similar with that in the historical run (Figure 3). In contrast, the −u ∂T ′ =∂x is greatly intensified in the Hindcast-F, with magnitude close to that in the observation. The stronger −u∂T ′ =∂x in Hindcast-F is associated with the stronger equatorial climatological subsurface zonal oceanic currents ( Figure 13). Compared with the Hindcast-A, the El Niño-like SST drift in Hindcast-F (Figures 8d-8f) enhances the eastward climatological equatorial subsurface zonal currents in the central eastern Pacific through enhancing the geostrophic currents under the geostrophic balance (Hirst, 1986;   Wang & Weisberg, 1994), which causes stronger mean zonal advective negative feedback and thus weaker El Niño amplitude.

Journal of Advances in Modeling Earth Systems
In summary, the differences in the prediction are associated with the El Niño intensity growth rate during the period of May to September, when Hindcast-F underestimates the warming tendencies in the central eastern Pacific. The underestimation is primarily contributed by the terms of surface shortwave radiation (Q ′ SW ), the anomaly thermal advection by the mean zonal currents (−u∂T ′ =∂xÞ, associated with the model drifts in the SST and zonal ocean currents in the CEP, respectively.

Discussion
The above results in section 3.4.3 indicate that two damping effects by the Q ′ SW and −u∂T ′ =∂x in the Hindcast-F are greatly improved relative to the Hindcast-A and much closer to those in the observation. However, its predictive skill for the ENSO is significantly lower than the Hindcast-A. ENSO is a complex air-sea interaction mode, which contains a variety of positive and negative feedbacks. The relative magnitudes of these feedbacks in the model are quite different from those in the real world (Figure 3). Although the model drift due to the FFI improves some aspects of these feedbacks, it breaks their relative roles during the ENSO development and thus reduce the ENSO predictive skill.
In this study, the AI is superior to the FFI for the ENSO prediction. It is worth noting that this result should be model-dependent. Supposed that a model has a very weak bias in the tropical Pacific and can simulate all

Journal of Advances in Modeling Earth Systems
ENSO-related positive and negative feedbacks close to those in the observation, the FFI should be a better option, because it avoids the biases introduced when calculating the anomalies. First, the assimilated ocean temperature anomalies are calculated relative to the climatology derived from the EN4 gridded objective analysis data over the period of 1950-2005. It is only an expedient to use the gridded data as the climatological component of the unstructured ocean temperature profiles. Second, the tropical Pacific shows a La Nina-like warming trend during the last 40 years (Luo et al., 2012;McPhaden et al., 2011). The change in the background mean state was not considered in the calculation of the climatology, which would introduce biases to the initial conditions.

Conclusions
When initiated from observational conditions, coupled climate models used in seasonal predictions are typically subject to climate drifts. This study investigates the impact of model drift on ENSO prediction by using seasonal hindcasts conducted with two distinct initialization approaches, anomaly, and FFIs, in a climate prediction system termed IAP-DecPreS. It is constructed based on a state-of-the-art fully coupled global climate model, FGOALS-s2, and a newly developed EnOI-IAU ocean data assimilation scheme. Analysis shows that the ENSO predictive skill is dynamically influenced by model drift, whose impact can be partly overcome by using AIs. The underlying mechanisms are shown in Figure 14. The major conclusions are summarized in the following.
Measured by the metric of the anomaly correlation coefficient, the ENSO predictive skill in Hindcast-A is higher than that in Hindcast-F, particularly at long lead times. The 1997/1998 and 2015/2016 El Niño cases are used to investigate the physical processes causing the differences in predictive skill. When initialized from May, the drifts are largely reduced in Hindcast-A, while Hindcast-F features a growing warming of the equatorial central eastern Pacific. The major difference between the two hindcast types emerges during the period of May to September, during which the ENSO intensity growth rate in Hindcast-A is stronger than that in Hindcast-F associated with the greater warming tendency over the central eastern Pacific. Oceanic mixed-layer heat budget analyses are performed in this area. The results show that the underestimation of ENSO intensity in Hindcast-F is primarily contributed by the advection of anomalous temperature by mean zonal currents (−u∂T ′ =∂xÞ and damping effects of the anomalous shortwave radiation (Q ′ SW ), as both are modulated by model drifts.
The model drift in the hindcast based on the full field initialization influences ENSO prediction through both dynamic and thermodynamic processes. On the one hand, the drift enhances the eastward climatological equatorial subsurface zonal currents in the central eastern Pacific, leading to stronger mean zonal advective negative feedback and thus weaker El Niño amplitude. On the other hand, the drift shifts convert the air-sea interaction regime in the CEP from the stratiform-related positive feedback to deep convectionrelated negative feedback. Neither processes are favorable to ENSO development.
One important implication of this study is for the empirical bias correction technique in climate prediction. In previous studies and operational climate prediction, the model climate drift induced by initialization has been typically removed by statistical bias correction based on retrospective predictions Jin et al., 2008;Stockdale, 1997;Tippett et al., 2017). Although this technique works to some extent, our results show that the impact of model drift cannot be thoroughly removed via empirical bias correction, because model drift can interact with the predictive anomalies in a nonlinear manner. The nonlinear interaction of model drift with the prediction calls for special attention from both research and operational seasonal forecast communities.
More importantly, the key physical processes for the formations of the model drifts should be understood through more analysis, which can offer model developers some implications about how to reduce the model biases.