Skill Assessment of Seasonal-to-Interannual Prediction of Sea Level Anomaly in the North Pacific Based on the SINTEX-F Climate Model

Extreme sea level rise seriously impacts habitation and is indicative of changes in primary production in the North Pacific. Because of its rising trend associated with global warming, skillful seasonal-to-interannual predictions have become increasingly valuable to guide the introduction of suitable adaptation measures that help us reduce the risks of socioeconomic losses. Here, we have used a dynamical coupled ocean–atmosphere model called “SINTEX-F” to revisit the potential predictability of sea level anomalies at a lead of up to about 2 years. Skillful prediction is found mainly in the tropical Pacific as shown in previous work. Here, we found a new skillful prediction region in the North Pacific (30°–40°N, 180°–150°W) at about 2 years’ lead time. We also analyzed the co-variability among ensemble members and found the long-lasting ENSO/ENSO-Modoki in the tropical Pacific seems to contribute to the predictability source. The result may be useful to develop systematic and synergistic attempts to predict marine ecosystem responses to regional and global climate variations.


INTRODUCTION
The North Pacific marine ecosystems are primary sources of ecosystem services (e.g., fishing, shipping, and recreation) for its surrounding countries including Canada, U.S., China, Russia, Korea, and Japan. Coastlines of the North Pacific are seriously damaged by extreme sea level rise (Nicholls et al., 2007). In particular, the coastal zones are immediately affected by submergence and increased flooding of coastal land, as well as saltwater intrusion of surface waters (Nicholls and Cazenave, 2010).
In addition to the rising sea level trend associated with the global warming, extreme sea level events occur in association with natural climate variability such as the Pacific Decadal Oscillation (PDO), the Interdecadal Pacific Oscillation (IPO), the North Pacific Gyre Oscillation (NPGO), and the El Niño-Southern Oscillation (ENSO) in the North Pacific (Mantua et al., 1997;Zhang et al., 1997;McGowan et al., 1998;Lombard et al., 2005;Di Lorenzo et al., 2008;Hamlington et al., 2019;Han et al., 2019). To address relatively short-term risks, stakeholders desire a forecast of monthly/seasonal rising or falling sea levels caused by those climate variabilities. Hinkel et al. (2019) analyzed user needs for sea level rise information, and how they are able to be met given the state-of-the-art of sea level forecast science. Jacox et al. (2020) reviewed statistical and dynamical marine ecosystem forecasting methods and highlighted examples of their application along U.S. coastlines for seasonal-to-interannual prediction. Payne et al. (2017) also reviewed the state of the art marine ecological forecasts and suggested forecasts ranging from seasonal to decadal time scales are now a reality. Tommasi et al. (2017a) evaluated the multiannual SST predictions over Large Marine Ecosystems (LMEs), a coastal scale relevant to managed fisheries stocks. Tommasi et al. (2017b) also highlighted advances in seasonal to decadal prediction of managing living marine resources in a dynamic environment. Those previous studies provide information relevant for supporting coastal adaptation decision making.
Although skillful predictions of SST have already proven useful for a number of marine resource applications (e.g., Hobday et al., 2014;Stock et al., 2015), further studies about sea level anomaly are necessary. Rebert et al. (1985) showed that the oceanic Kelvin and Rossby waves have a direct relation between thermocline depth and sea level, while they have only an indirect relation to SST. These ocean dynamics are responsible for the relatively high skill of sea level prediction relative to SST prediction (Miles et al., 2014). Zainuddin et al. (2017) found that SST was an important variable for detecting hotspots of skipjack tuna distribution, as they are sensitive to the changes in temperature. Sea level anomaly is related to the changes in the depth of the thermocline and mesoscale variability. They combined these variables to improve detection of potential pelagic hotspots for skipjack tuna.
To expand prediction of large-scale sea level anomalies into coastal areas and to further the understanding of its potential predictability, it is necessary to evaluate the lead-time and locations in which a dynamical, physics-based prediction system performs well. It might allow coastal communities to better adapt for the impacts of severe flooding and erosion driven by high sea levels.
Although decadal climate variation is more predictable than previously thought, it is still challenging (Meehl et al., 2014;Smith et al., 2019). Here, our focus is on seasonal-to-interannual prediction. Generally speaking, the most important potential source of seasonal-to-interannual predictability is often from ENSO events, which develop via air-sea coupled feedback. Therefore, application of an ocean-atmosphere coupled general circulation model (GCM) is naturally a proven approach to overcome shortcomings of stand-alone atmospheric/oceanic models. Miles et al. (2014) initially attempted to apply a coupled GCM to predict seasonal sea level anomalies, and assessed the skill globally for up to 7 months in advance. McIntosh et al. (2015) showed the prediction skill by dynamical GCMs is better relative to statistical approaches for coastal sea level. Polkova et al. (2015) used the decadal prediction system and found predictability in the subtropics. Roberts et al. (2016) assessed the predictability of large-scale dynamic sea level anomalies up to 15 months using a climate model and found that prediction of seasonalto-interannual sea level variability in the extratropics will be governed by the predictability of surface wind stress and modes of atmospheric variability.
This study is a follow-up study of those pioneering studies. Here, we have revisited the predictability of sea level up to about 2 years in advance by analyzing results of a coupled ocean-atmosphere general circulation model "SINTEX-F." Such a long lead time retrospective forecast is beyond most current operational capabilities, and hence a skill assessment of the model results is conducted here as a first attempt. We believe that the obtained result is useful to attempt systematic and synergistic prediction of marine ecosystem responses to regional and global climate variations.

