1 Introduction

The North Pacific is characterized by areas of high-nutrient but moderate-chlorophyll because of iron limitation in the subarctic (Martine and Fitzwater 1988), and of low-nutrient and low-chlorophyll due to strong stratification in the sun-rich subtropics (Karl and Church 2019). The boundary region between the subarctic and the subtropics is an area of strong CO2 uptake (Takahashi et al. 2009).

A long-term decrease and ENSO-related variation of surface nutrients (phosphate, nitrate, and silicate) and biological production in the subarctic have been reported at several regions in the North Pacific (Freeland et al. 1997; Ono et al. 2002). A DIC increase in the western subtropics is also observed (Midorikawa et al. 2010; Ono et al. 2019) Goes et al. (2004) and Yasunaka et al. (2013, 2014a, b, 2016) made nutrient and DIC gridded products, and presented their basin-scale distributions and variabilities. However, these products use empirical relationships between nutrient concentrations and environmental variables, such as sea surface temperature (SST). Yasunaka et al. (2016) presented nutrient gridded products by an optimal interpolation, but their climatologies mostly depended on the empirical estimates of Yasunaka et al. (2014b). Assuming empirical relationships could induce systematic biases in some cases and may cause uncertainties (Wanninkhof et al. 2013; Yasunaka et al. 2019). Also, the marginal seas, such as the Okhotsk Sea and the Bering Sea, were often masked out probably because the same relationships as the open ocean could not be applied there.

Although spatial coverage of surface nutrient data has been improved thanks to ship-of-opportunity observations (Whitney 2011; Yasunaka et al. 2014b), uncertainty still exists in their basin-wide distribution and long-term variability (Yasunaka et al. 2016). DIC variability has been investigated often using DIC concentration calculated from partial pressure of CO2 (pCO2), because ocean surface pCO2 data have been dramatically increased owing to the underway measurements (Takahashi et al. 2009; Bakker et al. 2016). However, uncertainty from the DIC calculation is inevitable (Yasunaka et al. 2014a). Synthesized analysis of basin-scale nutrient and DIC variabilities in a consistent way, which has not been done yet as far as the authors know, will improve our understanding of seasonal and interannual variability.

Here, we make gridded products of phosphate, nitrate, silicate, and DIC by an interpolation technique in the North Pacific and its marginal seas, north of 10 °N from 1961 to 2016. We then assess their seasonal and interannual variability.

1.1 Data

Data used in this study are listed in Table S1.

1.2 Discrete water measurements for nutrients

