Mechanism of skillful seasonal surface chlorophyll prediction over the southern Pacific using a global earth system model

This study investigates the physical mechanism involved in an Earth system model (ESM)-based global marine biogeochemical prediction system providing successful forecasts of surface chlorophyll concentrations over the southern Pacific. The significant correlation skill of the surface chlorophyll concentration over the south-central Pacific (SP region, 160°–110° W, 10°–5° S) appears up to a 15-month lead. In contrast to the previously known role of the vertical nutrient supply on the predictive chlorophyll concentration forecasts, the NO3 budget analysis indicates that this prediction skill over the SP region is mostly controlled by the meridional advection of nutrients. Further analysis indicates that the controlling mechanisms involved in chlorophyll variability over the SP region can be explained by atmospheric and oceanic dynamics during the ENSO events. During La Nina, equatorial NO3 anomalies are increased due to enhanced equatorial upwelling, and the climatological southward current then advects nutrient-rich waters from the equator to the SP region i.e.,\,positive-v¯∂NO3′∂y\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( {{\text{i.e.,\,positive}}\, - \overline{v}\frac{{\partial {\text{NO}}_{{3}}^{\prime } }}{\partial y}} \right)$$\end{document}. In addition, anomalous easterly surface winds blow over the SP region as a circulation response to atmospheric diabatic heating anomalies during La Nina, which leads to southward current anomalies over the surface-layer ocean. This advects high climatological NO3 over the tropics to the subtropical south Pacific, which increases the NO3 anomalies i.e.,\, positive-v′∂NO3¯∂y\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( {{\text{i.e.,\, positive}} - v^{\prime}\frac{{\partial \overline{{{\text{NO}}_{3} }} }}{\partial y}} \right)$$\end{document}. This positive NO3 advection over the SP region is realistically simulated only at lead times shorter than 9 month, and the multi-season persistency of the nutrients contributes to the surface chlorophyll bloom at lead times longer than 1 year.


