Strengthened ocean-desert process in the North Pacific over the past two decades

North Pacific ocean desert (NPOD) refers to the subtropical North Pacific Ocean of low chlorophyll-a (Chl-a) concentrations, as the largest ocean desert globally. Studies have suggested a development of NPOD over recent decades based on limited evidences from in-field measurements and yet elusive mechanism. In this study, we characterize intensity, area and position of the NPOD from year 1998 to 2018, and investigate its control by the coherent climate processes, based on an available, longest satellite observations of Chl-a concentration. Our results suggested that NPOD oligotrophication and expansion processes were correlated with warming upper oceans in most part of the NPOD, except for the SW NPOD area where the Chl-a variations were linked with regional change in sea surface heights. Moreover, based on our analysis, insignificant shift but only NW-SE variability of the NPOD mean position was likely controlled by the Pacific decadal oscillation processes.


Introduction
Satellite observations of surface-ocean chlorophyll-a (Chl-a) concentrations reveal the existence of five ocean deserts, globally including the North Pacific ocean desert (NPOD), South Pacific ocean desert, South Indian ocean desert (SIOD), North Atlantic ocean desert (NAOD) and South Atlantic ocean desert (figure S1(a) (available online at stacks.iop.org/ERL/16/024034/mmedia)). These ocean desert areas are characterized according to their oligotrophic feature, i.e. Chl-a concentrations less than 0.07 mg m −3 (Polovina et al 2008). Geographically, the five ocean deserts occupy ∼40% of the low-latitude ocean area (McClain et al 2004) and play the roles in significantly changing the Earth ecosystem and carbon cycles, thus to change atmospheric CO 2 and the ongoing global warming (Sarmiento et al 1998, Behrenfeld et al 2006, Doney et al 2014.
Moreover, modelling studies have projected a continuous decrease in the Chl-a concentrations and primary production (PP) commonly in ocean desert areas during the coming centuries (Henson et al 2010, Steinacher et al 2010, Hofmann et al 2011, Polovina et al 2011. Overall, studies have highlighted the emergency in understanding the linkage between the change in ocean desert and the global warming processes. As the largest ocean desert globally, NPOD has shown the decline in Chl-a concentrations, i.e. oligotrophic intensification process, during the periods of year 1997-2003(McClain et al 2004, Gregg and Rousseaux 2005, year 1998-2007 (Henson et al 2010), year 1998-2010 (Signorini and McClain 2012), year 1998-2012 (Gregg and Rousseaux 2014) and year 1998-2013 (Signorini et al 2015). In parallel, studies have suggested that the NPOD area becomes larger, i.e. expansion, during year 1997-2010 (Polovina et al 2008, Irwin andOliver 2009). All these NPOD variations have been linked to either climate natural variability, anthropogenic forcing or their combined control. As a most affective factor, the nutrients budget into the upper oceans regulates the Chl-a variations via several dynamic processes, including advection transport, vertical mixing in the ocean column and atmospheric deposition at the air-sea interface (Patra et al 2007, Xiu and Chai 2012, Liu and Levine 2016, Gittings et al 2018, Yang et al 2020. In addition, both the decrease of surface Chl-a concentration and the ocean desert expansion have been correlated to higher sea surface temperatures (SST, Gregg and Rousseaux 2005, Polovina et al 2008, Hoegh-Guldberg and Bruno 2010, Boyce et al 2012, Lewandowska et al 2014. In parallel, higher SSTs have been also suggested to cause enhanced stratification and thus nutrient-limitation in upper ocean during the climate warming. These processes determine the NPOD oligotrophication (Chl-a concentration declining process) together with NPOD area expansion (Polovina et al 2008, Lewandowska et al 2014. Moreover, decadal variability in the NPOD area and intensity have been observed (Henson et al 2010) and likely linked to the Pacific decadal oscillation (PDO, Irwin andOliver 2009, Martinez et al 2009) and the El Nino-Southern oscillation (ENSO, Chavez et al 2011, Racault et al 2017, Park et al 2018. On the other hand, the variations of NPOD apparently rely on the change in the ocean PP seasonality (Henson et al , 2017, the phytoplankton photoacclimation response (Behrenfeld et al 2016) and the phytoplankton community structure (Lewandowska et al 2014).
In this study, based on a combined, longest dataset of Chl-a concentration derived from the Sea-WiFS and MODIS satellite observations, we investigate the deseaonalized variations of NPOD during year 1998-2018. In particular, three key characteristics of NPOD, i.e. position, area and intensity, and their alignment via climate physics are discussed, meanwhile these aspects have been investigated only individually or in part by previous studies.

Dataset
We characterize the NPOD characteristics based on the surface-ocean Chl-a concentration derived from the ocean color data by SeaWiFS and MODIS-Aqua satellite observations by NASA Ocean Biology Processing Group (OBPG 2014). Level-3 global standard mapped images of daily Chl-a datasets of 9 km spatial resolution are used between January 1998 and December 2007 from the SeaWiFS dataset while for the MODIS data between January 2003 and December 2018 (figure S1(c)). A 13 d running-mean and spatially 100 km smoothing are applied to both datasets to fill in the cloud-covered blank and to avoid the impact of the mesoscale processes, respectively. Using the overlapping period of 2003-2007 in Sea-WiFS and MODIS-Aqua datasets, we have combined the two Chl-a concentration data into a continuous one from January 1998 to December 2018, by applying the linear regression approach on each paired data image pixel and at monthly resolution (figure S1(b)). Overall, the two datasets present an R-squared coefficient of 0.925 with offset up to 0.002 mg m −3 , suggesting a strong correlation between the SeaWiFS and the MODIS datasets (figure S1(b)).
To understand the mechanism in the NPOD variation due to climate physics, we have performed the analysis about the SST data from NOAA optimum interpolation SST data version 2 (Reynolds 1988), photosynthetically active radiation (PAR) data from MODIS (NASA OBPG 2014), aerosol optical depth (AOD) data from MODIS (Levy et al 2015), sea surface height (SSH) data from DUACS, 10 m wind speed (WS) data and precipitation rate (PRE) data from ECWMF (C3S 2017), nutrients (N + P) concentration data from WOA 2005 (Levitus 1983) and the mixed layer depth (MLD) data from SODA3 (Carton et al 2018). Moreover, A 1D K-profile parameterization (KPP) algorithm (Large and Gent 1999) is driven by the horizontal velocity components, temperature, salinity data from SODA3 reanalysis data (Carton et al 2018) to investigate the stratification conditions in NPOD area. More dataset details are shown in table S1.

Definition of NPOD indexes
Our definition of the NPOD is based on its Chl-a concentration less than 0.07 mg m −3 , as defined in McClain et al (2004) and Polovina et al (2008). In addition, the NPOD intensity is represented by the averaged value of Chl-a concentration in the core NPOD region determined by 0.07 mg m −3 contour of climatological mean Chl-a concentration (figure S1(a)). Moreover, we assess the core position of NPOD based on its zonal (α lon ) and meridional (δ lat ) extent using the following equations: where c chl-a is surface Chl-a concentration, n and m are the numbers of grid points where the c chl-a is below 0.07 mg m −3 , l lat and l lon are the latitude and longitude of each grid point, respectively. Here, our method to track the NPOD position is originally evoked by Hsin et al (2013) and has been used to obtain the central position of the North African dust transport route and central position of the Ocean Gyre (Meng et al 2017, Yang et al 2020. In the following, our analysis is based on the calculation results with the p-value less than 0.01 presented, for statistical robustness.

Granger causality test
Granger causality test is a statistical hypothesis test for examining whether a time series X meaningfully depends on another time series Y (Granger 1969). If Y can be better predicted using the histories of both X and Y than using the history of Y alone, namely lagged values of X contain information that helps explain current values of Y, then time series X is said to Granger cause Y (Kaufmann and Stern 1997). Its mathematical formulation is based on linear regression modeling of stochastic processes (Granger 1969). To test the null hypothesis that X does not Grangercause Y, one first does a regression of △X on lagged values of △Y. Once the appropriate lag interval for Y is proved significant through t-test or p-value, the regression is augmented with lagged levels of △X. If (1) it is significant according to a t-test and (2) it and the other lagged values of △X jointly add explanatory power to the model according to an F test, any lagged value of △X is retained in the regression. Then the null hypothesis of no Granger causality is accepted if and only if no lagged values of △X have been retained in the regression. In this study, the Granger causality test was used to examine the link between the Chl-a variation and SST, SSH, and before that detrending and unit root test were applied to ensure the validity of the result.

Results and discussion
Our results characterize NPOD by a continuous oligotrophication process at a rate of −2.1 × 10 −4 mg Chl-a m −3 yr −1 (around −0.43% yr −1 ) and its expansion at 13.5 × 10 4 km 2 yr −1 (around 0.9% yr −1 ) from 1998 to 2018, linked to higher SST-induced vertical stratifications and SSH-induced lower nutrients transport. Furthermore, we have tested the effects of SST, SSH (represent climate warming) and WS, AOD, PRE, PAR on the NPOD. Moreover, we find the NPOD position exhibits a northwest-southeast migration, following the PDO index.

Oligotrophication, expansion and migration of the NPOD between 1998 and 2018
As shown in figure 1(a), according to the linear-fitted regressions of deseaonalized Chl-a time serious, oligotrophication (Chl-a concentration reduce) mainly occurred in northwest and south NPDO area at a rate of 0.43% yr −1 during the years 1998-2018.
We notice that, the start year 1998 is a strong El Nino year, so we applied further sensitivity test with different start years (text S1(c) and figure S2). The oligotrophication speed in this study is lower than the value of 0.75% yr −1 in 1998-2010 and 0.71% yr −1 in 1998-2013, as proposed in previous studies McClain 2012, Signorini et al 2015). Here, one possible reason for the change in oligotrophication speed is that the period of reported negative trend of Chl-a concentration is upon the decreasing stage of the oscillation, likely reflecting the impact of nature decadal variability (Henson et al 2010). Following the method of Henson et al (2010), to analyze the oligotrophication rates in different stages, the time series of Chl-a concentration in the NPOD is split into overlapping 10 year sections ( figure S3(a)). In the first 10 year stage (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) and last stage (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), oligotrophication processes are faster compared to the middle stage, based on continuous decrease over the 21 years. The slow-down of the NPOD oligotrophication during year 2005-2014 could be attributed to the cooling of the North Pacific SST associated with warming 'hiatus' (Kosaka andXie 2013, Hu andFedorov 2017).
Along with the oligotrophication process, the NPOD have expanded in size over the past 21 years at a rate of 13.5 × 10 4 km 2 yr −1 (0.9% yr −1 , figure 1(e)). Spatially, the NPOD expansion has mainly occurred at the northwest, south and east borders of the NPOD (figure 1(b)), coherent with the declining in Chl-a concentrations in same region (figures 1(a) and (b)). Here, our calculation of NPOD expansion has applied a definition of NPOD using the Chl-a concentrations lower than 0.07 mg m −3 as the threshold (McClain et al 2004, Polovina et al 2008, meanwhile we also notice that the definition of NPOD and its size depend on this threshold value. Therefore, we have additionally applied sensitivity test with various thresholds in Chl-a concentrations to investigate the sensitivity of NPOD area expansion to these thresholds, as shown in figure S4. Using different thresholds between 0.05 and 0.11 mg m −3 in Chl-a concentration, the NPOD shows expansions at ratios between 35% (0.05 mg m −3 ) to 9%, thus consistently presenting the NPOD expansion and robustly independent from the selective NPOD definitions in our analysis. Moreover, the climatological mean position of NPOD, at around 18 • N, 177 • E, shows interannually north-south oscillations by 1.6-2 • and west-east shift of 7.8-8.4 • during years of 1998-2018 (figures 1(f) and (g)). At each 10 years (figures S3(b) and (c)), the sign of the trends in the NPOD position switches between northwestwards and southeastwards on decadal timescales, implying the potential effect of natural decadal variability on NPOD migration.

Linkage between NPOD and climate factors
The NPOD oligotrophication, expansion and northwest-southeast migration processes are correlated to the nutrient change within upper euphotic layer. Previous studies have linked these NPOD variations to climate processes based on observation datasets (Gregg and  SSTs enhance stratification and thus reduce nutrient upwelling, leading to the expansion of ocean desert in the coming century; (b) weaker wind stress causes less wind energy input and the resultant efficiency of vertical mixing; (c) deeper nutricline and less upwelling reduce the nutrient content in the upper ocean aligned with higher SSHs; (d) anthropogenic atmospheric nutrients deposition changes the nutrients structure in the surface ocean due to dry deposition (related to AOD) and precipitation; (e) the decrease of surface PAR contributes to the lowering in Chl-a concentration. Here, we detect the correlations between SST, SSH, WS, AOD, PRE, PAR and the NPOD variations, based on deseaonalized data. Overall, as shown in figure S5, the Chl-a concentrations show negative correlations with the SST and SSH in the NPOD region, more robust compared to other factors. This is thus in line with the previous illustration about the linkage between Chl-a and SST, SSH in North Pacific (Thomas et al 2012).
Moreover, we have made further analysis about the spatial characteristics of the correlations 'Chl-a vs SST' and 'Chl-a vs SSH' during year 1998-2018. As shown in figure 2(a), the NPOD region can be divided to NPOD SST and NPOD SSH, according to which factor of the maximal correlation coefficient to Chl-a. In most part of NPOD region, the Chl-a variation significantly correlates with the SST, meanwhile the southwest NPOD shows a dependence on SSH features ( figure 2(a)). In the NPOD SSH region, the Chl-a concentrations show variability without a significant trend throughout the 21 year interval. In parallel, the SSH variability in the NPOD SSH area is two times higher than that in the NPOD SST region, also higher than the value of global average (Kobashi andKawamura 2001, Zhang andChurch 2012). Statistically, both these two features act as reasons for SSH correlates well with Chl-a variation in the NPOD SSH area (r = −0.67, p < 0.01). In the NPOD SST region, the Chl-a concentrations show a decline significantly, correlated to the SST variations at r = −0.46 and p < 0.01 (figures 1(a), 2(a) and (c)). This result supports the studies which suggest that warmer upper layer causes more stratification and inhibits the upward nutrients transport in global scale (Polovina et al 2008, Lewandowska et al 2014. However, within the large area of NPOD SST , regional patterns of Chl-a may be also under the impact of other factors, e.g. nutrients lateral transport (Letscher et al 2016) and eddy pumping (Oschlies 2008). Therefore, the correlation between Chl-a and SST in NPOD SST (0.46) is weaker than the value between Chl-a and SSH in NPOD SSH (0.67). Furthermore, our granger causality test results show that the probabilities of the null hypothesis that the SST/SSH derived nutrients flux do not cause Chl-a variation in the NPOD SST /NPOD SSH area were nearly 100% at the 1% significance level when 1-3 months lags were assumed (table S2). Statistically, this can be a further support for the causality link between NPOD variation and SST/SSH.
In addition, our results show a negative correlation (figure S6) between the intensity and area of NPOD, suggesting the occurrence of NPOD expansion coherent with its oligotrophication process during year 1998-2018. In most part NPOD areas, both the oligotrophication process and the NPOD expansion correlate to the geographic part where SST correlates well with Chl-a variation (figures 1(a), (b) and 2(a)). This indicates the existence of a control of warmer upper ocean on the NPOD oligotrophication and coherent expansion.
Furthermore, the NPOD position exhibits northwest-southeast oscillation without a significant MLD (m) Figure 3. Impact of sea surface temperate (SST) variation on the NPOD oligotrophication. (a) Time series of PDO index (light blue, blue for 3 month running-mean) and the deseaonalized time series of SST in NPOD SST (light red, red for 3 month running-mean). (b)-(d) Time series of temperature anomaly (b), salinity anomaly (c) and density anomaly (d) in NPOD SST , and the black dashed lines represent the contours of monthly temperature, salinity and density from 1998 to 2016. (e) Deseaonalized time series of MLD (light blue, blue for 3 month running-mean) and mixture coefficients (KM, light red, red for 3 month running-mean) in KPP algorithm.
trend during year 1998-2018, at a negative correlation of −0.57 (p < 0.01) between the deseaonalized NPOD meridional position and the PDO (please see text S1 and figure S7 for the link between NPOD and ENSO), as well as a negative correlation of −0.63 (p < 0.01) between the deseaonalized NPOD zonal position and the PDO (figures 1(f) and (g)). Generally, NPOD moves southeastward by ∼2 • and ∼7.8 • in PDO warm phases, northwestward by ∼1.6 • and ∼8.4 • in PDO cold phases, respectively. As shown in figure 1(c), the distribution of the correlation coefficients reveals that Chl-a time series over NPOD are positively correlated with the PDO at the northwest boundary but negatively correlated with the PDO at the southeast boundary. Together with spatial pattern in the correlations of 'Chl-a vs SST' , our results indicate a process that lower SSTs increase the Chl-a concentrations at the northernwestern boundary and higher SSTs decrease the Chla concentrations at the southern-eastern boundary based on the PDO warm condition. Consequently, NPDO moves the southeastward in the PDO warm phases meanwhile northwestward in the PDO cool phases. Due to the 'horseshoe' pattern of SST anomaly induced by PDO process, a basin-specific response of ocean desert to large-scale climate oscillator is proved here. However, the Chl-a concentrations show a positive correlation with the PDO processes in the southwest tongue area of NPOD. This is in line with NPOD SSH features, thus following the previous findings that regional SSHs in western tropical pacific negatively correlate with the PDO (Moon et al 2013, Hamlington et al 2014, Chang et al 2015 which further influence the Chl-a variation in the NPOD SSH .

Linkage of NPOD to vertical structure of ocean physics
The Chl-a variation in NPOD relies on the upperocean vertical mixing process, determined by the vertical density structure and vertical mixing efficiency, i.e. mixture coefficients (KM, referring to the stratification feature in KPP algorithm). Physically, the water density profiles depend on the coherent characteristics in temperature and salinity. In the NPOD SST region, we assess the coherence between the MLDs and the Chl-a variation of NPOD, by analyzing the change in monthly SST, vertical profiles of temperature, salinity and density, as well as the linkage with PDO index, as shown in figure 3. Under PDO positive conditions, e.g. during year 2004-2006 and 2014-2018, SSTs and temperatures of the upper ocean of 5-60 m show warm anomalies compared to an average throughout the entire interval from year 1998 to 2018. This suggests a downward propagation of the surface warm signal, controlled by the advection and diffusion processes within MLD (Sedláček andKnutti 2012, Llovel et al 2013). In contrast, such warmer upper ocean is accompanied by cooler deeper layer of 60-250 m, because the enhanced thermal MLD stratification inhibits further heat flux to deeper oceans. On one hand, vertical salinity structure presents 7 year variations, distinct from the density variability. Overall, the temporal variations in density profiles show resembled patterns as in temperatures, suggesting a dominant, thermal impact on the vertical stratification in the NPOD SST region (figures 3(b)-(d) and S8). thermal impact on the Moreover, the KMs in 5-30 m is counted with 1D KPP algorithm driven by the temperature, salinity, horizontal velocity fields in SODA3 reanalysis data. The KM and MLD coherently represent the capability of turbulent mixing in ocean column (figure 3(e)). The nutrients flux can be estimated by the diffusion equations with KM and climatological mean nutrients field in WOA (figure S9). For instance, under the warmer upper-ocean conditions of year 2004-2006 and 2014-2016, the nutrients flux into the surface layer is inhibited by the relaxation of vertical mixing represented by shallower MLD and smaller KM (figure S9) and finally induces the oligotrophication in the NPOD SST . For the NPOD SSH region, the correlation of 'Chl-a vs SSH' suggests that the dynamic condition determined by SSH variation dominates the depth of nutricline and then leads to the change of upward transport of nutrients, overcoming the impact of higher SSTs on Chl-a concentrations during year 1998-2018. However, our findings further imply that future rising in SSH under climate warming (Schewe et al 2011) may aggravate NPOD oligotrophication.
Here, we notice that other physical factors, e.g. changing ocean circulation (Andreas 2005, Yang et al 2020, nutrient advection flux (Dave andLozier 2015, Letscher et al 2016) or eddies (Liu andLevine 2016, Huang andXu 2018) may also contribute to the transport of nutrients to the surface layer. Previous studies have made attempt to figure out the biological impact of the climate warming on the Chl-a concentrations in the ocean (Behrenfeld et al 2016), leading to further complexity in ocean Chl-a variation, and more concrete evidence of NPOD variations and their causes will require a longer-term time series and more focused analyses.

Conclusions
Our findings show continuous oligotrophication and expansion of NPOD, as well as NW-SE variability of its mean positions but insignificant shift, during year 1998-2018. In our analysis, the oligotrophication occurred in most part of the NPOD region. The variations in NPOD intensity and area were found to moderately correlate to SST of surface seawater. The higher surface SST can strengthen vertical density stratification, subsequently resulting in lower nutrients flux into surface ocean (figure 4). In comparison, SW part of NPOD was analyzed to be linked with the SSH-controlled vertical nutrient flux into the upper ocean. In parallel, the NW-SE oscillation of the NPOD position on decadal time scales correlated to PDO processes.
Here, we notice that the existence of multiinterannual scale variability coincides with the decadal-scale strengthening process in NPOD during year 1998-2018 ( figure 1(a)). Therefore, in line with previous studies (Henson et al 2010, we encourage continuous observation and longer dataset for global ocean deserts. In particular, the oligotrophication and expansion processes have been also characterized in the SIOD and NAOD (figures S10 and S11), meanwhile the mechanisms remain elusive. Our study has focused on the linkage between NPOD and climate physics in ocean vertical structure, we believe more physical processes should be included and further modelling work need to be done to confirm and quantify the upper-ocean nutrient transport driven by the potential controlling factors mentioned in this study. Meanwhile the marine Chl-a variation is also linked to some biogeochemical-related processes e.g. phytoplankton metabolic rates, phytoplankton photoacclimation responses, dinitrogen fixations (Sommer and Sommer 2006, Taucher and Oschlies 2011, Sommer et al 2012, Olonscheck et al 2013, Lewandowska et al 2014, Behrenfeld et al 2016, Shiozaki et al 2018. To identify the biogeochemical impact on NPOD separately from the climate physics is thus suggested for future studies, e.g. using Earth system models.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).