Dynamical Prediction System
The Scale Interaction Experiment-Frontier ver. 1 (SINTEX-F1) prediction system was used here, which is based on a fully coupled global ocean-atmosphere circulation model (CGCM) developed under the EU-Japan collaborative framework (Luo et al., 2003;Luo et al., 2005;Masson et al., 2005). The atmospheric component has a horizontal resolution of 1.125 • (T106) with 19 vertical levels. The oceanic component has a horizontal resolution of about 2 • × 2 • but with meridional refinement to 0.58 • in the tropics. It has 31 vertical levels from the surface to the bottom with a relatively finer resolution of 10 m from the sea surface to 110 m depth. This system adopts a relatively simple initialization scheme based only on the nudging of observed SST. In consideration of uncertainties of both initial conditions and model physics, it has nine ensemble members. More details about the prediction system are available in Luo et al. (2005). This system has demonstrated high performance for prediction of ENSO (Jin et al., 2008). In particular, Luo et al. (2008) showed that several ENSO events can be predicted at lead times of up to 2 years by this system, which can be a strong advantage in this study. The quasi real-time predictions are updated every month and made publicly available from 2005 (see http://www.jamstec. go.jp/aplinfo/sintexf/e/seasonal/outlook.html).
We have analyzed the reforecasting experiments for the 1993-2018 period issued on the first day of March, June, September, and December with about 2-year lead time. The prediction anomalies were determined by removing the model mean climatology at each lead-time over the same period. To evaluate the prediction results, we have used the multi-mission altimeter satellite gridded sea surface heights (SEALEVEL_GLO_PHY_L4_REP_ OBSERVATIONS_008_047; available from http://marine.cop ernicus.eu/services-portfolio/access-to-products/?option=com_ csw&view=details&product_id=SEALEVEL_GLO_PHY_L4_RE P_OBSERVATIONS_008_047) for sea level, the NOAA OISSTv2 (Reynolds et al., 2002) for SST, and the NCEP/NCAR reanalysis data (Kalnay et al., 1996) for atmospheric variables anomalies. The monthly climatologies of these datasets are also calculated by averaging monthly data over the same period, and then anomalies are derived through deviations from those mean climatologies. All anomalies are linearly detrended, which can

Skill Assessment up to 2-Year Lead
The correlation coefficient (Pearson's "r") between two time series of observed and predicted anomalies for each grid points is used as a deterministic prediction skill score of the phase variation. The persistence method is the simplest way of producing a forecast; it assumes that the conditions will not change. The persistence method works well when anomalies vary very slowly. Therefore, the correlation of the persistence is generally used to assess the advantage of prediction models. Skillful prediction of sea level is found mainly in the tropical Pacific (Figure 1).
It drops outside of the oceanic Kelvin, Rossby, and coastally trapped waveguides in the tropical Pacific region. The correlation often exceeds the skill of persistence and 0.6 in many regions within 20 • of the equator at the first season (0-2 months lead). The correlation decreases at longer lead times but generally remaining above 0.5 in the waveguides up to 11 months lead times. The advantage of the SINTEX-F prediction relative to the persistence increases at longer lead times. This suggests that the skill is mainly derived from the ability to predict ENSO accurately as expected from the previous works (Miles et al., 2014). It is also found skillful prediction regions off the west coast of Australia and California, which may be related to the successful predictions of the Ningaloo Niño/Niña (Doi et al., 2013(Doi et al., , 2015a and the California Niño/Niña (Doi et al., 2015b); some of those events are strongly linked with coastally trapped ocean waves forced by ENSO events. We can also find the seasonal dependence of the correlation, likely because of the so-called spring predictability barrier in forecasting the development of ENSO (e.g., Latif et al., 1998). Prediction of sea level anomalies in the 3-5 months lead average issued on March 1 shows low correlation in many regions relative to that issued on other seasons. Also, we see a quick decrease in the correlation from prediction of the December-January-February (DJF) average into that of the March-April-May (MAM) average.
At the extratropical latitudes in 20 • -30 • N, the skill still remains in some regions in the North Pacific up to about 2-year lead. This is likely due to the slowly propagating Rossby wave features and some stationary anomalies. In the extratropics of the North Pacific, which is often defined the latitude bands of 30 • -60 • N, interestingly, we found the predictability of a region in the North Pacific (30 • -40 • N, 180 • -150 • W) to be skillful up to 2 years ahead. This is not yet discussed by the previous work. Figure 2 shows its time series for the MAM seasonal average. The time series at first glance shows the presence of a decadal variability. The correlation of the 2-year lead prediction is 0.67 for 26 samples (1993-2018 years). After removing the 5-years running mean, it is 0.65 for 22 samples (1995-2016 years). The spread of the ensemble of prediction provides information about the uncertainty inherent. The large uncertainty suggests low potential predictability of sea level anomalies here. Interestingly, however, the sea level anomaly in 2000 exhibits relatively high predictability. We will discuss the details later.
How about the other oceans? We can find some skillful prediction regions in the southern Indian Ocean beyond 1-year lead (Figure 1). This may be related to the successful prediction of the Subtropical Indian Ocean Dipole (Beherea and Yamagata, 2001;Yuan et al., 2014). In the Atlantic, prediction beyond 1-year lead is relatively challenging. For example, the pattern associated with the North Atlantic Oscillation, which is the dominant climate mode in the Atlantic Ocean, is not represented well by the model. We need further analysis to understand similarities and differences among the ocean basins.

Inter-Ensemble Members Relationship
Why is the skillful prediction found about 2-year lead in a region in the North Pacific (30 • -40 • N, 180 • -150 • W)? Investigating co-variability of inter-member anomalies (defined as deviations from the ensemble mean) may provide useful insights into possible precursors and teleconnection patterns related to a climate event considering the intrinsic variability (Ma et al., 2017;Ogata et al., 2019;Doi et al., 2020a,b). Figure 3A shows the correlation coefficients among the inter-ensemble members of the reforecast for the March-May average of 1995-2016 (198 sample: 9 members times 22 years after removing 5-years running mean) at 2-year lead. In this analysis, the conventional time dimension could be enlarged by the ensemble dimension. The covariability between the sea level anomaly in that region and the tropical Pacific condition shows a pattern resembling a mixture of the Modoki-type and the canonical-type of ENSO (Ashok et al., 2007;Karnauskas, 2013). Also, a similar co-variability is seen between the sea level and local wind-induced Ekman downwelling in that region ( Figure 3B). Since the similar features are able to be captured by liner regression analysis to ENSO (Vimont, 2005;Zhang and Church, 2012;Han et al., 2019), the successful prediction of ENSO and/or ENSO-Modoki in the tropical Pacific may be related to the success in predicting sea level anomaly in that region at about 2-year lead. Note that a corresponding co-variability with the surface heat flux was not found in that region ( Figure 3C). This may suggest that the dynamic process is more important in that region than the thermodynamic process.  We see a similar relation with SST anomalies when an ensemble mean of the prediction was considered (Figure 3D), which shows the horizontal map of the correlation coefficients between the ensemble mean of the sea level anomalies averaged in the box and the ensemble mean of the SST fields (22 sample: 22 years of the reforecast). However, the correlation with local wind-induced Ekman downwelling does not show a clear relation ( Figure 3E). Since the signal-to-noise ratio is relatively low in the mid-latitude atmosphere (Scaife and Smith, 2018), the sample size of 22 may not be enough to capture the signal reasonably in the ensemble mean.

Case Study for the 2000 Event
The successful prediction of the high sea level in 2000 (Figure 2) demonstrates the model's high skill to predict such events. As in the observations, the sea level anomaly developed from boreal summer of 1998 and reached about 7 cm during March-May 2000. The prediction issued on June 1, 1998, captured the subsequent development in the tropical Pacific and the target region (30 • -40 • N, 180 • W-150 • W), albeit a bit weaker than in observations (Figures 4A,B). This is supported by the pattern correlations shown in the upper right corner of each panel of Figure 4B, which are calculated after interpolating the horizontal distributions of the observational data to those of the prediction output. At 3-5 months lead, the pattern correlation for the sea level prediction is 0.70. At this time, a La Niña Modoki was observed (Figures 4C,D) in the tropical Pacific. Then, the pattern correlation reduced at longer lead time, however it is still 0.47 at 21-23 months lead time. Local processes seem to contribute to the variability in 30 • -40 • N, 180 • -150 • W relative to remote processes such as a propagation of Rossby waves. Dynamic process associated with the local wind-driven Ekman downwelling may be responsible for that (Figures 5A,C), while the heat flux anomaly acted as the damping of the anomaly in the reanalysis ( Figure 5B) and showed very weak values in the model (Figure 5D). Those features are consistent with the results shown by the previous subsection. We note that the signal in the Kuroshio Extension region was not represented well in this prediction system. Nonaka et al. (2016) revealed that stochastic variability in that region limits deterministic potential predictability of its interannual variability through three-member ensemble simulations with an eddy-resolving ocean model. Even if the spatial resolution of the SINTEX-F is enhanced, it might be intrinsically difficult to improve the prediction skill in that region.

DISCUSSION
Sea level anomalies in the region of (30 • -40 • N, 180 • -150 • W) may be related to the PDO. The PDO is now interpreted as an empirical mode, which includes teleconnection from ENSO and stochastic atmospheric/oceanic fluctuations (Schneider and Cornuelle, 2005;Newman et al., 2016). Decadal or longer timescale signals appear also to be important for the 2000 event. About 50% of the sea level anomaly averaged over 30 • -40 • N, 180 • -150 • W during March-April-May 2000 is due to the decadal signal in the prediction (Figure 2). Actually, 2000 is an extreme year for decadal variations in the IPO index and also basin-wide sea level in the Pacific (see Figure 3 in Lyu et al., 2017), which is closely related to the ENSO-like low-frequency variability. Although low-frequency sea level variations with periods longer than interannual time scales are interesting, it might be difficult to clearly separate the interannual variations from the decadal and longer timescale variations mainly due to limitation of the sample size and the lead time of the reforecast experiments. Although the focus of this study is on seasonal-to-interannual scale prediction, in the future, we may need to develop skillful seamless prediction abilities from seasonal-to-decadal (S2D) timescale.
Tropical and North Pacific processes are interlinked, which means that the North Pacific processes might also contribute to ENSO predictability (e.g., Ogata et al., 2019). The build-up of subsurface ocean heat content in the tropical western Pacific as well as the northeastern subtropical Pacific is identified as ENSO precursors (Capotondi et al., 2015;Yu and Paek, 2015). Chikamoto et al. (2015) also showed that the low-frequency trans-basin tropical climate variations between the Pacific and the other two adjacent ocean basins can be predicted up to 3 years ahead. Further studies are necessary to estimate the role of the inter-basin coupling on multi-year predictability of the tropical and North Pacific using partial assimilation reforecast experiments.
Since our results are based on a single-model system, we need to check them by a multi-model ensemble system (e.g., Kirtman et al., 2014;Tompkins et al., 2017;Widlansky et al., 2017).
A noble path to systematic and synergistic prediction of marine ecosystem variations may be to develop an earth system model, to incorporate biogeochemical processes into a climate model to represent the interacting physical, chemical, and biological processes. It can provide outlooks for marineresource-relevant changes beyond physical variables. Along this line, Park et al. (2019) showed that an earth system model can skillfully predict seasonal to multiannual chlorophyll fluctuations in many regions.
Although predictability of open-ocean anomalies was the focus in this study, its connection to coastal sea level is also important. However, it is still challenging to resolve the complicated topography near the coastal regions for the resolution used in current climate models. Therefore, downscaling techniques are helpful to capture the open-ocean and coastal region connections (e.g., Jacox et al., 2020) in a manner similar to successful examples for atmospheric downscaling (e.g., Ratnam et al., 2016Ratnam et al., , 2017. Enhancement of the relatively coarse ocean model grid will help to resolve more accurately some islands and narrow upwelling regions. Higher resolution in the atmospheric model may also help to improve winds that are an important component of the ENSO teleconnection. In addition, the accuracy should be improved by better initial conditions by explicit use of altimeter data and in situ subsurface ocean temperature and salinity observation from the expendable bathythermographs (XBTs), mooring buoys, sea stations, Argo floats, etc. Increasing the ensemble size may be beneficial for improving prediction of the extratropics, where the signal-to-noise ratio is relatively low. Actually, we have been developing the new version of the SINTEX-F prediction system called as SINTEX-F2 based on a high-resolution model by updating the initialization scheme and increasing the ensemble size (Doi et al., , 2020aMorioka et al., 2019). However, because the computational cost is expensive, the SINTEX-F2 mainly targets for prediction up to 11-month lead time at this stage. We are now extending the lead time up to 23 months because its benefit was shown in this paper.

CONCLUSION
We assessed the prediction skill of sea level anomaly up to 23 months in advance by the SINTEX-F system and found a skillful prediction region in the North Pacific (30-40 N, 180-150 W) at about 2-year lead. The successful prediction of the long-lasting ENSO/ENSO-Modoki in the tropical Pacific seems to contribute to that sea level predictability. The result may be useful to attempt systematic and synergistic prediction of marine ecosystem responses to regional and global climate variations.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the online repositories. The names of the repository/repositories and accession number(s) can be found below: https://fbox. jamstec.go.jp/public/3vnQgARS2EDAjpEBzL10R3qP8jcy2W-xWLfXy5mE7u8Z/m/3x5AD2lA.

AUTHOR CONTRIBUTIONS
TD performed the seasonal prediction experiments, and analyzed the observation data and model prediction outputs. All authors contributed to designing the research, interpreting results, and writing the manuscript.

ACKNOWLEDGMENTS
The SINTEX-F seasonal climate prediction system was run by the Earth Simulator at JAMSTEC (see http://www.jamstec.go. jp/es/en/index.html, for the system overview). We are grateful to Wataru Sasaki, Jing-Jia Luo, Sebastian Masson, Antonio Navarra, and Toshio Yamagata and our European colleagues of INGV/CMCC, L'OCEAN, and MPI for their valuable contributions in developing the prototypes of the systems. We would also like to thank Shoshiro Minobe for his helpful comments and suggestions. The GrADS software was used to create the figures and maps.