Oceanic Rossby waves drive inter-annual predictability of net primary production in the central tropical Pacific

In the Pacific Ocean, off-equatorial Rossby waves (RWs), initiated by atmosphere-ocean interaction, modulate the inter-annual variability of the thermocline. In this study, we explore the resulting potential gain in predictability of central tropical Pacific primary production, which in this region strongly depends on the supply of macronutrients from below the thermocline. We use a decadal prediction system based on the Max Planck Institute Earth system model to demonstrate that for the time period 1998–2014 properly initialized RWs explain an increase in predictability of net primary productivity (NPP) in the off-equatorial central tropical Pacific. We show that, for up to 5 years in advance, predictability of NPP derived from the decadal prediction system is significantly larger than that derived from persistence alone, or an uninitialized historical simulation. The predicted signal can be explained by the following mechanism: off-equatorial RWs are initiated in the eastern Pacific and travel towards the central tropical Pacific on a time scale of 2–6 years. On their arrival the RWs modify the depths of both thermocline and nutricline, which is fundamental to the availability of nutrients in the euphotic layer. Local upwelling transports nutrients from below the nutricline into the euphotic zone, effectively transferring the RW signal to the near-surface ocean. While we show that skillful prediction of central off-equatorial tropical Pacific NPP is possible, we open the door for establishing predictive systems for food web and ecosystem services in that region.