Surface nutrient (phosphate, nitrate, and silicate) measurements obtained from water depths shallower than 10 m were extracted from the Global Ocean Data Analysis Project version2 (GLODAPv2; Olsen et al. 2016; https://www.glodap.info/) and the World Ocean Database 2018 (WOD 2018; Boyer et al. 2018; https://www.nodc.noaa.gov/OC5/WOD/pr_wod.html). We also used ship-of-opportunity (SOO) surface measurements of nutrients collected by the National Institute for Environmental Studies, Japan (NIES_SOO; Yasunaka et al. 2014b) and the Institute of Ocean Science, Canada (IOS_SOO; Whitney 2011). We also included nutrient concentrations measured by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) surface nutrient monitoring system (Yasunaka et al. 2016). Furthermore, we included nutrient concentrations from bottle samples collected by JAMSTEC (https://www.godac.jamstec.go.jp/darwin/datatree/e) and by IOS along Line P (Fisheries and Oceans Canada 2012; https://www.waterproperties.ca/linep/index.php) which have not been archived yet in GLODAPv2 and WOD2018. Duplicates with measurements archived in GLODAPv2 were excluded from WOD2018. The NIES_SOO, IOS_SOO, and JAMSTEC monitoring data are not included in WOD2018. The total number of measurements is about 21,000 for phosphate, 14,000 for nitrate, and 17,000 for silicate. Detection limits of NIES and IOS measurements are 0.02 μmol/l for phosphate, 0.1 μmol/l for nitrate, and 0.6 μmol/l for silicate (Yasunaka et al. 2014b).

The locations of the nutrient measurements are spread widely over the North Pacific, although many of them are concentrated around Japan, near the North American coast, and along several fixed lines repeatedly observed by research vessels and main routes of volunteer ships (Fig. 1a‒c). The numbers of phosphate and silicate measurements are similar throughout the analysis periods, but the number of nitrate measurements is low before the 1980s.

Fig. 1
figure 1

Total number of measurements in each 1° grid cell between 1961 and 2016 (left) and the total number of measurements each month in the North Pacific and its adjacent seas (10 °N–60 °N, 120 °E–90 °W) (right) for a phosphate, b nitrate, c silicate, and d pCO2

Full size image

1.3 Underway pCO2 measurements

We used surface water pCO2 measurements (converted from the fugacity of CO2 values; a correction of < 1%) from the Surface Ocean CO2 Atlas version 6 (SOCATv6; Bakker et al. 2016; https://www.socat.info/). The total number of measurements in the North Pacific (10 °N–60 °N, 120 °E–90 °W) from 1961 to 2016, is about 270,000 (Fig. 1d). Routine observations were started in the early 1970s along Line P (around 50 °N, 150 °W–125 °W; Wong et al. 2010), and then in the early 1980s along 137 °E (Inoue et al. 1995) and between North America and Australia (Wong et al. 1993). The number of the measurements increased in the mid-1990s when routine observations between Japan and North America were started (Murphy et al. 2001), and then in the mid-2000s, when the observations between Japan and Oceania region were started (Nakaoka et al. 2013).

1.4 Other gridded products

SST data were obtained from the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST) with a 1° × 1° monthly resolution (Rayner et al. 2003; https://www.metoffice.gov.uk/hadobs/hadisst/). Sea surface salinity (SSS) data were obtained from the quality controlled ocean temperature and salinity monthly objective analyses version 4 by the Met Office Hadley Centre (EN4; Good et al. 2013; https://www.metoffice.gov.uk/hadobs/en4/). Satellite chlorophyll-a concentration (Chl-a) with a 1° × 1° monthly resolution since September 1997 were downloaded from the website of the GlobColour project (GlobColour_R2018; https://hermes.acri.fr/index.php; Maritorena et al. 2010). Zonal mean data for the atmospheric CO2 mixing ratio (xCO2) were retrieved from the NOAA Greenhouse Gas Marine Boundary Layer Reference data product (NOAA_MBL; Conway et al. 1994; https://www.esrl.noaa.gov/gmd/ccgg/mbl/index.html) and were interpolated into 1° × 1° × 1 month grid-cells. Sea level pressure, 6-hourly 10-m wind speed data, and zonal wind data were obtained from the U.S. National Centers for Environmental Prediction–Department of Energy Reanalysis 2 (NCEP2) (Kanamitsu et al. 2002; https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.html). We used climatological means of mixed layer depths (MLD) produced by JAMSTEC (MILA_GPV; https://www.jamstec.go.jp/ARGO/srgo_web/argo/?page_id=223&lang=en/; Hosoda et al. 2010; missing data were interpolated with data in the surrounding grids). Surface current velocity data came from the monthly means satellite-derived Ocean Surface Current Analysis (OSCAR; Bonjean and Lagerloef 2002; https://www.oscar.noaa.gov/).

2 Methods

2.1 Gridded procedure

To derive gridded products of surface nutrients and pCO2, we performed statistical quality control, calculated long-term monthly means, and adapted optimal interpolation.

First, we eliminated measurements in each grid-cell that differed by more than three standard deviations from the long-term mean in ± 2° of latitude, ± 5° of longitude, and ± 1 month (regardless of the year). The window sides were set to include more than 100 measurements. These excluded data may have been erroneous, but would also include small-scale anomalies such as mesoscale eddies (Whitney and Robert 2002). This procedure identified in total about 0.5% of the data as outliers. For pCO2, we first normalized the data to the reference year 2005 based on changes in the atmospheric CO2 (Takahashi et al. 2002a, b; Nakaoka et al. 2013) to avoid the influence of the long-term increase of pCO2. Then, the remaining measurements were binned into 1° × 1° × 1 month grid-cells from 1961 to 2016.

We then calculated long-term monthly means and smoothed them by the exponential weighted average with space and time constants of 2° of latitude, 5° of longitude, and 1 month. We use these products (i.e. monthly climatologies) for detecting seasonal variation, and also as reference values for covariance analysis and optimal interpolation below.

Next, using anomalies of the quality controlled 1° × 1° × 1 month data from the monthly climatology, we calculated the data covariance at each 1° × 1° grid cell separated by a lag increment in the zonal, meridional, and temporal directions, and averaged them in the whole analysis area in each direction. Refer to Yasunaka et al. (2019) for detailed procedure. We defined the decorrelation radius as the space and time constants of the best fitted exponential decay line on the lag covariance in each direction. Since the difference in covariance between lag 0 and very short lag approximately represents the noise variance, the signal-to-noise ratio is obtained from the difference in covariance between lag 0 and the minimum lags (i.e. 1° in zonal and meridional directions here; the difference in covariance between lag 0 and a lag of 1 month cannot be considered as the noise variance since the decorrelation radius in a temporal direction is comparable). The resultant constant values of the best fitted exponential decay line on the lag covariance and differences in covariance between lag 0 and a lag of 1° are listed in Table S2. As a result, we set the decorrelation radius to 7° both in the zonal and meridional direction and 1 month in the temporal direction for nutrients, and 9° in the zonal direction, 5° in the meridional direction and 1 month in the temporal direction for pCO2. Signal-to-noise ratios were set to 1.7 for nutrients and 3.6 for pCO2. Since the lag covariance of nutrients is noisy in some cases, the same decorrelation radius and signal-to-noise ratio are used for phosphate, nitrate and silicate.

Finally, we applied optimal interpolation to the quality controlled 1° × 1° × 1 month data. Detailed equations of optimal interpolation presented by Hosoda et al. (2008). We used the monthly climatology as the first guess, and the decorrelation radius and signal-to-noise ratio obtained from the covariance analysis above. Window sizes for the interpolation were set to be twice the size of the decorrelation radius. For pCO2, we then added the atmospheric CO2 change initially subtracted (see above). As a result, we obtained phosphate, nitrate, silicate, and pCO2 values and their interpolation errors at each 1° × 1° × 1 month grid-cells in 10 °N–60 °N, 120 °E–90 °W from January 1961 to December 2016. In following analyses for interannual variability, we use the optimal interpolation results but only where the interpolation square error ratios are less than 0.9 to avoid reduction of the signals. The total number of data points where the error ratio is less than 0.9 during the analysis period is more than 100 north of 25 °N, but less than 100 in many parts south of 20°N and in the Bering and Okhotsk seas (Figure S1). That integrated over the North Pacific is generally more than 1000 in all analysis years for phosphate and silicate, after the late 1980s for nitrate, and after the mid 1990s for pCO2, but quite limited before the 1990s for pCO2.

2.2 Calculation of dissolved inorganic carbon

First, we used two estimates of total alkalinity (TA). One is from the Japan Meteorological Agency (TA_JMA; Takatani et al. 2014; https://www.data.jma.go.jp/gmd/kaiyou/english/co2_flux/co2_flux_data_en.html). The other was calculated from SSS and SST using the simple empirical function, TA = a + b (SSS − 35) + c (SSS − 35) 2 + d (SST − 20) + e (SST − 20)2, of Lee et al. (2006) (TA_Lee). The different coefficients “a” to “e” were specific to two regions in our analysis area: the subtropical region where SST is higher than 20 °C and the remaining North Pacific region. Since Takatani et al. (2017) pointed out a systematic bias in the North Pacific in Lee’s formulation, we used TA_JMA to examine seasonal variation. However, the Bering Sea, the Okhotsk Sea, and the Japan Sea are out of the coverage of TA_JMA, therefore, we used TA_Lee in those regions. For the interannual variation, we used TA_Lee because of their better coverage both in time and space. We confirmed the results are quite similar with both TA in terms of the interannual variation. Uncertainty of the TA estimate is 8 μmol/kg (Lee et al. 2006; Takatani et al. 2014).

From pCO2 gridded products and total alkalinity (TA) estimates, we calculated sea surface DIC concentration, using the CO2SYS program (Lewis and Wallace 1998; van Heuven et al. 2009) with the dissociation constants of Lueker et al. (2000). Uncertainties of the calculation (0.5%; Lueker et al. 2000) and TA estimation (8 μmol/kg ≈ 0.4%) derive uncertainty of DIC estimates of 12 μmol/l (sqrt (0.5%2 + 0.4%2) = 0.6%).

We also calculated salinity normalized DIC to a constant salinity of 35 (nDIC = 35 × DIC/SSS).

2.3 Calculation of seasonal drawdown

Following Yasunaka et al. (2014b), we calculated carbon drawdown from month m1 to month m2 using concentration changes, subtracting the effects of dilution/condensation, horizontal advection, and air-sea flux:

$${\text{C}}\_{\text{m1 }}{-}{\text{ C}}\_{\text{m2 }} = {\text{ nDIC}}\_{\text{m1}}{-}{\text{nDIC}}\_{\text{m2 }} + {\text{ }}({\text{u}}\_{\text{m1}}\cdot\nabla {\text{DIC}}\_{\text{m1 }} + {\text{ u}}\_{\text{m2}}\cdot\nabla {\text{DIC}}\_{\text{m2}}){\text{ }}/{\text{ 2 }} + {\text{ }}\left\{ {\left( {{\text{Cflux}}\_{\text{m1 }} + {\text{ Cflux}}\_{\text{m2}}} \right){\text{ }}/{\text{ 2}}} \right\}{\text{ }}/{\text{ MLD}}\_{\text{m2}},$$

where u is the horizontal velocity, ∇DIC is the horizontal gradient of DIC concentration, and Cflux is air-sea CO2 flux (positive value indicates efflux from sea to air). For the calculation of air-sea CO2 flux, we used the gas exchange coefficient of Sweeney et al. (2007) but optimized for NCEP2 winds following Takahashi et al. (2009), and the Schmidt number of CO2 in seawater at a given SST based on Wanninkhof (2014). We also assumed that vertical diffusion was negligible during the mixed layer shoaling period (Lee 2001).

In the similar way, we calculated nutrient as follows:

$${\text{Nut}}\_{\text{m1 }}{-}{\text{ Nut}}\_{\text{m2 }} = {\text{ nNut}}\_{\text{m1 }}{-}{\text{ nNut}}\_{\text{m2 }} + {\text{ }}({\text{u}}\_{\text{m1}}\cdot\nabla {\text{Nut}}\_{\text{m1 }} + {\text{ u}}\_{\text{m2}}\cdot\nabla {\text{Nut}}\_{\text{m2}}){\text{ }}/{\text{ 2}},$$

where nNut salinity-normalized nutrient concentrations (to S = 35), and ∇Nut is the horizontal gradient of nutrient concentration. Atmospheric input of nutrients was not assessed, since its contribution was estimated to be smaller than 10% in the open ocean (Krishnamurthy et al. 2010).

2.4 Decomposition of interannual variation

To detect nutrient variations associated with the Pacific Decadal Oscillation (PDO), the North Pacific Gyre Oscillation (NPGO), both of which have been known as the dominant large-scale climate variabilities in the North Pacific (Mantua et al. 1997; Di Lorenzo et al. 2008), and the linear trend, we multiply regressed the PDO index, the NPGO index, and the linear trend time series on the nutrient anomalies in each grid as in Yasunaka et al. (2016). A multiple linear regression equation with three variables is represented as follows:

$${\text{Y}}\, = \,{\text{b1}}\cdot{\text{X1}}\, + \,{\text{b2}}\cdot{\text{X2}}\, + \,{\text{b3}}\cdot{\text{X3}}\, + \,\varepsilon ,$$

where Y is the nutrient or DIC concentration anomaly from long-term monthly mean, X1, X2, and X3 are the PDO index, NPGO index, and linear trend time series respectively. Values of b1, b2, and b3 (the partial regression coefficients) are found by minimizing the sum of the squares of the ε values.

3 Results and discussion

3.1 Seasonal variation

3.1.1 Climatological mean state

Surface nutrient and DIC reached maximum concentrations in February or March, and minimum concentrations in August or September in almost all of the North Pacific, north of 30 °N (Figure S2a‒d), as Yasunaka et al. (2013; 2014b) showed for area averaged nutrients. Nutrient and DIC concentrations were high in the subarctic especially in winter (> 1 μmol/l for phosphate, > 10 μmol/l for nitrate, > 20 μmol/l for silicate, and > 2050 μmol/l for DIC) and low in the subtropics (< 0.5 μmol/l for phosphate, < 5 μmol/l for nitrate, < 10 μmol/l for silicate, and < 2050 μmol/l for DIC; Fig. 2a‒d). Standard error of the climatological means is much less than the mean concentration, although it is sometimes large in the Okhotsk and Bering seas due to the large small-scale variance and limited number of the observations (Figure S3). Compared to our previous studies (Yasunaka et al. 2013, 2014b), noisy features were suppressed mainly because more data were used.

Fig. 2
figure 2

a Phosphate, b nitrate, c silicate, and d DIC concentrations in March (left) and August (right)

Full size image

Marine biological production is often limited by the supply of nutrients. The limitation occurs not abruptly at a certain concentration, but gradually with decreasing of nutrient concentration. The point at which limitation starts varies across the North Pacific due to environmental condition such as sunlight intensity or iron supply. A half-saturation point in the Michaells–Menten equation is used as an indicator of the limitation criteria especially in the numerical models. Here, we show the lines of 0.2 μmol/l phosphate, 2 μmol/l nitrate, and 4 μmol/l silicate in Fig. 2 based on a synthesized study of diatom limitation by Sarthou et al. (2005). In the central region, the northern edge of the limitation is around 33°N in winter and 40°N in summer for all nutrients (Fig. 2a‒c). Phosphate concentration is higher than the criteria in the eastern subtropics, and silicate concentration is higher than the criteria in the western subtropics in both seasons. Those results mean nutrients are sufficient for the growth of diatoms north of 40°N through out the year, but are limiting growth south of 40°N in summer.

N* (= [nitrate] − 16 × [phosphate]) is an indicator of deviations from typical oceanic nutrient ratios (N:P ~ 16:1) caused by nitrogen fixation, denitrification, and atmospheric nitrogen deposition (Gruber and Sarmiento 1997). Although a broad perspective can be derived from existing datasets such as the World Ocean Atlas (Garcia et al. 2018), our map provides better regional detail (Fig. 3). N* is relatively high in the western and central North Pacific south of 40°N (> − 3), and low in the subarctic and the eastern subtropics (< − 5). High N* in the subtropical gyre can be accounted for by nitrogen fixation, as shown by a numerical model (Deutsch et al. 2007) and by direct observations (Shiozaki et al. 2010). N* in the subarctic and east of Japan is quantitatively consistent with Yoshikawa et al. (2018); − 8 to − 5 along 49°N, and − 3 to − 1 east of Japan. In the Bering Sea, N* is low in the eastern shelf region (< − 7), and high in the western basin area (> − 7). These values are consistent with Granger et al. (2011). Low N* in the eastern Bering Sea shelf is caused by intense denitrification (Koike and Hattori 1979; Devol et al. 1997; Tanaka et al. 2004). In the Okhotsk Sea, N* is higher to the north than the south. This is consistent with Yoshikawa et al. (2006), who attributed the higher N* to atmospheric nitrogen deposition and river discharge. N* in the Bering Sea basin and around the Aleutian Islands shows large seasonal variation, high in winter and low in summer (not shown here). Vertical profiles of N* in this area in summer (Lehmann et al. 2005) showed lower values in surface waters than were found in the subsurface. Since subsurface water in summer is the water in the mixed layer in winter, our seasonal variation of N* is consistent with the N* profiles.

Fig. 3
figure 3

Annual mean of N* (= nitrate − 16 × phosphate)

Full size image

3.1.2 Amplitude and phase

Standard deviations of seasonal variations of nutrients and DIC are large in the northwestern part of the subarctic including the Okhotsk and Bering seas (> 0.2 μmol/l for phosphate, > 3 μmol/l for nitrate, > 4 μmol/l for silicate, and > 30 μmol/l for DIC; Figure S4a‒d). They also show local maxima in the transition zone between the subarctic and subtropics. Standard deviations of seasonal nutrient variations are very small in the subtropics (< 0.1 μmol/l for phosphate, < 1 μmol/l for nitrate, and < 2 μmol/l for silicate), while DIC shows substantial variation in northwestern and the southeastern subtropics (10‒30 μmol/l), half of which is accounted for by the salinity variation (Figure S4e). In the North Pacific, the seasonal variation of SSS is largest in northwestern and southeastern subtropics due to large seasonal variation of evaporation and precipitation, respectively (Bingham et al. 2010).

Nutrient and DIC seasonal variability north of 30°N shows positive correlation with MLD seasonal patterns, and negative correlation with SST (Fig. 4). MLD north of 30°N reaches its maximum in February or March and minimum in July or August, and SST is maximum in August and minimum in March (Figure S2e‒f and S5a‒b). The seasonal cycle of nutrient and DIC concentrations results from entrainment of nutrient-rich subsurface water in late winter and biological consumption from spring to summer. Although the seasonal variation of SST also affects DIC levels due to temperature-dependency of pCO2, the effect north of 30°N is quite limited (not shown here).

Fig. 4
figure 4

Correlation of seasonal variation with MLD (left), SST (middle) and Chl-a (right) for a phosphate, b nitrate, c silicate, and d DIC. Shaded areas indicate correlation coefficients significant at p < 0.05. Values are omitted in grids where the minimum concentration is smaller than 0.02 μmol/l of phosphate, 0.1 μmol/l of nitrate, and 0.6 μmol/l of silicate, and standard deviation of DIC seasonal variation is smaller than 8 (= 16/sqrt(12/3)) μmol/l

Full size image

The correlation between nutrient and DIC seasonality and Chl-a variation is positive between 30°N and 40°N, and negative north of 40°N (Fig. 4). Between 30°N and 40°N, Chl-a is maximum in March or April, and minimum in July or August (Figure S2g and S5c). In this area where sunlight is abundant throughout the year, biological production is controlled by nutrient supply (Polovina et al. 1995). North of 40°N, away from the coast, Chl-a is maximum from June to September, and minimum in March or April. Biological production starts to be active when the ocean is warm and sunlight is sufficient, resulting in nutrient consumption. High levels of Chl-a continue until summer and autumn in this region as Goes et al. (2004) and Sasaoka et al. (2011) showed. The boundary between positive and negative correlations of nutrients and Chl-a seasonal variations is the same latitude (40°N) where the nutrients limit diatom growth in summer (see Sect. 4.1.1). In the coastal area, Chl-a is maximum in May, which is not precisely opposite to nutrient and DIC variations.

3.1.3 Seasonal drawdown

Seasonal drawdown from March to August of nutrients and DIC is large in the northwestern part of the subarctic (> 0.4 μmol/l for phosphate, > 12 μmol/l for nitrate, > 15 μmol/l for silicate, and > 100 μmol/l for DIC; Fig. 5a‒d), and it is especially large in the Okhotsk and Bering seas and around the Kuril Islands (> 1.0 μmol/l for phosphate, > 15 μmol/l for nitrate, > 25 μmol/l for silicate, and > 125 μmol/l for DIC). Local maxima of the seasonal drawdown in the north–south direction can be seen around 40°N, which is most prominent in DIC. Net primary production in the mixed layer (drawdown integrated over MLD) from March to July is more than 6 gC/m2/mo, i.e. 24 gC/m2/4mo, in the boundary region between the subarctic and subtropics, the western subarctic, the Gulf of Alaska, and the Bering Sea (Fig. 6). Our estimates are consistent with other regional estimates; 9 gC/m2/mo in Oyashio-region (38 °N–43 °N, 141 °E–148 °E; Ono et al. 2005), and 5 gC/m2/mo at station P (50 °N, 145 °W; Lockwood et al. 2012 and references therein). Compared to our previous studies (Yasunaka et al. 2014b), the inclusion of more data and the suppression of noisy features enable us to discuss the stoichiometric relationship of seasonal drawdown in more detail as follows.

Fig. 5
figure 5

Seasonal drawdown from March to August for salinity-normalized a phosphate, b nitrate, c silicate, and d DIC, and their stoichiometric ratios of seasonal drawdown of e phosphate, f silicate, and g DIC to nitrate. Effect of air-sea CO2 flux was eliminated. Values are omitted in grids where concentration in March or August is smaller than 0.1 μmol/l of nitrate in (e–g)

Full size image
Fig. 6
figure 6

Net community production in mixed layer from March to July. Colors are omitted in grids where MLD does not monotonically decrease from March to July

Full size image

A relatively large seasonal nitrate-to-phosphate drawdown ratio (ΔN: ΔP > 18) is found in the central subarctic and the southern Bering Sea (Fig. 5e). Nitrate-deficient water from the Bering Sea shelf region spreads to the surface layer of its basin (Lehmann et al. 2005). In the Bering Sea basin, nitrate-deficient water from the shelf would remain at surface in summer, whereas it mixes with nitrate rich subsurface water in winter. Seasonal drawdown of nitrate to phosphate in the other region is similar to the Redfield ratio (16:1).

Drawdown of silicate is more than 1.5 times that of nitrate in the Okhotsk Sea, the Bering Sea, and the western subarctic (Fig. 5f). Drawdown of silicate is smaller than 1.5 times of nitrate drawdown in the coastal region of the Alaskan gyre and along 40°N. Relative large drawdown of silicate in the Bering Sea and the western subarctic, and small silicate drawdown in the eastern subarctic are also pointed out by Koike et al. (2001) using the profiling data in summer at several points. Since diatoms assimilate silicate to grow, large silicate drawdown relative to nitrate indicates rich abundance of diatom. Although the Si:N ratio of diatom’s chemical composition is normally 1.12 ± 0.33 (Brzezinski 1985), that under stress condition, which often appear in the North Pacific, is higher by a factor of about 2–3 (Harrison et al. 1976, 1977; Takeda 1998; Takeda et al. 2005). The spatial pattern of drawdown ratio of silicate to nitrate is quite similar to diatom distribution estimated by Hirata et al. (2011); rich abundance of diatoms in the northwestern North Pacific, the Okhotsk Sea, and the Bering Sea. East–west contrast of phytoplankton community composition in the North Pacific was known from stationary observations (Harrison et al. 2004, and their references). Diatoms account for most of the large production in the Okhotsk Sea (Nakatsuka et al. 2004), and are dominant in the spring bloom in the Bering Sea (Takahashi et al. 2002a, b).

Seasonal drawdown of carbon to nitrate is larger than the Redfield ratio (106:16) in most regions, and becomes large to the south (Fig. 5g). The C:N drawdown ratio in excess of the Redfield ratio is called as “carbon overconsumption” (Toggweiler 1993). Biological uptake in warm regions results in a C:N ratio of > 30 in some cases (Taucher et al. 2012). The C:N drawdown ratio is relatively low in the western subarctic and the Bering Sea. The C:N drawdown ratio tends to be smaller in iron-rich waters (Takeda 2011), and the iron supply to surface water is relatively abundant in the western subarctic and the Bering Sea (Uchida et al. 2013; Nishioka and Obata 2017). Wong et al. (2002) attributed the variation of the C:N drawdown ratio in the subarctic North Pacific to calcium carbonate formation. Coccolithophores can be a dominant phytoplankton group in the eastern subarctic, whereas diatoms are dominant in the western subarctic and in the Bering Sea (Harrison et al. 2004, and their references).

3.2 Interannual variation

3.2.1 Amplitude and phase relationship with SST and Chl-a

Interannual variability of surface phosphate, nitrate and silicate concentration are large north of 40 °N, as much as 50% of the seasonal change (standard deviations are > 0.1 μmol/l for phosphate, > 1 μmol/l for nitrate, and > 2 μmol/l for silicate; Figure S6a‒c). Interannual variation of DIC has a standard deviation of about 10‒15 μmol/l in most of the area, which is 50% of the seasonal variation north of 40 °N and larger than the seasonal variation south of 30 °N (Figure S6d). Relative amplitude difference between seasonal and interannual variations shows the same tendency when the unbiased estimate of population variance, or difference between 10 and 90 percentiles is used.

Interannual variation of nutrient and DIC concentrations shows negative correlation with SST interannual variation north of 30 °N (Fig. 7). These relationships are the same as in seasonal variation. Nutrient and DIC concentrations are high when SST is low. The low SST relates to southward Ekman advection, entrainment of subsurface water, and surface heat flux (Mochizuki and Kida 2003), first two of which have an effect on increasing of nutrient and DIC concentrations. As with seasonal variation, the correlation between nutrients and Chl-a is positive south of 40 °N and east of 160 °E and negative north of 40 °N, although the significant area is limited. Biological production is controlled by nutrient availability in the subtropics, whereas nutrients remain abundant in the subarctic throughout the year. Although the noisy features in correlation maps with Chl-a might be caused by the non-linear and/or lag relationship, data are not enough to examine it and these topics are left for future studies.

Fig. 7
figure 7

Correlation of interannual variation with SST (left) and Chl-a (right) for a phosphate, b nitrate, c silicate, and d DIC. Shaded areas indicate correlation coefficients significant at p < 0.05. Values are omitted in grids where the minimum concentration is smaller than 0.02 μmol/l of phosphate, 0.1 μmol/l of nitrate, and 0.6 μmol/l of silicate, and standard deviation of DIC interannual variation is smaller than the error (i.e. 16/sqrt(data number/3) μmol/l)

Full size image

3.2.2 Variation related to PDO and NPGO

The PDO and the NPGO affect nutrient concentrations as shown by Yasunaka et al. (2016). When the PDO is in its positive phase, nutrient concentrations in the western North Pacific are significantly higher than the climatological means, and those in the eastern North Pacific are significantly lower (Fig. 8). When the NPGO is in its positive phase, nutrient concentrations in the subarctic are significantly higher than the climatological means. Nutrient signals in the subtropics are not affected, while DIC does vary. Changes in horizontal advection and vertical mixing associated with the PDO and the NPGO are correlated to changes in nutrients; mixed layer deepening induces an increase of nutrient concentrations via enhanced entrainment of subsurface, nutrient-rich water, and stronger westerly winds force more nutrient-rich water southward from higher latitudes (Yasunaka et al. 2016). A PDO-related DIC signal in the subtropics is associated with the dilution/condensation effect by the flesh water flux variation (Yasunaka et al. 2014a). Heavy rainfall around the Hawaiian Islands and higher salinity surface water in the western subtropics have been observed when the PDO is in the positive phase (Minobe and Nakanowatari 2002; Chu and Chen 2005; Delcroix et al. 2007), which would induce negative and positive DIC anomalies, respectively. North–south salinity contrast in the eastern Pacific associated with NPGO was shown by De Lorenzo et al. (2009). PDO- and NPGO-related nDIC signals in the subtropics are weak (not shown here).

Fig. 8
figure 8

Multiple regression maps onto the PDO (left), the NPGO (middle) and liner trend (right) for a phosphate, b nitrate, c silicate, and d DIC. Right panels show linear trends per decade for a phosphate, b nitrate, and c silicate, and that per year for d DIC. Shaded areas indicate regression coefficients significant at p < 0.05. Values are omitted in grids where the minimum concentration is smaller than 0.02 μmol/l of phosphate, 0.1 μmol/l of nitrate, and 0.6 μmol/l of silicate, and standard deviation of DIC interannual variation is smaller than the error (i.e. 16/sqrt(data number/3) μmol/l)

Full size image

3.2.3 Long-term trend

The nutrient change associated with the long-term trend averaged over the analysis area is − 0.010 ± 0.002 μmol/l/decade for phosphate, − 0.01 ± 0.11 μmol/l/year for nitrate, and − 0.64 ± 0.05 μmol/l/decade for silicate (average ± standard error), which are consistent with Yasunaka et al. (2016). As a result of warming and freshening, surface ocean density is decreasing in the large area of the North Pacific (Yasunaka et al. 2016). The decrease of surface ocean density suggests an increase of the upper ocean stratification, and therefore, a declining supply of nutrient to surface waters. In addition, anthropogenic nitrogen deposition increases on surface nitrate supply as discussed in Yasunaka et al. (2016). However, spatial patterns of nutrient long-term trends are not consistent with each other (Fig. 8a–c). Longer time-series are still needed to clarify the spatial patterns of nutrient trends and their interpretation.

The DIC change associated with the long-term trend is positive all over the analysis area (Fig. 8d). The DIC trend averaged over the analysis area is 0.77 ± 0.03 μmol/l/year (0.75 ± 0.02 μmol/kg/year), but shows spatial variation. It is larger than 1.0 μmol/l/year in the western and central subtropics. The nDIC trend along 137°E is quite similar to that reported by Ono et al. (2019) (0.88–1.10 μmol/kg/year in this study, and 0.89‒1.20 μmol/kg/year in Ono et al.).

Here, we calculated the Revelle factor (R), which is the ratio of instantaneous change in pCO2 to the change in DIC (Revelle and Suess, 1957);

$${\text{R}}\, = \,\left( {\Delta p{\text{CO}}_{{\text{2}}} /p{\text{CO}}_{{\text{2}}} } \right){\text{ }}/{\text{ }}\left( {\Delta \left[ {{\text{DIC}}} \right]{\text{ }}/{\text{ }}\left[ {{\text{DIC}}} \right]} \right).$$

Δ[DIC] at ΔpCO2 = 1 is calculated using CO2SYS with the assumption of Δ[TA] = 0 (Broecker et al. 1979). Climatological means are used for pCO2, [DIC] and other parameters. High value of the Revelle factor indicates a low capacity of the ocean to absorb atmospheric CO2. The Revelle factor is low in the western and central subtropics and higher to the north (Fig. 9) as Sabine et al. (2004) showed. It generally corresponds in the opposite sign with the spatial contrast of the long-term trend of DIC (Fig. 8d). The long-term trend of nDIC is also large in the subtropics and the smaller to the north, although it is relatively small in the southwest subtropics as Ono et al. (2019) discussed (Figure S7).

Fig. 9
figure 9

Revelle factor (= (ΔpCO2 / pCO2) / (Δ[DIC] / [DIC])

Full size image

4 Concluding remarks

In this study, we made gridded products of surface water phosphate, nitrate, silicate, and DIC using a technique to interpolate measured values across the North Pacific and its marginal seas, then assessed their seasonal to interannual variability. We characterized these variabilities and quantified the seasonal drawdowns and stoichiometric ratios. Compared to our previous studies (Yasunaka et al. 2013, 2014a, b, 2016), more data were used, and analysis periods were extended. They enabled us to present the monthly maps without assuming empirical relationships with environmental variables. As a result, noisy features were suppressed especially in the monthly climatologies. Seasonal drawdowns and their stoichiometric relationships were discussed in more detail, and the spatial pattern of the DIC trend was presented. The gridded products presented here are available at https://www.jamstec.go.jp/res/ress/yasunaka/nutrient/.

Increasing the number of measurements is crucial to clarify the variabilities, especially in the large spatial and long temporal scales. Continuing SOO measurements of nutrients and pCO2 is desirable. Installing nitrate sensors as part of underway observing systems and increasing the coverage of biogeochemical floats are promising ways to increase surface nutrient data. Most surface nutrient concentrations in the subtropics are below the limit of detection of conventional measurements, but recent high-sensitivity analyses have revealed non-zero concentrations in surface waters (Hashihama et al. 2013). Accumulation of such high-quality data is needed to clarify the nutrient variability throughout the subtropical Pacific Ocean. Data synthesis is as important as the observation itself. Data archiving and quality control efforts such as SOCAT, GLODAP, and WOD should continue to be encouraged.