Introduction
Earth system models (ESMs) are the most complex among all interrelated system models, and they can be used to represent associated physical, chemical, and biological processes (Collins et al. 2011;Watanabe et al. 2011;Dufresne et al. 2013;Dunne et al. 2013;Lindsay et al. 2014). ESMs offer great potential to move beyond being physical descriptors of the Earth's system and become biogeochemical simulators of systems such as the global carbon cycle, primary productivity, nutrient cycling, and air quality (Anav et al. 2013;Stock et al. 2014b;Lin et al. 2015). With the advent of ESMs in recent decades, it is now understood that variability and future changes of physical climate system should be studied in an integrated process of physical, biological, chemical and socioeconomic components in Earth system.
Many studies have been conducted to investigate the interaction of marine biogeochemistry and physical climate system using ESMs together with satellite-retrieved ocean color data (Ballabrera-Poy et al. 2003, 2007Lin et al. 2008;Kang et al. 2017;Zhang 2015;Zhang et al. 2015Zhang et al. , 2018. Chlorophyll and related pigments in phytoplankton modify the vertical profile of shortwave absorption in the upper ocean, and this in turn affects upper ocean dynamics, largescale circulation, and climate variability. Tropical Pacific is the most active research area in terms of this bio-climate interaction study because of high primary productivity, strong biological shortwave heating, and high sensitivity to atmospheric feedback in this area, all of which can enable substantial impacts of biological feedback on climate variability, such as El Niño-Southern Oscillation (ENSO). A series of modelling studies have shown that biologicallyinduced climate feedback substantially changes the mean state of equatorial Pacific and consequently ENSO variability (Murtugudde et al. 2002;Lengaigne et al 2007;Loptien et al. 2009;Jochum et al. 2010;Lim et al. 2018), although different modeling studies showed different magnitudes of biological feedbacks potentially due to discrepancies in their experimental designs.
One of the most untapped potentials of ESMs is the ability to predict marine biogeochemistry. Oceanic biogeochemical cycles, including primary production, nutrient composition, carbon sequestration, oxygen concentration, and ocean acidity, are constantly changing due to the anthropogenic and natural variability of Earth's physical climate system. Thus, extending seasonal to decadal predictions of physical climate systems to those of marine biogeochemistry can enhance the ability to manage marine ecosystems and resources (such as controlling fishery harvests) (Brasseur et al. 2009;Tommasi et al. 2017a, b).
Given the passive behavior of biogeochemical variables to physical states, skillful marine biogeochemistry predictions require ESMs that accurately depict both physical and biogeochemical ocean states. However, while there is a clear need to predict marine biogeochemistry, developing biogeochemistry forecasting systems is challenging due to the model uncertainty and problems relating to biogeochemical initialization. In this respect, model uncertainty includes imperfections in the model structure and relates to parameters in both physical and biogeochemical models, whereas the biogeochemical initialization problem relates to the relatively sparse number of global biogeochemistry observations, the large number of biogeochemical tracers in ESMs, and the extremely high sensitivity of ocean biogeochemistry to dynamical imbalances induced during data assimilation, all of which hamper the development of fully coupled physical-biogeochemical initialization systems (Ford et al. 2012;Raghukumar et al. 2015;Song et al. 2016). These challenges inhibit extending current physical climate predictions to include those of global marine biogeochemistry predictions and integrating them within an operational climate data assimilation system. Therefore, no modeling group has yet reached the stage of enabling biogeochemical predictions using ESMs, and this currently remains in the developmental stage (Seferian et al. 2014;Ford and Barciela 2017;Rousseaux and Gregg 2017).
A recent study tackled the above-mentioned challenges and showed the feasibility of enabling ESM-based marine biogeochemistry predictions and its associated application (Park et al. 2019). The study used the extensive retrospective forecast experiments with an ESM and showed skillful marine ecosystem predictions in many ocean basins. It also showed the possibility of making skillful predictions up to 2 years in advance in some areas, as the ESM is capable of successfully simulating the evolution of marine biogeochemical responses to physical climate perturbations. The study further showed that ESM predictions can be utilized to anticipate observed fish catch fluctuations over the past 30 years in several coastal large marine ecosystems, suggesting that ESM-based predictions on seasonal to inter-annual time-scales will provide significant opportunities for society to adapt to future marine ecosystem changes.
Following the previous study conducted by Park et al. (2019), this study examines the detailed mechanisms involved in the predictive skill of marine biogeochemistry. While the previous work briefly suggested that the prediction skill arises from the longer memory of the marine biogeochemistry response to physical perturbations (such as ENSO), a detailed mechanistic understanding of the biogeochemical prediction skill on a regional basis has not yet been obtained. We focus here on the southern subtropical Pacific region, where the primary production level is intermediate compared to the equatorial Pacific to the north and the oligotrophic southern Pacific to the south. This region has been recently identified as an area where a longer-lead biogeochemical prediction skill is possible and the predictability of net primary production is substantially higher than that of SST (Seferian et al. 2014;Park et al. 2019;Taboada et al. 2019), while there has been little attempt to investigate what processes control the prediction skill. The high biogeochemical predictability and its underlying mechanism will be of great interest for the management of marine living resources, given the high number of tuna fisheries and aquaculture activities in the tropical Pacific (Lehodey et al. 1997).

Biogeochemical prediction system
The global marine biogeochemical prediction system used in the present study builds on the seasonal climate-prediction system developed at the Geophysical Fluid Dynamics Laboratory (GFDL). The prediction system is based on GFDL's ESM, which is a fully coupled atmosphere-land-sea iceocean model integrated with a marine ecosystem model, the Carbon, Ocean Biogeochemistry and Lower Trophics (COBALT) (Stock et al. 2014a, b). COBALT considers 33 tracers to resolve global-scale cycles of dissolved organic/ inorganic species with three explicit phytoplankton and three zooplankton groups. The horizontal resolution of the ESM is 2.5° longitude × 2° latitude for the atmosphere and land, and 1° × 1° for sea ice, the ocean, and marine biogeochemistry. The data assimilation system in this ESM-based biogeochemical prediction system is the ensemble coupled data assimilation (ECDA), which employs an ensemble Kalman filter assimilation scheme (Zhang et al. 2007). The ocean model in the ECDA is constrained by in situ ocean temperature and salinity observations, including oceanic profiles from the World Ocean Database (WOD), Argo profiles since 2000, and global temperature-salinity profile program (GTSPP) data sets since 2009. The daily SST field obtained from the National Oceanic and Atmospheric Administration's (NOAA) optimum interpolation SST v2 high resolution dataset (Reynolds et al. 2007) is also assimilated in the ECDA system. The atmosphere in the ECDA is constrained by the National Centers for Environmental Prediction, Department of Energy (NCEP-DOE) Reanalysis 2 (Kanamitsu et al. 2002), and 6-h mean atmospheric winds and temperature are assimilated in the ECDA. More details on the prediction and data assimilation systems are found in Chang et al. (2013) and Park et al. (2019).

Data
For the reference variations of the surface chlorophyll concentration, the satellite chlorophyll was obtained from the two ocean color sensors, the Sea-viewing Wide Fieldof-view Sensor (SeaWiFS) and the Moderate Resolution Imaging Spectroradiometer (MODIS) (Esaias et al. 1998;McClain et al. 1998). Further details about the model and the experimental design of the forecast run can be found in Park et al. (2019). The analyzed period is from 1998 to 2017.
For the reference variation of three-dimensional atmospheric and oceanic variables, the retrospective atmosphere/ ocean state estimates (such as 500 hPa vertical velocity and wind stress in the atmosphere, and nutrients, temperature, and currents in the ocean) produced from a historical run with the ECDA system are used. This run assimilates observed physical variables in the ocean and atmosphere but does not assimilate biogeochemical variables; therefore, the marine ecological model, COBALT, is integrated with the ECDA to drive ocean biogeochemistry in this physicallyconstrained ocean simulation. Thus, the simulation of biogeochemical variables benefits from improved physical states constrained by observations. A previous work has already shown that the global ocean biogeochemical reanalysis produced by this assimilation strategy is very capable of capturing observed large-scale biogeochemical patterns and variability (Park et al. 2018b). The data assimilation experiment was integrated for the period 1990-2017 and we analyzed the last 20 years of data to compare with observed satellite chlorophyll (which is available from 1998).
The forecast runs are initialized by data assimilative hindcasts from the ECDA system integrated with COBALT (which is the first dataset described above). These retrospective predictions are started at the 1st day of every calendar month from 1991 to 2017 with a set of 2 year-long, 12-member ensembles. In this study, we utilized the hindcasts starting from June 1, and note that the hindcast skill of any specific target month are not sensitive to the slight change of the selected forecast start month (see Fig. 1 of Park et al. 2019). Here we used ensemble mean predictions, and the predicted anomalies were calculated based on the lead-dependent monthly mean drift for each initialization month. The prediction skill was measured by the temporal anomaly correlation coefficient between the anomalies of predicted and satellite-retrieved chlorophyll concentrations (Park et al. 2019;Ham et al. 2019).

The climatology and variation of the surface chlorophyll concentration over SP region
Here, we first investigate the dominant interannual variability of surface chlorophyll over the southern subtropical Pacific in the satellite-based observation. The first empirical orthogonal function (EOF) of satellite-observed surface chlorophyll over the southern subtropical Pacific (180°-90° W, 5°-30° S) during boreal summer (June-July-August, JJA) exhibits a maximum amplitude at around 10°-20° S and has a tilted structure in a northwest-southeast direction (Fig. 1a). The center of the first EOF of the chlorophyll concentration anomaly is located over the southern edge of the equatorial chlorophyll concentration climatology (Fig. 1b), which implies the significance of the first EOF as an indicator of the southward extension of the equatorial chlorophyll-rich waters to the south Pacific. This tilted structure has been previously reported as being a southern branch of major tropical (30° S-30° N) chlorophyll variability, and its existence is apparent in the first EOF of summertime tropical surface chlorophyll (Behrenfeld et al. 2001;Yoder and Kennelly 2003). The previous studies also found that this leading mode of tropical chlorophyll is correlated with the ENSO with a time lag of 6-months (Park et al. 2011). This implies that the leading mode of the JJA-season surface chlorophyll over the southern Pacific would indicate chlorophyll variability during the decaying phase of ENSO. The processes how the ENSO leads the surface chlorophyll concentration variations over the equatorial Pacific is relatively well known; the enhanced upwelling due to the stronger trade wind during the La Nina increases the nutrient supply and the chlorophyll concentration over the surface layer (Chavez et al. 1999). On the other hand, the process how the ENSO-related equatorial chlorophyll variability is linked to the SP surface chlorophyll variation is not clear.
The index used in this study is defined by the areaaveraged surface chlorophyll concentration anomaly over the area 160°-110° W, 15°-10° S, where its leading mode has the maximum amplitude (denoted by the black box in Fig. 1). We refer to this area-averaged surface chlorophyll concentration anomaly as the SP chlorophyll index, and the target region is referred to as the SP region. It is of note that the chlorophyll forecast skill shown in this study is not sensitive to slight changes in the definition of SP region (not shown).
Among several nutrients, the surface chlorophyll variation over the SP region is mostly controlled by the variation of NO 3 . Figure 2 shows a scatter diagram between the monthly-mean nitrate (or iron) and surface chlorophyll anomalies over the SP region. For the comparison, same scatter diagram for the equatorial Pacific (160°-160° W, 5° S-5° N) is also shown. While the surface chlorophyll concentration over the equatorial Pacific is strongly coupled to the variation of the iron (Fig. 2b), those over the SP region is dominated by the variation of the NO 3 (Fig. 2c). This can be understood by the physiological characteristics of the phytoplankton; among several possible nutrients, the growth rate of the phytoplankton is mainly controlled by a most limiting one. That is, the phytoplankton is modulated by whatever nutrient is most limiting. As the climatological NO 3 over the SP region is much lower than that over the equatorial Pacific (i.e., the SP region is a NO 3 -limited region), the surface chlorophyll concentration is well correlated with the NO 3 .   (1), and e JJA(1) seasons. f Correlation skill of the area-averaged surface chlorophyll concentration index over the south Pacific region (i.e., 160°-110° W, 15°-10° S) with respect to forecast lead months (red line). The correlation skill of the persistent forecasts also shown (black line). Note that the most recent monthly-averaged value before the forecast start date (i.e., May-averaged value) is used as a persistent forecast result predicted in most of the tropical Pacific except for the east of the SP region. For example, based on the formula for the two-tailed t-test for the significance of the correlation coefficient (i.e., t = r √ n−2 √ 1−r 2 , where t is the P-value for the t-test, r and n are correlation skill and the number of samples, respectively), the correlation skill with a value above the 95% confidence level (i.e., 0.45) is apparent over most of the equatorial and southern subtropical Pacific, including the defined SP region (Fig. 3a). During September-October-November [SON(0)] season (i.e., averaged from 4-to 6-month lead time), although the skillful prediction skill area is reduced compared to that of the JJA(0) season, surface chlorophyll over the equatorial Pacific, subtropical northern, and SP regions are still predictable (green dotted area in Fig. 3b).

The prediction skill of the surface chlorophyll concentration over SP region
The surface chlorophyll concentration anomalies over the equatorial Pacific, southern subtropical, and northern off-equatorial Pacific are predictable up to December-January-February [D(0)JF(1)] season (i.e., 7-9 month lead average). However, for the forecast target season during the subsequent March-April-May [MAM(1)] season (i.e., 10-12 month lead average), surface chlorophyll concentration anomalies are only predictable over the south of the SP region. The correlation skill of the SST anomalies shows similar behavior that the equatorial SST is successfully predicted (i.e. correlation skill > 0.5) up to D(0) JF(1) season, however, the forecast skill rapidly drops during MAM(1) season (not shown). The abrupt drop of the forecast skill during boreal spring season is a common problem for most of model-based forecasts so called 'spring predictability barrier' (Zhang et al. 2005;Barnston et al. 2012). The area of the significant surface chlorophyll correlation skill expands again at JJA(1) season (i.e., 13-15 month lead average).
The correlation skill of the SP chlorophyll index exhibits a similar time-evolution (Fig. 3f); it is predictable (i.e., correlation skill > 95% confidence level) up to FMA(1) season, and reaches a minimum in the AMJ(1) season, which denotes that the correlation skill drop during boreal spring season is also evident not only over the equator but also over the SP region. Then, the correlation skill increases again and a second peak occurs in the JJA(1) season. This result is consistent with recent findings that the surface chlorophyll concentrations over the south subtropical Pacific is predictable more than 1-year lead time (Park et al. 2019). When the correlation skill of the climate model-based forecasts is compared to those of the persistent forecasts, it is clearly demonstrated that the model forecast has advantages to accurately predict the surface chlorophyll concentration over the SP region for lead months more than 1-year lead. On the other hand, the model-based forecast skill is similar or slightly lower than the persistent forecast skill for the lead month earlier than 9-month lead, which further emphasizes the necessity to focus on the model's surface chlorophyll predictability more than 1-year lead.
The high correlation skill of the surface chlorophyll concentration might be possible with aids of the multi-year persistent of the nutrients (Seferian et al. 2014), which is known to be predictable for longer lead time than SST. To examine the biogeochemical process behind the skillful forecast with a lead of more than 1 year, we analyzed the temporal evolution of subsurface biogeochemical anomalies related to the SP chlorophyll concentration in the JJA(1) season (i.e., 13-15 month lead average). Figure 4 shows the lag-regression coefficients between the observed SP chlorophyll index in the JJA(1) season and subsurface nitrate (NO 3 ) anomalies from the reanalysis (left panels) (or prediction (right panels)) data from the MJJ(0) to OND(1) seasons. Note that in contrast to the ironlimited equatorial Pacific the SP region is nitrate-limited, thus nitrate is considered to be a controlling nutrient in the present study as shown in Fig. 2. In the reanalysis, NO 3 anomalies in the upper-ocean down to a depth of 100 m gradually increase with the lead time, and a maximum regression coefficient is seen after a lead of 1-year (Fig. 4a). As the NO 3 signal in the surface layer is much larger than that in the subsurface layer below a depth of 50 m, and surface NO 3 is likely to be propagated to the subsurface from SON(0) to MAM(1) seasons, this implies that surface chlorophyll during the JJA(1) season may not simply be driven by the nutrient supply from nutrient-rich subsurface waters.
In the forecasts, the amplitude of the regression coefficient is largest at a depth above the 50 m during the D(0)JF(1) season (Fig. 4b). This indicates that the successful forecast of SP surface chlorophyll variations in the JJA(1) season is closely linked with mixed-layer NO 3 variations during the previous lead months. The detailed timing of the predicted NO 3 anomalies coupled to the SP chlorophyll index in JJA(1) season is different from that of the reanalysis, and the vertical propagation of NO 3 anomalies from the subsurface layer is also weak, as seen in the reanalysis.
It is unclear how SP surface chlorophyll concentration anomalies are well predicted in the model even though it is controlled by nutrient anomalies at different timings from those of the reanalysis. In next section, we will examine the detailed process of the time evolution of the SP chlorophyll index and the NO 3 anomalies through the NO 3 budget analysis.

Physical mechanism involved in the SP
surface chlorophyll variations Figure 5 shows the composited mixed-layer (0-45 m average) NO 3 anomalies and NO 3 budget terms over the SP region during the years when the satellite-observed SP chlorophyll index in the JJA(1) season was greater than 0.5 standard deviation. The composited cases represent forecasts initialized on June 1st of 1998, 1999, 2000, 2008, 2010, 2011, and 2012. The composited SP chlorophyll index from the satellite-observation and the forecast is also shown. The time evolution of NO 3 anomalies from the reanalysis indicates that the positive NO 3 anomalies in the JJA(1) season are mostly contributed to by the increase of the NO 3 anomalies in the JJASON(0) and AMJJA(1) seasons (Fig. 5a). The time-variation of the satellite-observed SP chlorophyll index exhibits a continuous increase in time, denoting that the chlorophyll variation is well explained by those of the NO 3 anomalies.
The forecast well mimics the increase in NO 3 anomalies during the JJASON(0) season, however, the second increase during the AMJJA(1) season is poorly simulated. As a result, the predicted SP surface chlorophyll anomalies rapidly increase during the JJASON(0) season and then continuously decay since then. Nevertheless, the predicted SP chlorophyll index is still positive up to JJA(1) season due to the strong increase from JJA(0) to OND(0) season, denoting that the predicted SP chlorophyll concentration is also higher than normal when the satellite-observed SP chlorophyll concentration is higher than normal. This means that the interannual variation of the SP chlorophyll index is successfully simulated in the forecast, which is consistent with the predictive correlation skill over 1-year lead. This implies that successful forecasts in the model mainly originate from a realistic simulation of the time evolution of NO 3 anomalies during the JJASON(0) season. As the second increase in the NO 3 anomalies during the AMJJA(1) season is not simulated in the model, the SP surface chlorophyll anomalies in the JJA(1) season are strongly coupled to the NO 3 anomalies during previous seasons as shown in Fig. 4b.
Prior to analyzing the mechanism involved in the time variation of the NO 3 anomaly using a budget analysis, we firstly determined whether the time tendency of the NO 3 anomaly can be well explained by the sum of linearized NO 3 advections. The time variation of the total NO 3 advection (black line in Fig. 5b) is well matched with that of the time tendency of the NO 3 anomalies (green line in Fig. 5b); the time tendency of the NO 3 anomalies exhibits a relatively large value during the JJASON(0) and AMJJA(1) seasons, and it is well explained by total NO 3 advection. Also, the  (Table 1). Therefore, an analysis of the NO 3 budget can be used to diagnose why a strong increase in the NO 3 anomaly occurs during the JJASON(0) and AMJJA(1) seasons than other seasons. Note that total NO 3 advection is systematically overestimated compared to the NO 3 time tendency as the biological sink of the NO 3 (i.e., the decrease in the NO3 anomalies due to the consumption by the chlorophyll) is not considered.
The budget analysis clearly shows that the meridional advection of NO 3 is responsible for the increased NO 3 anomalies over the SP region in both the reanalysis and forecasts (Fig. 5b). In general, the positive advection of anomalous NO 3 due to the climatological meridional current  Table 1 The square of the correlation coefficient (r 2 ) between the sum of linearized NO 3 advections and the time-tendency of the NO 3 anomalies (2nd column), and each linearized advection term (from 3nd to 8th columns) averaged over the SP region in the reanalysis during 1998-2017 the positive NO 3 anomalies in the forecast (Fig. 5c). This indicates that the horizontal advection of NO 3 is responsible for the successful forecast of surface chlorophyll anomalies over the SP region. It should be noted that the vertical advection terms, which are reported to be a key process in the modulation of surface nutrients over the equator and coastal regions (Messié and Radenac 2006;Park et al. 2018a), are not important in the time evolution of NO 3 anomalies over the SP region. The role of meridional NO 3 advection due to the climatological currents is robust not only for the composited year, but for all-year data (Table 1). Among all linearized NO 3 advection term, the meridional advection of the NO 3 due to the climatological current explains 85% of total NO 3 advections, supporting the critical role of the meridional transport of the NO 3 anomalies on the surface chlorophyll variability over the SP region. While the meridional NO 3 advections are well simulated in forecasts during the early forecast lead seasons, there are systematic errors in the NO 3 forecasts after the D(0)JF(1) season: the amplitude of the simulated meridional advection terms is underestimated, and the forecasts shows the excessive negative zonal NO 3 advection i.e., −u NO � 3 x from the D(0)JF(1) to MJJ(1) season (cf. Fig. 5b, c). Due to the weak positive meridional advection and excessive negative zonal advection, total NO 3 advection exhibits negative values in the forecasts of the JFM(1) to MJJ(1) seasons. This contributes to lower mixed-layer NO 3 and the surface SP chlorophyll index in the forecasts compared to the reanalysis. Nevertheless, surface chlorophyll anomalies over the SP region is still positive for the forecast longer than 1-year lead, which denotes that a realistic simulation of the meridional advection of the NO 3 during the relatively early stage of the forecast is responsible for the successful forecast of surface chlorophyll anomalies over the SP region at a lead time of more than 1-year.
The composited NO 3 for the negative SP chlorophyll case (i.e., when the satellite-observed SP chlorophyll index in the JJA(1) season was smaller than negative 0.5 standard deviation) generally shows a mirrored time-evolution of the NO 3 (Fig. 6), however, there is also few systematic differences. In the forecast, the positive −v NO � 3 y is nearly zero for early lead time [i.e., from MJJ(0) to ASO(0) seasons]. This leads less accumulation of the NO 3 anomalies compared to the reanalysis, as a result, the amplitude of the forecasted SP chlorophyll index is systematically smaller than the satelliteobservation (Fig. 6a). For example, the satellite-observed SP chlorophyll index shows − 0.4 μmol∕kg during D(0)JF(1) season, while the forecasted value is nearly zero. As the NO 3 accumulation is quite insufficient, the NO 3 advection at forecast lead time longer becomes crucial for the variation of the SP chlorophyll index during JJA(1) season.
During the MAM(1) season, the wrongly simulated positive −u NO � 3 x exhibits peak value (Fig. 6c), which leads a spurious positive surface chlorophyll anomalies over the SP region (red dotted line in Fig. 6a). As the simulated sign of the SP chlorophyll index is opposite to the satellite-observation during MAM(1) and AMJ(1) seasons, this contributes to drop the correlation skill during MAM(1) and AMJ(1) season. On the other hand, the amplitude of the arbitrary positive −u NO � 3 x becomes weak during JJA(1) season, the forecasted SP chlorophyll index exhibits negative value as shown in the satellite-observation. This is consistent with the recovery of the correlation skill during JJA(1) season, which will be discussed in Sect. 6 in more detail.
To understand the process for positive SP chlorophyll index case during JJA(1) season, the meridional NO 3 advection −v NO � 3 y for the positive SP chlorophyll cases (Fig. 5a) is further decomposed into the climatological current and NO 3 anomalies. Figure 7 shows the climatological mixedlayer (0-45 m average) meridional current during 1998-2017 and the composite of the anomalous mixed-layer NO 3 anomaly averaged during the JJASON(0) season for the positive SP chlorophyll index cases (i.e., in the forecast cases where the satellite-observed SP chlorophyll index in the JJA(1) season was greater than 0.5 standard deviation).
In the reanalysis, the NO 3 anomalies have a maximum value over the equatorial eastern Pacific, which leads to a positive meridional gradient for the anomalous NO 3 NO ′ 3 y > 0 over south of the equator. Therefore, the southward climatological current ( v < 0) induces positive meridional NO 3 advection −v NO � 3 y > 0 in the reanalysis. As the anomalous NO 3 in the forecast has a maximum value over the northern edge of the SP region, the southward climatological current also leads to positive meridional NO 3 advection. It is interesting to note that the positive NO 3 anomalies over the equatorial eastern Pacific are not well simulated in the forecasts; this is consistent with the shorter predictable lead time of surface chlorophyll concentration anomalies over the equatorial Pacific compared to those over the SP region (Fig. 3d, e).
The composite of the anomalous mixed-layer meridional current v ′ and climatological mixed-layer NO 3 ( NO 3 ) during the JJASON(0) season are also well captured in the forecasts (Fig. 8b, d). The anomalous current exhibits negative values, and the meridional gradient of climatological NO 3 is positive in the SP region both in the reanalysis and the forecasts. This leads to positive meridional NO 3 advection by the anomalous meridional current −v � NO 3 y > 0 over the SP region.
While the spatial distribution of mixed-layer NO 3 and the current anomalies are relatively well simulated during the JJASON(0) season, those during the JJA(1) season are problematic; the amplitude of the simulated NO 3 anomalies over the SP region is only about 30% of the that in the reanalysis (Fig. 9). This is consistent with the weak surface chlorophyll anomalies over the SP region in the forecasts than that in the reanalysis during JJA(1) season (Fig. 5a). This is also It is shown that the spatial distribution of the anomalous mixed-layer NO 3 and the currents in early forecast lead seasons are well simulated in the forecasts, which gives the NO 3 anomalies a realistic time evolution. The realistic nutrient supply at early forecast lead can lead to the successful surface chlorophyll prediction over the SP region more than 1-year lead. This might be possible due to the multi-season persistency of the nutrients, which is much prominent than those of physical variables (Seferian et al. 2014). The cause of these realistic spatial distributions of anomalous currents and NO 3 thus needs to be determined, therefore, the following section examines the large-scale features inducing NO 3 and the current anomalies associated with the surface chlorophyll anomalies over the SP region.  Figure 10 shows the anomalous sea surface temperature (SST) and atmospheric vertical pressure velocity at 500 hPa in the reanalysis and forecasts in the positive SP chlorophyll index cases. In JJA(0) season, negative SST anomalies are clearly observed over the equatorial Pacific. The amplitude of equatorial negative SST anomalies in the reanalysis is intensified up to the D(0)JF(1) season, which represents the La Nina event in its mature phase during the boreal winter season (Tziperman et al. 1998;Ham and Kug 2014). The overall amplitude of the negative SST anomalies over the equatorial eastern Pacific is then weakened in MAM(1) season; in particular, the equatorial SST anomalies to the east of 140° W are almost diminished. In contrast, the negative SST anomalies over the equatorial central Pacific and the subtropical centraleastern Pacific around 180°-120° W are maintained until MAM(1) season. The tropical cold SST anomalies (i.e., La Nina event) indicate the occurrence of enhanced equatorial oceanic upwelling over the equatorial central-eastern Pacific (Jin 1997), and the increased equatorial mixed-layer nutrients as a result (Gorgues et al. 2010). It is also evident in the composite analysis that the increased NO 3 anomalies are robust over the equatorial Pacific (Fig. 7c) is linked to the local cold SST anomalies (Fig. 10e). As previously mentioned, the equatorially increased NO 3 anomalies are During the La Nina event, the negative SST anomalies over the equatorial central-eastern Pacific suppress atmospheric upward motions (contours in Fig. 10) and reduce convective activities; this results in negative precipitation anomalies over the equatorial central Pacific (Philander 1985;Kang and Kug 2002;Okumura and Deser 2010;Cai et al. 2015). As the amplitude of the negative SST anomalies are strongest during the D(0)JF(1) season, the downward vertical velocity anomaly in the atmosphere over the equatorial central Pacific also has a maximum amplitude during this season.
In the forecasts, the negative SST anomalies and the associated downward motion anomalies are well simulated in the JJA(0) season. However, the forecast is problematic that the negative SST anomalies are dampened faster than the reanalysis. For example, there is a continuous weakening of the amplitude of negative SST anomalies with an increased forecast lead time, and the amplitude is reduced to approximately 20% that of the reanalysis in the D(0) JF(1) season. This is more evident in the MAM(1) season, where that the negative SST anomalies are barely shown over the entire tropical Pacific. As the SST anomalies in the forecast are systematically weaker than those of the observed anomalies, the amplitude of the atmospheric downward vertical velocity over the equatorial central Pacific is also weaker than that in the reanalysis.
The fast decay of the La Nina SST signal in the forecast might be related to unrealistic El Nino-Southern Oscillation (ENSO) phase-locking in the GFDL models (Ham et al. 2013;Ham and Kug 2014). Although the observed ENSO SST is in a mature phase during the boreal winter season, the GFDL-ESM2M simulates a summer ENSO peak because the air-sea coupling strength is too weak during the boreal summer season (Ham and Kug 2014). As a result, the observed La Nina signal continuously grows until the D(0)JF(1) season, while the negative SST anomalies related to the La Nina event cannot be intensified.
Even though the overall amplitudes of SST and atmospheric circulation anomalies are weaker in the forecast, the negative SST anomalies over the subtropical south Pacific are relatively well maintained up to the D(0)JF(1) season. In addition, atmospheric downward motion anomalies appear around the dateline of the subtropical southern Pacific in the forecast (Fig. 10d, f). This is different from the anomalies over the equatorial eastern Pacific SST and vertical atmospheric velocity anomalies, which completely disappear at this lead time. This (i.e., longer preservation of the negative SST anomalies over the SP than the equator) can be possible due to the local air-sea coupled process , or remote atmospheric teleconnection excited by the wave source over the equatorial western Pacific (Folland et al. 2002;Rodrigues et al. 2015), however, the detailed examination of the possible mechanisms is out of a main research scope of this study.
When the anomalous atmospheric downward motion becomes dominant over the equatorial central Pacific, it leads to low-level atmospheric divergence: easterly anomalies over the western-central Pacific and westerly anomalies over the eastern Pacific (Fig. 11a, anomalies for JJASON(0) season). Simultaneously, the lower-level divergence leads to low-level equatorward atmospheric flow through the Sverdrup balance: northerly to the north of the equator and southerly to the south of the equator (Clarke et al. 2007). As a result, cyclonic (i.e., clockwise) flow is induced over the subtropical south-eastern Pacific.
The surface wind response over the equator in the forecasts is fairly weak, and it is shifted to the west (Fig. 10b); this is consistent with the significantly weak negative SST anomalies over the equatorial Pacific during the JJASON(0) season. However, the surface wind responses over the subtropical south Pacific are relatively well simulated with stronger amplitudes. In particular, easterly wind stress anomalies with amplitudes of 0.003 N/m 2 are located along 15° S, and these are also shown in the reanalysis. This indicates that in both the reanalysis and forecasts, the SP region is dominated by easterly wind stress anomalies, which in turn generate a southward current over the mixed-layer through Ekman transport (Wickett 1967;Price et al. 1987). It is of note that the spatial distribution of the southward current ( Fig. 8a, b) is quite consistent with that of the easterly wind stress anomalies (Fig. 11a, b). The southward current advects high NO 3 from the equatorial regions; therefore, it is responsible for the positive meridional NO 3 advections over the SP region i.e., positive − v � NO 3 y .

Recovery of surface chlorophyll forecast skill
This section investigates the recovery of the SP chlorophyll index forecast skill: the correlation skill is minimum during AMJ(1) season, then it is recovered to some extent during JJA(1) season (Fig. 3f). As it is appeared that the recovery of the correlation skill is much clear for the negative SP chlorophyll index composite (Fig. 6, with selected cases when the satellite-observed SP chlorophyll index during JJA(1) season is lower than − 0.5 standard deviation) than the positive case composite, we firstly examine the negative SP chlorophyll case. In Sect. 4, we argued that the following two points are responsible for the minima of the correlation skill during AMJ(1) season for negative SP chlorophyll index cases. (1) In the forecasts, −v NO � 3 y is too weak during JJASON(0) season, this leads insufficient accumulation of the NO 3 anomalies over the SP region at early forecast lead time. As a result, the simultaneous NO 3 advection becomes crucial for the variation of the SP chlorophyll index at longer than 1-year lead. (2) In the forecast, −u NO � 3 x is too strong during AMJ(1) season, which leads opposite sign of the forecasted SP chlorophyll index to that of the satellite-observation. The amplitude of the error in −u NO � 3 x is weakened during JJA(1) season, the forecast skill of the NO 3 and SP chlorophyll anomalies is recovered to some extent.
To examine point (1), Fig. 12 shows the negative-case composited SST and NO 3 anomalies during JJASON(0) season. The El Nino-related SST anomalies over the equatorial Pacific is clear in the reanalysis. The reduced upwelling, which is closely linked to the positive SST anomalies, is responsible for the less-than-normal NO 3 anomalies over the equatorial Pacific. On the other hand, the positive SST anomalies are barely remained over the equatorial Pacific except for the off-equatorial western Pacific in the forecast. As a result, the negative NO 3 is only shown over the offequatorial western Pacific.
The fast damping of the El Nino SST anomalies in the forecast is a mirror image of the fast decaying of the La Nina SST anomalies shown in the Fig. 10, denoting that both the forecasted El Nino and La Nina SST damps systematically faster than the reanalysis. On the other hand, the degree of excessive SST damping is much severe during the El Nino in the negative SP chlorophyll composite than during the La Nina in the positive SP chlorophyll composite. As a result, the amplitude of the forecasted NO 3 anomalies over the north of SP region during the El Nino (Fig. 12d) is only half of those during the La Nina (Fig. 7d) during JJA-SON(0) season. The leads weaker meridional gradient of the NO 3 anomalies, and weakens the meridional NO 3 advection (a) (b) (c) (d) Fig. 12 The SST (upper panels) and mixed-layer (0-45 m average) NO 3 (lower panels) anomalies during JJASON(0) season in reanalysis (left panels) and forecasts (right panels) for the negative SP chlorophyll cases at early forecast lead month in the negative SP chlorophyll cases.
To investigate point (2), Fig. 13 shows the NO 3 anomalies during AMJ(1) season for the negative SP chlorophyll cases. In the reanalysis, the negative NO 3 anomalies, whose amplitude was maxima over the equatorial Pacific during previous seasons, shifts to the south during AMJ(1) season. The amplitude of the NO 3 anomalies over the SP region exhibits the pronounced difference in the meridional direction, while the difference is quite small in the zonal direction. As a result, the meridional gradient of the NO 3 anomalies is robust than the zonal gradient in the reanalysis. On the other hand, the forecasted NO 3 anomalies during AMJ(1) season shows the robust zonal gradient over the SP region; the negative anomalies over the western part of the SP region and the positive anomalies over the eastern part. The excessive positive NO 3 anomalies centered at 110° W, 20° S might be related to a clockwise surface circulation and the resultant upwelling excited by the remote atmospheric forcing (Rodrigues et al. 2015) (not shown). Due to the strong positive zonal gradient of the NO 3 anomalies and the stronger climatological easterly in the forecasts than that in the reanalysis (Fig. 14a, b), spurious positive NO 3 advections can be induced in the forecasts, which is never shown in the reanalysis.
The amplitude of the positive zonal advection of the NO 3 anomalies is weakened to some extent during JJA(1) season. It is appeared that the different amplitude in the climatological zonal current is responsible for this amplitude change in the zonal advection (Fig. 14). That is, the simulated climatological easterly currents during JJA season is weaker than that during AMJ season, and this leads weaker zonal NO 3 advection during JJA(1) season in the forecasts. The spatial distribution of the NO 3 anomalies is not much changed from AMJ(1) to JJA(1) season (not shown). Then, one can ask whether the recovery of the correlation skill is totally due to the events selected for negative SP chlorophyll cases. It is found that the some of the positive SP chlorophyll cases also largely contributes to the recovery of the correlation skill. Figure 15a shows a normalized difference in the root-mean square error (RMSE) between the reanalysis and the predicted SP chlorophyll index in the AMJ(1) season compared to that in the JJA(1) season. To enable a good comparison of the RMSE differences to be made between target years with different initial forecast errors, the RMSE difference is normalized by the RMSE in the JJA(1) season. It is clear that degradation of the correlation skill in the AMJ(1) season is most pronounced for forecasts initialized on June 1, 1998. In other forecast cases (i.e., forecasts initialized on June 1 in 1999, 2000, 2008, 2010, 2011, and 2012), the difference in the forecast skill between the two seasons is relatively small. We thus focused on the forecast initialized on June 1, 1998 to understand the recovery of the forecast skill in the JJA(1) season.
The time evolution of NO 3 anomalies and a SP chlorophyll index for the forecasts initialized on June 1, 1998 shows quite similar feature emphasized for the negative SP chlorophyll composite (Fig. 15b). That is, the NO 3 anomalies are negative during early forecast lead month (roughly up to 6-months), which indicates the insufficient NO 3 accumulation. The forecasted NO 3 exhibits negative values during MAM(1) and AMJ(1) seasons, while the reanalysis NO 3 exhibits the positive values during the whole time. The excessive negative zonal NO 3 advection is responsible for the negative NO 3 anomalies during those seasons. During JJA(1) season, the forecasted SP chlorophyll anomalies becomes positive as shown in the reanalysis. This supports our previous arguments that the excessive zonal NO 3 advection during MAM(1) and AMJ(1) season is responsible for the low correlation skill during those lead times.

Summary and discussion
This study examines the physical mechanism involved in using an ESM to make successful surface chlorophyll concentration forecasts over the south-central Pacific with a lead time of up to more than 1-year. Interannual variability of the surface chlorophyll concentration during boreal (a) (b) Fig. 15 a Difference in root-mean square error (RMSE) between satellite-observed and predicted SP chlorophyll index in AMJ(1) season and JJA(1) season divided by RMSE in JJA(1) season. Note that only target years used in the positive SP chlorophyll composite are shown. b Time-evolution of reanalyzed (black solid line with circle) and predicted (red solid line with circle) mixed-layer (0-45 m aver-age) NO 3 anomalies over the SP region for 1998 forecast case. Timeevolution of satellite-observed (black dotted line) and predicted (red dotted line) surface chlorophyll concentrations over the SP region is also shown. The zonal NO 3 advection term due to the climatological zonal currents in the reanalysis (black solid line without circle), and the forecasts (red solid line without circle) are also shown summer is the strongest over the south-central Pacific (i.e., the SP region) and the significant correlation skill for the surface chlorophyll concentration appears up to JJA(1) season. The NO 3 budget analysis indicates that nutrients required to accumulate the surface chlorophyll concentration over the SP region are mostly provided by the meridional advections.
The composite of SST anomalies shows that cases with higher-than-normal surface chlorophyll concentrations over the SP region occur during La Nina events. The anomalous easterly surface wind blows over the southcentral Pacific as a circulation response to atmospheric diabatic heating anomalies during La Nina, which leads to southward current anomalies over the mixed-layer ocean. This also advects high climatological NO 3 over the tropics to the subtropical south Pacific, which increases the NO 3 anomalies (i.e., positive −v � NO 3 y ) (Fig. 16a). In addition, the equatorial NO 3 anomalies are increased due to the enhanced equatorial upwelling during La Nina years, and the climatological southward current then advects nutrient-rich waters from the equator to the SP region i.e., positive − v NO � 3 y (Fig. 16b). This positive NO 3 advection is strongest during the JJA-SON(0) and AMJJA(1) seasons in the reanalysis. However, only the early increase during the JJASON(0) season is well simulated in the model. In addition, negative zonal NO 3 advection over the SP region is excessively simulated in the forecasts from the MAM(1) and AMJ(1) seasons, which prevents the increase of the surface chlorophyll concentration over the SP region. This means that realistic simulations of NO 3 anomalies at early forecast lead months shorter than 6-months are responsible for successful forecasts of the surface chlorophyll anomalies over the SP region at lead times of up to more than 1-year.
If the NO 3 anomalies are insufficiently accumulated during early forecast lead months, as in the negative SP chlorophyll composite cases or 1998 case, the NO 3 anomalies at lead times of longer than 1-year lead are mainly determined by the simultaneous NO 3 advection. In negative SP chlorophyll cases (i.e., cases when SP chlorophyll index during JJA(1) season is smaller than − 0.5 standard deviation), the positive zonal NO 3 advections in the MAM(1) and AMJ(1) season can contribute to degrade the correlation skill of the surface chlorophyll concentration over the SP region. As the positive zonal NO 3 advection during JJA(1) season is not excessive much as that during AMJ(1) season, the correlation skill recovered during the JJA(1) season to exhibit a second peak.
The supply of nutrient-rich subsurface water from below the thermocline is considered to be key to surface Fig. 16 The schematic diagram of successfully forecasts of the chlorophyll blooming over the SP region during JJA(1) season. The positive meridional advection of the NO 3 during JJASON(0) as in a, b leads SP chlorophyll blooming in the forecasts, then, the positive chlorophyll and NO 3 anomalies are persisted for multi-seasons to successfully mimics the chlorophyll blooming during JJA(1) season in the reanalysis chlorophyll blooming over many areas. To explain the variability in concentrations of surface chlorophyll, many studies have been conducted to diagnose the sea-level anomalies (Murtugudde et al. 1999;Wilson and Adamec 2001), or wind speed is an indicator of vertical mixing (Chiswell et al. 2013), to explain the variability of the surface chlorophyll concentration. In contrast, this study demonstrates that the horizontal advection of nutrients mainly modulates the surface chlorophyll concentration over the SP region.
Not only over the SP region, it is also possible that the similar mechanism can work for other regions where the horizontal gradient of the nutrient is robust. In particular, the subtropical north Pacific exhibits similar degree of the correlation skill to those over the SP region with similar characteristics of the horizontal nutrient distributions (i.e., the equatorial maxima of the nutrient and its abrupt decrease to the pole), and the budget analysis for the subtropical north Pacific region showed the similar mechanism would be responsible for the interannual variation of the surface chlorophyll variations over the subtropical north Pacific (not shown). Therefore, conducting a three-dimensional budget analysis is beneficial to be utilized for the future studies for gaining a full understanding of the regional nutrient supply.