Introduction
In Earth's climate system, dynamic interaction of atmosphere and ocean leads to the initiation of oceanic waves in a wide range of spatio-temporal scales (Gill 1982). Depending on their characteristics, these waves may carry memory within the Earth's climate system, and can therefore be used in skillful climate predictions. Off-equatorial oceanic planetary waves (Rossby 1939(Rossby , 1945, commonly known as oceanic Rossby waves (RWs) (RWs; Platzman 1968), could carry memory on the interannual to decadal time scales (2-10 years). Oceanic RWs occur with typical phase speeds of 0.03-0.3 m s −1 (approx. 1000-10 000 km/year) and wavelengths larger than 500 km (Chelton andSchlax 1996, Killworth andBlundell 2003). The vertical displacement due to oceanic RWs varies from only a few centimeters at the ocean's surface (Chelton and Schlax 1996) to several tens of meters at thermocline depth (Kessler 1990). Satellite altimetry confirmed the abundance of oceanic RWs in all of Earth's major ocean basins (Chelton and Schlax 1996). However, the longest memory resulting from oceanic RWs can be expected to be found in the Pacific, where they cross the ocean on time scales of 5-30 years (Capotondi andAlexander 2001, Galanti andTziperman 2003). Here, oceanic RWs are predominantly initiated in the eastern portion of the Pacific Ocean, involving atmospheric and oceanic processes, e.g. Ekman pumping connected with surface winds (Capotondi and Alexander 2001), or modulations of the pycnocline in connection with El Niño-Southern Oscillation (ENSO, Vega et al 2003). Propagating RWs may have an impact not only on the interannual variability of the ocean's dynamical state, such as temperature, salinity, sea surface height, but also on nutrient availability, and related to that, biological primary productivity (Siegel 2001, Uz et al 2001.
At ocean depths below 100m, where the absence of Sun light prohibits biological production, variability of nutrient concentration is almost negligible on inter-annual time scales (Uz et al 2001). At shallower depths, in the euphotic zone, the availability of nutrients varies more than at larger depths, and may determine the strength of biological net primary productivity (NPP) (Barber 1992), in particular in tropical oceans. Vice versa, biological production leads to a depletion of nutrients in the euphotic zone and as a result a steep gradient in nutrient concentration with depth is formed, the nutricline. Oceanic RWs modulate the depth of the nutricline by a few tens of meters, with corresponding impacts on the nutrient availability above the nutricline, and on NPP and the resulting chlorophyll concentration (Uz et al 2001. Compared to climatology, a shallower (deeper) than normal nutricline leads to larger (smaller) than normal primary production. Oceanic RWs impact surface chlorophyll concentration also by a vertical displacement of the chlorophyll maximum (Kawamiya and Oschlies 2001). Signatures of oceanic RWs have therefore not only been observed in sea-level anomalies, but also in ocean color attributed to oceanic chlorophyll (Uz et al 2001, Gutknecht et al 2010. In addition, RWs could impact the variability of surface chlorophyll and organic particles by lateral advection, without a concurrent change in nutrient availability (Dandonneau et al 2003, 2004, Killworth 2004. While there is a general consensus that advection due to oceanic RWs directly modulates upper ocean chlorophyll concentrations on the interannual time scale, the degree of the impact of oceanic RWs on nutrient availability above the nutricline and a co-varying NPP remains unclear , Charria et al 2008. It still remains to be shown how the inter-annual memory potentially carried by oceanic RWs can be used for skillful decadal predictions of nutrient availability and NPP. Here, we set out to precisely fill in this gap, by investigating RWs in a decadal predictions system. On the basis of comprehensive Earth system models (ESMs), decadal prediction systems have been developed (see, e.g., Boer et al 2016) with the aim to predict Earth's climate, in particular the dynamical state of atmosphere and ocean, on inter-annual to decadal time scales. In the slipstream of maturing decadal prediction systems, some studies assessed the inter-annual prediction of biogeochemical parameters. In a groundbreaking study, Séférian et al (2014) found predictability of NPP in the tropical Pacific with an ESM for up to three years in advance, attributed, amongst other drivers, to the poleward advection and persistence characteristics of nutrients, driven directly by ENSO. Ham et al (2021) found that in the southern tropical Pacific, ENSO variability has a major impact on the advection of nutrients and thus on surface chlorophyll concentration, limiting skillful predictions to the seasonal time scale. For the North Atlantic sub-polar gyre, Li et al (2016) show that the CO 2 uptake by the ocean could be predicted 4-7 years in advance. On the global scale, Park et al (2019) found skill in their ESM-based seasonal predictions of chlorophyll concentrations, also mainly connected to ENSO, and Fransner et al (2020) showed that oceanic biochemical predictability is closely tied to the representation of the upper ocean dynamical state in their ESM.
Despite the success story of decadal prediction systems, the role of oceanic RWs in the predictability of tropical nutrient availability and NPP has not been comprehensively studied so far. Bridging the gap between dynamical and biogeochemical prediction remains an ongoing challenge. In this study, we use a comprehensive Earth system model to pick up this challenge with the goal to successfully predict NPP in the tropical Pacific several years ahead and to pinpoint the role of off-equatorial oceanic RWs in carrying and distributing inter-annual memory in the tropical Pacific.

Model simulations, observational references and prediction skill
We use the Max Planck Institute (MPI) for meteorology ESM at low resolution (MPI-ESM-LR, Giorgetta et al 2013) with the atmospheric component ECHAM6 (Stevens et al 2013) in T63L47 configuration, equivalent to 2.8 • horizontal resolution, and the oceanic component Max Planck Institute Ocean Model (MPIOM) (Jungclaus et al 2013) in GR15L40 configuration with 1.5 • nominal horizontal resolution. In MPI-ESM, oceanic biogeochemistry is represented by the Hamburg ocean carbon cycle model (HAMOCC, Ilyina et al 2013). All simulated values of NPP and nitrate concentration used in this study are calculated by HAMOCC, sea surface temperatures (SST) are taken from MPIOM. We use the vertical mass transport simulated in MPIOM, averaged over 200-450 m depths, to analyze the characteristics of oceanic RWs.
Our prediction system consists of 10-year 8-member ensemble hindcasts with MPI-ESM-LR, initialized from a weakly coupled data assimilation experiment at the 1st of January every year from MPI-ESM-LR assimilation covers the time period 1958-2014, where we incorporate ERA-40/ERA-Interim re-analyses of vorticity, divergence, temperature and sea level pressure (Uppala et al 2005, Dee et al 2011 via atmospheric nudging and observational profiles of oceanic temperature and salinity (EN4, Good et al 2013) using a non-localized 8-member oceanic ensemble Kalman filter (EnKF, Brune et al 2015). In this weakly coupled assimilation setup, the ocean responds mainly to the atmospheric nudging with weak changes added by the oceanic EnKF. We compare the initialized ensemble hindcasts with an 8-member ensemble of uninitialized MPI-ESM-LR historical simulations. Both the initialized and uninitialized simulations are run under identical prescribed external conditions, e.g. green house gas concentration, solar radiation, and volcanic aerosols, according to CMIP5 (Taylor et al 2012 (HadISST, Rayner et al 2003). For supplemental analysis (not shown) we use AVISO as our reference data for sea surface height, this altimeter products were produced by Ssalto/Duacs and distributed by Aviso+, with support from Cnes (www.aviso.altimetry.fr).
Inter-annual prediction skill is measured using the correlation coefficient and the root mean square error between the respective yearly means of simulations and reference data. Prior to skill analysis, linear trends are removed. Significance of correlations is calculated at the 95% level by analyzing both tails of the 1000-member bootstrapped probability distribution (not shown) of the respective simulations with a p-value of 0.025 (Efron 1979, Wilks 2011. In the central tropical Pacific, ENSO has a large direct impact, also on decadal time scales (C-Mode in Takahashi et al (2011)), which may mask the oceanic RW signal in oceanic surface quantities. In an effort to carve out the oceanic RW signal, we therefore remove the direct impact of ENSO using linear regressions first on the central Pacific Niño3.4 (5 • N-5 • S, 170 • W-120 • W) SST index and then on the eastern Pacific Niño1+2 (0 • -10 • S, 90 • W-80 • W) SST index (Trenberth and Stepaniak 2001).

Representation of off-equatorial oceanic Rossby waves
In general, Pacific off-equatorial RWs are abundantly simulated in the hindcasts of our MPI-ESM prediction system (figure 1). In accordance with observations (e.g., Chelton and Schlax 1996), simulated RWs emanate in the tropical Pacific from two distinct regions North and South of the equator: in the North Pacific around the South-Eastern edge of the sub-tropical gyre, in the South Pacific around the Peru-Chile current off-shore South America ( figure 1(a)). Two distinct RW trains extend from these source regions to the West: in the North Pacific along 15 • N, and in the South Pacific along 10 • S. Around these latitudes, the RWs modulate the depth of the nutricline (figures 1(b) and (c)) and the thermocline (not shown). The longitudinal extension of source regions and RW trains along these latitudes does not change over the entire time period 1998-2014 (figures 2(a) and (b)). At 15 • N (figure 2(a)), the characteristic period of RWs is about 2 years, the travel time from the source region around 140 • W across the basin towards the maritime continent is about 2-4 years. At 10 • S ( figure 2(b)), the characteristic period is 2-6 years, while the travel time from the source region off-shore South America towards the central tropical Pacific is about 2-6 years. Thus, we determine our target regions in the tropical Pacific as those regions, which show a large impact of oceanic RWs with a characteristic time period of 2-6 years.
Filtering the vertical mass transport with a characteristic RW period of 2-6 years reveals in the Pacific two bands on each side of the equator with large variance due to RWs (figure 2(c)). In each of these bands, the easternmost part represents the source region of RWs. They travel westward and thus defining the target regions for our analysis of prediction skill. In the tropical North Pacific, the source region can be broadly defined as a box at 120 • W-140 • W, 10 • N-20 • N, and the target region as a box 150 • E-170 • W, 10 • N-20 • N (henceforth referred to as North Pacific box). In the tropical South Pacific the source region is located at around 85 • W-100 • W, 20 • S-10 • S, and the target region at 160 • W-110 • W, 15 • S-5 • S (henceforth referred to as South Pacific box). In both source regions, RWs in our hindcasts are properly initialized in time and space by the assimilation, indicated by high correlations of the simulated sea surface height anomaly with AVISO observations (not shown, r = 0.80 and r = 0.85 for the Northern and Southern Pacific source regions, respectively). Hence, in the target regions we expect a significant contribution of RWs not only to the variability of the vertical mass transport around the thermocline, but also to the modulation of the nutricline

RW impact on prediction of nutrient availability and NPP
The predicted nutrient availability for up to 4 years in advance is modulated by the inter-annual variability of the strength of predicted off-equatorial RWs in the target regions. Variability of nitrate concentration is strongly influenced by RWs not only at depths below 100 m (figures 3(a) and (b)), but the RW signal is also transferred from these depths to the euphotic zone, i.e. the signal can be traced over all depths and often reaches the surface (figures 3(a) and (b)).
In both target regions, the year-to-year and the multi-year variability in nutrient concentration and NPP fit well with the RW signal seen below 100 m depth. In the euphotic zone, there is significant correlation between the predicted nutrient concentration and predicted NPP 4 years in advance (detrended and direct ENSO impact removed) in both the North Pacific box (figure 3(c), r = 0.57) and the South Pacific box (figure 3(d), r = 0.75). Simultaneously, the predicted NPP is significantly correlated with the observed NPP (detrended and direct ENSO impact removed) for 4 years in advance in both the North Pacific (figure 3(c), r = 0.68) and South Pacific box (figure 3(d), r = 0.76).
The correlation between predicted and observed NPP depends on the lead year. For the North Pacific box, we find significantly positive correlation skill for NPP in our prediction system for lead times up to 4 years, and for the South Pacific box for lead times up to 5 years (see filled circles in figure 4). In this time frame, the time-lag auto-correlation of observed NPP alone (NPP persistence) hardly leads to any significantly positive correlation skill, indicating a clear benefit of using MPI-ESM over persistence for NPP prediction. Correlation skill of NPP in a corresponding historical simulation, i.e. without transferring observed states of atmosphere and ocean to MPI-ESM, is negative in both Pacific boxes. The evaluation of the root mean square error (RMSE) against reference corroborates the findings of our correlation analysis. In the North Pacific box, the RMSE is for all lead years smaller for the predicted NPP than for NPP persistence, and for up to lead year 4 the RMSE for predicted NPP is also smaller than for the historical simulation. In the South Pacific box, RMSE for the predicted NPP is for all lead years smaller than for NPP persistence of NPP from the historical simulation. For predictions up to 5 years in advance, there is a substantial gain in prediction skill when assimilating observed atmospheric and oceanic states into MPI-ESM prior to starting the hindcasts. Thus the proper representation of oceanic RWs, allowed by the assimilation of observed atmospheric and ocean states at initialization, increases the predictability of NPP at the inter-annual time scales regions with a large impact of RWs.

Comparison of NPP and SST predictability
In our prediction system, correlation skill of NPP in both the North Pacific and South Pacific boxes is consistently higher than that of SST for lead years 2-5 ( figure 4). An important reason for this lies in different susceptibilities of NPP and SST to the local atmospheric surface conditions in the tropical Pacific. Near-surface ocean temperatures are directly coupled to the local atmospheric conditions, which actively  impact SST and add a chaotic component to its prediction, eventually lowering prediction skill. In contrast, NPP in the tropical Pacific is not as strongly affected by local atmospheric conditions as SST, since temperature is not a dominant factor to NPP in the tropical oceans in general (passive impact). In turn, inter-annual predictions of NPP in the North Pacific and South Pacific boxes benefit more from the inter-annual variability represented by properly initialized oceanic RWs than SST.
Looking beyond the selected regions, NPP is predictable up to lead year 4 over a substantial portion of the off-equatorial tropical Pacific ( figure 5). There is a strong indication that the mechanism proposed in our study is generally complemented by other processes, which eventually lead to predictability of NPP on the inter-annual time scale on a much larger geographic scale. In particular, our prediction system can skillfully predict NPP along the equatorward edges of the Pacific subtropical gyres. While the gyre boundary may act as a wave guide to oceanic RWs, which eventually modulate NPP variability, the correct prediction of the contrasting dominating effects of lateral and vertical nutrient advection outside the gyres, and strong stratification and nutrient recycling inside the gyres (Séférian et al 2014, Karl andChurch 2017) could also be the reason for skillful NPP predictions in our decadal prediction system.
Our results also reveal that oceanic RWs do not always imply skillful NPP predictions. Oceanic RWs traveling within the Pacific subtropical gyres do not seem to, if at all, positively impact inter-annual NPP predictions (e.g. figures 5(e) and (g)), most probably due to the dominance of the nutrient recycling regime within the gyres. Also, prediction skill decreases beyond lead year 5, indicating an increase in the unpredictable components in the climate system, counteracting the memory potentially carried by oceanic RWs.

Discussion and conclusion
The potential for inter-annual predictability of biogeochemical parameters in the Pacific has been evaluated previously by Séférian et al (2014). While using an early-state decadal prediction system (e.g. assimilating only SST), they found that prediction skill for NPP could be high for up to 3 years in advance. In our study, we use a more comprehensive decadal prediction system, assimilating 3D atmospheric and ocean states in a model-consistent approach, to specifically examine one potential dynamical source of inter-annual predictability in the tropical Pacific, oceanic RWs, while removing the direct impact of another potential source, ENSO. We find that in regions with a potentially large impact of oceanic RWs, the prediction horizon of NPP could be extended by 1-2 years over the 3 years found in Séférian et al (2014) for the tropical Pacific, with a major contribution by the successfully timed initialization and propagation of oceanic RWs. Although our study highlights the impact of the vertical supply of nutrients in connection to oceanic RWs, other dynamical processes are known to impact nutrient supply and primary production. Killworth et al (2004) discuss that horizontal advection of nutrients, also in connection with oceanic waves, may have a strong impact on both nutrient supply and chlorophyll production. Also, nutrients could be transported with oceanic currents, in particular along gyre boundaries, and thus impact NPP on the seasonal to inter-annual time scale. In our study we did not directly account for the effect of nutrient recycling, which has been discussed to impact primary production, in particular inside oceanic gyres (Séférian et al 2014, Karl andChurch 2017).
Our initialization approach relies on a rather simple weakly coupled assimilation of both atmospheric and oceanic observational products. It remains to be shown, if other state-of-the-art initialization approaches, such as discussed in Brune and Baehr (2020), Merryfield et al (2020) or Fransner et al (2020), could lead to further improvements in the prediction of NPP.
In our analysis, we rely on satellite-based ocean color, i.e. chlorophyll, observations as our reference for NPP observations. Earlier studies pointed out that although oceanic RWs could account for up to 30% of the modulation of regional NPP (Kawamiya and Oschlies 2001), the associated impact on satellitebased chlorophyll measurements could be smaller, Uz et al (2001) indicate a range of 5%-20%. In addition to in-situ primary production, visible chlorophyll concentrations are also directly dynamically related to transports within the oceanic mixed layer (Kawamiya and Oschlies 2001). We acknowledge that satellite-based instruments not only changed over the time period 1998-2014, but also may encounter sensor decay over their life-time. This contributes to an increased uncertainty on the reference side in our NPP predictions. To tackle this in future studies it could be advisable to revisit NPP references with Sentinel-3 data (Donlon et al 2012).
Although we removed the direct impact of ENSO prior to our analysis of NPP predictability, it should be noted that ENSO is indirectly involved in predicting NPP with MPI-ESM. ENSO determines, among other processes, the atmosphere-ocean interaction in the source regions of oceanic RWs. Therefore, the assimilation of the correct ENSO state at the start of the prediction is a key factor to take advantage of the multi-year memory offered by oceanic off-equatorial RWs in the tropical Pacific, and is the main source of our increased prediction skill in NPP in the initialized predictions when compared to uninitialized historical simulations.
From our analysis of tropical Pacific inter-annual prediction skill of NPP for the time period 1998-2014 in MPI-ESM, we conclude: • In the tropical Pacific, we find high prediction skill of NPP for up to 5 years in advance, with a large contribution from properly initialized and simulated off-equatorial oceanic RWs. • The prediction skill of NPP represents a significant increase on prediction skill of SST, which is limited to 1-3 years. We argue that this difference arises from the different impact of near surface atmospheric dynamics on nutrient concentration (passive impact) and near surface ocean temperature (active impact).
The connection of dynamical and biogeochemical predictions is encouraging for further studies on the predictability of the food web in the tropical Pacific ocean.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://cera-www.dkrz.de/WDCC/ui/cerasearch/ent ry?acronym = DKRZ_LTA_1075_ds00006 (Brune et al 2021).