Climatic Change(2021) 168:29
https://doi.org/10.1007/s10584-021-03239-1

Long‑term variations in ocean acidification indicesin the Northwest Pacific from 1993 to 2018
Miho Ishizu1,2· Yasumasa Miyazawa3 · Xinyu Guo4

Received: 27 December 2020 / Accepted: 2 October 2021
© The Author(s) 2021

AbstractLong-term variations in ocean acidification indices in the Northwest Pacific were examinedusing observational data and a biogeochemical model with an operational ocean model
product for the period 1993–2018. The model and observational data for the surface ocean(< 100-m depth) exhibit consistent patterns of ocean acidification in the subtropical andKuroshio Extension regions and relative alkalinization (i.e., reduced acidification) in thesubarctic region of the Northwest Pacific. Below 100-m depth, acidification dominated in
the subtropical regions and alkalinization in the subarctic regions. We attribute the excess
acidification in the subtropical and Kuroshio regions to the vertical mixing of dissolved
inorganic carbon (DIC) exceeding the DIC release by air–sea exchange. These regional
differences in acidification and alkalinization are attributed to spatially variable biologicalprocesses in the upper ocean and horizontal and vertical physical redistribution of DIC.
Our model and observational results have implications for the spatial extent and pattern ofocean acidification, along with the strength of the ocean carbon sink, which are key aspectsof global climate change.
Keywords Biogeochemical model · Northwest Pacific · Ocean acidification · Ocean
alkalinization · JCOPE* Miho Ishizu
mishizu@pusan.ac.kr
Yasumasa Miyazawa
miyazawa@jamstec.go.jp
Xinyu Guoguoxinyu@sci.ehime-u.ac.jp
1

Center for Climate Physics, Institute for Basic Science, Busan 46241, Republic of Korea

2

Pusan National University, Tonghapgigyegwan Bldg 2 Busandaehak‑ro, 63 beon‑gil,Geumjeong‑gu, Busan 46241, Republic of Korea

3

The Japan Agency for Marine‑Earth Science and Technology, Environmental Variability
Prediction and Application Research Group, Yokohama Institute for Earth Sciences, 3173‑25
Showa‑machi, Kanagawa‑ku, Yokohama 236‑0001, Japan

4

Center for Marine Environmental Studies, Ehime University, 2‑5 Bunkyo‑cho,Matsuyama 790‑8577, Japan

13

Vol.:(0123456789)

29

Page 2 of 20

Climatic Change(2021) 168:29

1 Introduction
The atmospheric partial pressure of ­CO2 ­(pCO2, air) has been increasing at a rate of ~ 1.8
parts per million by volume (ppmv) per year in recent decades, which is a consequence ofhuman activities such as fossil fuel burning, deforestation, and cement production (Takahashi et al. 2009; IPCC 2013). The contemporary global ocean now acts as a net sink of­ O2 outgasanthropogenic ­CO2, but this net flux reflects an imbalance between regions of C
sing and uptake. In general, the subtropics act as a source for C­ O2 and the subarctic is asink (Takahashi et al. 2009), but there is strong regional and temporal heterogeneity in theair–sea ­CO2 flux, because the sea surface ­CO2 concentration is affected by ocean circulation (subduction and obduction; Toyama et al. 2017; Ono et al. 2019; Ishii et al. 2020)and biological activity (Takahashi et al. 2002; Yasunaka et al. 2013). As such, the sea surface ­CO2 concentration can become higher than that of the atmospheric ­CO2 concentration, and vice versa depending on the regional internal oceanic processes. Interannual anddecadal variations in the climate system can also alter the patterns of mean-state ­CO2 flux
between ocean and atmosphere, and on longer time scales, anthropogenic perturbations to
the ocean’s physical and chemical state may adjust the spatial patterns of C
­ O2 flux.
The uptake of anthropogenic ­CO2 by seawater changes the natural chemical equilibrium of the seawater C­ O2–carbonate system, resulting in ocean acidification (Bates et al.
2014). In the Northwest Pacific, Midorikawa et al. (2010) estimated surface p­ H25 trends
from a 25-year time-series along the 137°E line south of 35°N using observational datafrom the Japan Meteorological Agency (JMA). The results revealed that ­pH25 decreased
by 0.0013 ± 0.0005 ­yr–1 in summer and 0.0018 ± 0.0002 ­yr–1 in winter from 1983 to 2007.
Wakita et al. (2013) investigated variations in ocean acidification during winter in the subarctic Northwest Pacific from 1997 to 2011 (K2 = 47°N, 160°E; KNOT = 44°N, 155°E) andfound that ­pHin situ in the winter mixed layer decreased at a rate of − 0.0011 ± 0.0004 ­yr–1,
which was lower than the average rate in the subtropical region (− 0.002 ­yr–1; Midorikawa
et al. 2010). Their results revealed a decline in p­ Hin situ at a rate of − 0.0051 ± 0.0010 ­yr–1
between 200 and 250 m below the mixed layer, which corresponds to the density zone of26.9σθ kg ­m–3, and is greater than any previously reported ­pHin situ decrease in the surface
layer in the open North Pacific.
Monitoring for biogeochemical variables including ocean acidification is importantfor understanding future ocean chemistry and the health of marine organisms and ocean
ecosystems. Chemical variables related to the oceanic carbonate cycle are currently being
monitored in several ongoing projects to determine whether the rate of ocean acidification
can be identified from pH changes and other variables in the open ocean (Bates 2007; Gonzalez-Davila et al. 2007; Dore et al. 2009; Olafsson et al. 2009; Midorikawa et al. 2010;
Bates et al. 2014; Carter et al. 2017; Chen et al. 2017; Wakita et al. 2013). However, making synoptic- to basin-scale inferences from such ocean acidification studies is challenging
due to the sparsity of p­ Hin situ measurements ­(pH25) in time and space.
Our current understanding of ocean acidification at present is based on time-series pH
measurements at several fixed monitoring stations (Bates 2007; Gonzalez-Davila et al.
2007; Dore et al. 2009; Olafsson et al. 2009; Midorikawa et al. 2010; Bates et al. 2014;
Wakita et al. 2013). As such, a robust three-dimensional overview of basin-scale ocean
acidification is currently unavailable.
The Japan Coastal Ocean Predictability Experiment (JCOPE; http://​www.​jamst​ec.​go.​jp/​jcope/) model was developed to integrate the sparse observational constraints and simulations of ocean physics. JCOPE is an operational eddy-resolving ocean physical model

13

Climatic Change(2021) 168:29

Page 3 of 20

29(Miyazawa et al. 2017, 2019) that incorporates observational data from floats, satellites,and shipboard measurements. Several studies have successfully used the JCOPE model
product to study ocean dynamics, including 10-yr oceanic variations (Miyazawa et al.
2009, 2017; Saeyanto et al. 2014; Chang et al. 2015; Miyama et al. 2021). Using the
JCOPE reanalysis model products as the physical background, Ishizu et al. (2019, 2020)
recently developed an offline biogeochemical model (JCOPE_EC). The model includes the
carbon cycle and state variables related to ocean acidification indices. The results showedthat the model reliably reproduces the distributions of the biogeochemical state variables(DIN, DIC, and ALK) on a seasonal scale (Ishizu et al. 2019, 2020).
In this study, we examined the recent trends of long-term (i.e., multi-decadal) variationsin the inorganic marine carbon cycle, using JCOPE_EC. In Sect. 2, details of the modeland data are described. Section 3 presents the results of the multi-decadal simulations of
the ocean acidification indices, and a comparison between the available observational data
and model outputs. In Sect. 4, air–sea C­ O2 fluxes and potential mechanisms of the multidecadal variations in ocean acidification are discussed. Finally, the conclusions are presented in Sect. 5.

2 Model and data
JCOPE_EC (Ishizu et al. 2019, 2020) is an off-line tracer model driven by physical processes simulated by an operational eddy-resolving ocean general circulation model(JCOPE2M; Miyazawa et al. 2017, 2019), which is based on the Princeton Ocean Modelwith a generalized sigma coordinate (Mellor 2001). The model is a three-dimensional
regional model covering the Northwest Pacific (10.5°–62°N, 108°–180°E) with a horizontal resolution of 1/12 of a degree (4.4–9.1 km) and 46 active vertical levels. The set-up of
the carbon cycle in the model is the same as the general form of the OCMIP protocol (Orret al. 1999). The JCOPE_EC model calculates ­pCO2 based on temperature, salinity, DIC,
and ALK (Ishizu et al. 2019, 2020). The atmospheric ­CO2 values are obtained from anempirical formula based on observational data (Ishizu et al. 2019). The reproducibility ofthe model output above 500-m depth has been demonstrated for the seasonal carbon cycle(Ishizu et al. 2020).
The model was driven by forcing from daily oceanic (JCOPE2M) and six-hourly atmospheric (NCEP/NCAR) reanalysis data from 1993 to 2018, yielding an analysis period of
26 years. The initial concentrations of phytoplankton were set to 0.1 and 0.0 mmol ­m–3
for depths above and below 150 m, respectively. The initial zooplankton concentrations
were set to 10% of the phytoplankton concentrations. The initial detritus concentration was
set to 0.0 mmol ­m–3. Dissolved inorganic nitrogen (DIN) and phosphate (DIP), and DIC
were initialized using the monthly climatological data for January, and ALK was initializedusing the annual climatology at the beginning of the model run (Ishizu et al. 2019, 2020).
The boundary conditions were set using the monthly climatology for DIN, DIP, and DIC,
and the annual climatology for ALK (Ishizu et al. 2019, 2020). JCOPE_EC is described in
more detail in Ishizu et al. (2019, 2020).
Temperature, DIC, and ALK drive the annual variations in ocean acidification indices, but DIC exerts the strongest influence. We used observational profiles to confirm thatchanges in ­pHin situ ­(pH25 and Ωarg) are primarily linked to changes in DIC (not shown).
To validate the reproducibility of the long-term variations in ocean acidification from themodel, we compared observations of DIC from 1999 to 2018 with the model outputs along

13

29

Page 4 of 20

Climatic Change(2021) 168:29

the same section (line 165°E) in Sect. 3.2. There are no available observations of ALK thatencompass at least a 10-year interval within the target period from 1999 to 2018.
The observational data for 2017 and 2018 along line 165°E were obtained from theJapan Meteorological Agency (JMA; https://​www.​data.​jma.​go.​jp/​kaiyou/​db/​vessel_​obs/​
data-​report/​html/​ship/​ship.​php), and the data for 1992, 2000, and 2006 were obtained from
PACIFic ocean Interior Carbon (PACIFICA; https://​clima​tedat​aguide.​ucar.​edu/​clima​te-​data/​pacif​i ca-​pacif​i c-​ocean-​inter​ior-​carbon) dataset. Monitoring data from line 165°E in
the Northwest Pacific were used for the analysis described in Sect. 3.2.
The temporal variations in DIC calculated from the model were relatively large between
1993 and 1998 (Appendix), but became stable after 1999. Therefore, we regard 1993–1998
to be a spin-up period, and only use the model outputs after 1999 for the analysis (Figs. 2,
3 and 6, 7, 8, 9, 10, 11; Table 1). We use the following equation to quantify the rate ofchange in state variables C (x, y, z, and t) from 1999 to 2018 in Sect. 3.1:
C (x, y, z, 2018 − 1999) = (< C(x, y, z, 2018) > − < C (x, y, z, 1999) >)∕(2018 − 1999 + 1)

(1)
where <  > denotes the annual mean of the state variable for the specified year. We also calculated linear trends for DIC, p­ Hin situ, ­pH25, and Ωarg at each depth in each region (Table 1;
Appendix Fig. 13), using the annual mean values for every year, which is different from thecalculation method of Eq. (1).

Table 1  Averaged trends of model outputs in DIC, ­pHinsitu, ­pH25 and Ωarg at each depth in the subarctic
region (155°–175°E, 40°–50°N), Kuroshio Extension (142°–175°E, 30°–40°N) and subtropical region(130°–175°E, 15°–30°N), respectively. Colors denote the strength of the values, which are determined by
the maximum and minimum values in each variable13

Climatic Change(2021) 168:29

Page 5 of 20

29

Fig. 1  a–d Horizontal distributions of the rate of modeled ­pHin situ change from 1999 to 2018 at 0, 100, 200,
and 500 m depth. The annual values in 1999 were subtracted from those in 2018. e Vertical sections of therate of modeled p­ Hin situ change along 165°E from 1999 to 2018

3 Results
3.1 Long‑term variations in ocean acidification indices obtained from JCOPE_EC
The rates of change in the ocean acidification index ­(pHin situ) obtained from JCOPE_EC
are shown in Fig. 1. The surface ­pHin situ declined by 0.002 ­yr–1 over most of the study
area, which means that ocean acidification is occurring over most of the surface Northwest Pacific (Fig. 1a). However, their rate changes are locally different. Prominent areas ofocean acidification (rate of change of p­ Hin situ = –0.002 to –0.004 ­yr–1) exist off the Honshu
and Ryukyu islands along Kuroshio (25°–38°N, 125°–145°E) and at the polar front region(38–42°N, 150–175°E).
Acidification is more significant at increasing depth and occurs further south at depths
below 100 m (Fig. 2b–d). The rate of change in p­ Hin situ is ≤  − 0.002 ­yr–1 in the surface
layer and − 0.002 to − 0.004 ­yr–1 at a depth of 100 m (Fig. 1b). ­pHin situ increases (i.e.,alkalinization) are also evident in some regions (Fig. 1b–d), although the areas showinga ­pHin situ increase at the surface are very limited in extent (Fig. 1a). Large increases in­pHin situ occur in the northern Japan Sea and subarctic region (35°–40°N, 140°–150°E;
Fig. 1b). The boundary between ocean acidification and alkalinization occurs at around
35°N, corresponding to the Kuroshio Extension area at a depth of 100 m (Fig. 1b), but
the boundary is further south at around 15°–20°N at depths of 200 and 500 m (Fig. 1c–d).
There are some prominent areas of alkalinization (140°–175°E, 38°–50°N) in the northwestern Japan Sea and northern subarctic region off the northeast coasts of Hokkaido andKuril Islands (Fig. 1c–d).
Vertical sections of ­pHin situ along 165°E (Fig. 1e) reveal that these patterns of acidification and alkalinization extend vertically, as well as horizontally, with rates of –0.008 to
0.004 ­yr–1. Ocean acidification is occurring at rate of ≤ 0.002 ­yr–1 at the surface, and strong
acidification is occurring in the subtropical region south of 25°N at depths of 200–800 m.

13

29

Page 6 of 20

Climatic Change(2021) 168:29

Fig. 2  Same as for Fig. 1 but for modeled DIC. Density contours in (e) are for 1999

Alkalinization is occurring in the area north of 25°N at depths below 100 m, especiallyin the subarctic region (Fig. 1e). Strong acidification or alkalinization is occurring below
200-m depth. Areas of enhanced acidification and alkalinization are found vertically in theregions from 15°–25°N to 40°–45°N, respectively (Fig. 1e). Although not shown here, theother ocean acidification indices ­(pH25 and Ωarg) also exhibit similar patterns of acidification and alkalinization as the p­ Hin situ values (Fig. 1).

3.2 Comparison between the observational DIC data and model results
A primary controlling factor on ocean acidification indices is DIC (Ishizu et al. 2019).
Therefore, we also show the rate of change of the modeled DIC (Fig. 2). The surface ocean
exhibits a large-scale increase in DIC (Fig. 2a), but the DIC changes are not spatially uniform (Fig. 2a–b). The patterns are similar to those of the ocean acidification indices shownin Fig. 1, which indicates that the DIC increase and decrease reflect the pattern of ocean
acidification and/or alkalinization.
In our study region, we have obtained DIC measurements along 165°E. Although the
DIC observations exhibit seasonal, interannual, and bi-decadal variability due to temporalvariations in the collection of the DIC measurements (Cook et al. 1997; Biondi et al. 2001),
the observational DIC data (Fig. 3a–f) show increases and/or decreases that are consistent
with the simulation results, including: (1) the DIC increase to be more significant withincreasing depth in the subtropical region (15°–30°N); (2) the DIC decrease that spreads
partially and/or to all depths in the subarctic region (35°–48°N); (3) some of the areas witha DIC increase being in the subarctic region (Fig. 3a, b, d, e); (4) the DIC increase in thesubarctic region being located near the surface layer or close to the Kuroshio Extension(35°–40°N; Fig. 3a–f).

13

Climatic Change(2021) 168:29

Page 7 of 20

29

Fig. 3  Vertical sections of the rate of observed DIC change from (a–c) 2018 and (d–f) 2017 to 1992, 2000,
2006, and 2011 along 165°E. (g) Same as for Fig. 2e, but for the change in the rate of modeled DIC from1995 to 2018

The vertical distributions of the simulated DIC and observational DIC changes (Figs. 2,3a–f) show that the area of DIC decrease below the subsurface is smaller in the model thanin the observations. The vertical section of the DIC differences between 1995 and 2018
may include effects from the spin-up process (Fig. 3g). However, the agreement between
the DIC changes obtained from the model and observations is better for 1995–2018 than
1999–2018 (Figs. 2e and 3).

4 Discussion
4.1 Comparison with observational data
Midorikawa et al. (2010) found that ocean acidification is occurring in the surface watersof the subtropical region, and estimated p­ H25 changes of 0.0013 ± 0.0005 ­yr–1 for summer and 0.0018 ± 0.0002 ­yr–1 for winter along 137°E. Our model estimates that the surface­pH25 changes at a rate of 0.0019 ­yr–1 at 30°N and 137°E, which is consistent with theestimate of Midorikawa et al. (2010). In the northern part of our study area, Wakita et al.
(2013) estimated that ocean acidification was occurring at a rate of 0.0011 ± 0.0004 ­yr–1
in the surface mixed layer, and 0.0051 ± 0.0010 ­yr–1 in the 200 − 250-m depth range at the

13

29

Page 8 of 20

Climatic Change(2021) 168:29

K2 station (47°N, 160°E). However, our simulation results in the North Pacific suggest
alkalinization is occurring at depths greater than 100 m, and that even stronger alkalinization is occurring below 100 m in the northern Northwest Pacific and Japan Sea (Fig. 1).
The differences between the observations and modeling studies could be due to the timerequired to confidently detect (anthropogenic) trends (McKinley et al. 2016; Schlunegger
et al. 2019). In general, modeling studies indicate that an ~ 20-year time-series of carbonchemistry related to sea-surface properties is needed to robustly identify anthropogenic
trends (McKinley et al. 2016; Schlunegger et al. 2019).
Several research programs are measuring state variables related to ocean acidification
indices. The World Ocean Circulation Experiment (WOCE) includes cruise lines thatobtain annual observations from some of these programs (Carter et al. 2017; Ono et al.
2019). The Ocean Carbon and Acidification Data Portal is currently organizing a system
for gathering observational data and making it publicly available (https://​www.​nodc.​noaa.​
gov/​oads/​stewa​rdship/​data_​portal.​html). The additional data from these programs may provide an improved three-dimensional view of acidification and alkalinization in the global
oceans.

4.2 Air–sea ­CO2 flux
The horizontal distribution of annual air–sea C­ O2 fluxes (Fig. 4) shows that most of the
Pacific plays a role in outgassing C­ O2 to the atmosphere, except for marginal seas (i.e.,
the Japan and East China seas). The polar frontal area (40° − 45°N, 140° − 180°E) also
appears to be undergoing active uptake of C
­ O2. Several historical studies have noted thepoor agreement between measured and modeled p­ CO2 values which is calculated from
DIC and ALK (Orr and Epitalon 2015; Valsala and Maksyutov 2010; Ishizu et al. 2019,2020). In our study region, the calculated oceanic ­pCO2 value is larger than the climatological value by –70 to 70 μatm (Valsala and Maksyutov 2010). The modeled absorption
rate in the Kuroshio–Oyashio transition area therefore may be smaller than the observed
absorption rate and then the area of the DIC increase in the model tends to be smaller thanthat suggested by the observational data (Fig. 2).

Fig. 4  Annual mean air–sea C­ O2flux of the modeled outputs from1999 to 2018. Positive and negative values represent ­CO2 uptake
and release, respectively

13Climatic Change(2021) 168:29

Page 9 of 20

29

Fig. 5  Time-series of annual
air–sea ­CO2 fluxes of the model
outputs from 1999 to 2018 in thesubarctic region (blue), Kuroshio
Extension (red), and subtropicalregion (gray). The annual air–sea­CO2 flux is divided into thesubarctic region (155°–175°E,
40°–50°N), Kuroshio Extension(142°–175°E, 30°–40°N), andsubtropical region (130°–175°E,
15°–30°N). Positive and negativevalues represent ­CO2 uptake andrelease, respectively

The annual air–sea C­ O2 flux from the ocean to atmosphere (Fig. 5) gradually decreased
from 1999 to 2018 in the subarctic region, Kuroshio Extension, and subtropical region.
The statistical regression analyses of the subarctic region, Kuroshio Extension, and subtropical region yielded trends of 0.0015, 0.0015, and 0.0006 Pg C ­yr–1, respectively. The
p-values for each region were 6.9 × ­10–5, 5.5 × ­10–2, and 1.4 × ­10–1, respectively. Since all
the regressions are positive, the C­ O2 release in the subarctic region, Kuroshio Extension,and subtropical region decreases with increasing atmospheric C­ O2 concentrations. Only the
p-values for the Kuroshio Extension and subtropical region are > 0.05, indicating that thestatistical regression in the subarctic region is not reliable. Large variability in the air–sea­CO2 fluxes in the Kuroshio Extension (Fig. 5) is likely driven by strong eddy activity.

4.3 Possible mechanisms of long‑term changes in ocean acidification
and alkalinization in the Northwest Pacific represented by JCOPE_EC
We further examined the processes leading to either acidification or alkalinization between
the subtropical and subarctic regions of the Northwest Pacific (Figs. 1, 2). To do this, we
calculated the mean balances of the DIC variation terms generated by physical and biological processes, as follows (Ishizu et al. 2019):
𝜕[DIC]
𝜕[DIC]
𝜕[DIC]
𝜕[DIC]
𝜕[DIC]
𝜕[DIC]
=[
] +[
]
+[
]
+[
] +[
]
,
𝜕t
𝜕t
𝜕t
𝜕t
𝜕t
𝜕t
A
xy_difz_dif
Bio
air−sea(2)where the subscripts A, xy_dif, z_dif, Bio, and air–sea represent the time derivatives of DIC
induced by advection, horizontal diffusion (i.e., horizontal mixing), vertical mixing, biological processes, and air–sea ­CO2 exchange, respectively (positive air–sea values indicate
a transfer of ­CO2 from air to sea). We refer to these as the “DIC variation terms”, and the
DIC variations induced by air–sea ­CO2 exchange were only included for the surface level.
The total DIC variation term on the left-hand side of Eq. (2) is the sum of all terms on theright-hand side of the equation. The averaged terms for the period from 1999–2018 were(left-hand side of Eq. (2)) isanalyzed. The temporally averaged total variation term 𝜕[DIC]𝜕tthus not the same as the change rate term (Eq. (1)), which was calculated from the difference between the annual mean for 1999 and 2018.
The horizontal distributions of the DIC variation terms in Eq. (2) are shown in Figs. 6,
7, 8, while Fig. 9 displays the depth-dependent variation in the DIC variation terms

13

29

Page 10 of 20

Climatic Change(2021) 168:29

Fig. 6  Horizontal distributions of annual modeled DIC variations generated by (a, g) advection, (b, h) horizontal diffusion, (c, i) vertical mixing, (e, j) biological processes, and (e) air–sea C­ O2 exchange, and (f, k)
total modeled DIC temporal variations at (a–f) 0 m and (g–k) 50 m depths. Positive and negative values
represent increases and decreases that each term produces on the modeled DIC variations, respectively. The
modeled DIC variation due to air–sea ­CO2 exchange was only considered in surface waters

along the 165°E line at 20°, 30°, 40°, and 50°N. The balances at the surface (Fig. 6a–f)
indicate that air–sea exchange is generally balanced by vertical mixing (Figs. 6c, e and9a–d), and that biological processes drive the DIC increase in the Kuroshio Extension
and subtropical region, as well as the DIC decrease in the subarctic region north of 50°N
(Figs. 6d and 9a–d). This DIC increase in the Kuroshio Extension and subtropical region
can be attributed to remineralization, which occurs due to the spreading of detritus originating from zooplankton (not shown). The DIC decrease in the subarctic region is due to13

Climatic Change(2021) 168:29

Page 11 of 20

29

Fig. 7  Same as for Fig. 6, but the horizontal distributions of annual modeled DIC variations generated by(a, e) advection, (b, f) horizontal diffusion, (c, g) vertical mixing, (e, h) biological processes, and (d, i) totalmodeled DIC temporal variations at (a–d) 100 m and (e–i) 200 m depths

photosynthesis. The balance of the annual DIC variations terms along the 165°E line in
the subtropical region and surface region of the Kuroshio Extension (Fig. 9a, b, e) suggests that the DIC supply, which is primarily controlled by vertical mixing, surpasses the
DIC release by air–sea exchange. As the depth increases (Figs. 6g–k, 7, 8, and 9e–s), the
DIC temporal variations are primarily driven by vertical mixing, biological processes,
and/or advection. Vertical mixing transports DIC from deeper to subsurface waters(above 50 m) in the subtropical region and Kuroshio Extension (Figs. 6c, j and 9e, f). At
100-m depth, the positive DIC values from the subarctic to subtropical regions (Fig. 7d)
indicate production of DIC by remineralization. Such latitudinal differences between thesubtropical and subarctic regions due to biological processes (Fig. 7d, i) can be linked tolatitudinal differences in the depth of the chlorophyll maximum in the Northwest Pacific(Sauzede et al. 2015; Ishizu et al. 2020). Below 200 m depth (Figs. 7f–j and 8), biological processes become less important in driving DIC variations. Instead, advection

13

29

Page 12 of 20

Climatic Change(2021) 168:29

Fig. 8  Same as for Fig. 6, but the horizontal distributions of annual modeled DIC variations generated by(a, e) advection, (b, f) horizontal diffusion, (c, g) vertical mixing, (e, h) biological processes, and (d, i) totalmodeled DIC temporal variations at (a–d) 500 m and (e–i) 1000 m depths

is largely responsible for DIC variations via vertical mixing (Figs. 7f, h, j and 9m–s).
Horizontal mixing does not appear to have an important role in driving lateral or depthdependent DIC variations across the study area, especially in the surface and subsurface(Figs. 7, 8).
DIC observations (Fig. 3) suggest that ocean acidification in the subtropical region
at the surface, as well as in deeper layers, has progressed more rapidly than in thesubarctic region. Results from the simulation suggest that vertical mixing balances theemission of DIC at the surface, despite ­CO2 release to the atmosphere. However, in thesubtropical region, the positive DIC values due to remineralization also occur at thesurface (Fig. 9a–b), because the photosynthesis zone occurs to depths of 150–200 m(Ishizu et al. 2019). Below 100-m depth, the model results suggest that the supply from
DIC transport due to vertical mixing and/or advection in the subarctic region below

13

Climatic Change(2021) 168:29

Page 13 of 20

29

Fig. 9  Histograms of annual modeled DIC variations at 165°E. Partitioning of DIC variations into contributions from advection, horizontal diffusion, vertical mixing, biological processes, and air–sea ­CO2 exchange,
and total modeled DIC temporal variations, at (a–d) 0 m, (e–h) 50 m, (i–l) 100 m, (m–p) 200 m, (q–t)
500 m, and (u–s) 1000 m depths at the (a, e, i, m, q, u) 20°E, (b, f, j, n, r, v) 30°E, (c, g, k, o, s, w) 40°E,
and (d, h, l, t, s) 50°E. The annual mean values were calculated within five grids (22.0–45.5 km) of the target region. The colored bars have the same colors as the processes in shown in Figs. 6, 7 and 8

100 m is less than that of the upper layers (Figs. 6, 7, 8 and 9k–l, n–p, r–t, v–s). This
process is responsible for the DIC decrease in the subarctic region (Figs. 2b–d and3), and likely leads to the alkalinization below the subsurface layers (Fig. 1b–e). DIC

13

29

Page 14 of 20

Climatic Change(2021) 168:29

Fig. 10  Schematic illustrationof the formation of the main
water masses and near-surface
currents in the Northwest Pacific.
Red, green, and blue ellipses
denote the formation regions of
the STMW, S-CMW, D-CMW,
OSMW, and DWS, respectively.
Colored arrows show the current
paths of the Kuroshio, OSMW,
and DSW. Colored dashedarrows show the spreading of
the water masses after subduction. STMW = SubtropicalMode Water; S-CMW = Shallower Central Mode Water;
D-CMW = Deeper Central Mode
Water; OSMW = Ohotsk Sea
Mode Water; DSW = Dense Shelf
Water; NPIW = North PacificIntermediate Water

Fig. 11  A schematic of theacidification and alkalinization
mechanisms across different latitudes in the Northwest
Pacific, which were obtained
from the simulation results. Red
and blue colored areas represent
acidification and alkalinization,
respectively. The green coloredarea shows the phytoplankton
zone, where modeled DIC isconsumed in the surface and subsurface layers. The blue arrows
indicate the movement of DIC
and carbon input due to air–sea­CO2 exchange

is transported to the Kuroshio Extension and subtropical region below 200-m depths
and is involved in the formation and transport of mode waters in the Northwest Pacific(Figs. 1 and 10), such as the Central Mode Water (S-CMW, 26.0–26.4 σθ; D-CMW,
26.4–26.7 σθ; Suga et al. 1997), Subtropical Mode Water (STMW, 25.2–25.8 σθ; Masuzawa 1969), and NPIW (26.7–27.0 σθ; Sverdrup et al. 1942; Talley 1993; Yasuda 1997,
2004; Fujii et al. 2013), which occur in the 200–600-m depth range (Figs. 2c–e, 3).
The density contours (Fig. 2e) imply that the DIC in the subarctic region is transported
to the subtropical region by the water masses of mode waters. Figure 11 illustrates themechanisms of acidification and alkalinization in a vertical section in the Northwest
Pacific.

13

Climatic Change(2021) 168:29

Page 15 of 20

29

5 Conclusions
A biogeochemical model with an operational ocean product was used to examine longterm three-dimensional variations in ocean acidification indices and state variables
related to ocean acidification from 1993 to 2018. The model results reveal an increasein DIC and ocean acidification over most of the surface layer. The Northwest Pacificbecame less of a C­ O2 source to the atmosphere from 1999 to 2018. Our results reveal
differences between the subtropical and subarctic regions, particularly at increasing
depths. Ocean acidification is occurring in the subtropical region, while alkalinization
is occurring in the subarctic region. The model indicates that ocean acidification represented by our model is strongly controlled by the physical transport of DIC by ocean
circulation and biological activity. Historical DIC observations in the Northwest Pacificshow similar patterns across the Kuroshio Extension region. At depths below 100 m, we
identified ocean acidification in the subtropical region and alkalinization in the subarctic region, where advection, horizontal transport, and vertical diffusion are important
processes. Biological processes also modulate the physical processes above depths of
150–200 m in the subtropical region.
Although our model experiment does not distinguish between DIC changes due toocean circulation and those due to increased sequestration of atmospheric ­CO2 as a consequence of increasing atmospheric ­CO2, our model does allow us to examine the mechanisms responsible for the stronger ocean acidification in the subtropical gyre south of theKuroshio Extension relative to the subarctic region. C­ O2 may be sequestered by the Kuroshio Extension in winter by subduction processes (not shown), and the injected ­CO2 in thesubsurface and deeper layers may gradually reach the surface, with possibly more rapidtransport to the surface in the subtropical region. Our simulation results show a sustained
amplification of DIC and acidification in the subtropical region and along the coastal area
south of 35°N in the Northwest Pacific, and it is uncertain if this trend will persist into the
future. The possible sustained amplification of spatially heterogeneous acidification rates
has important consequences for marine life and mitigation strategies.

Appendix Spin‑up processes
The averaged DIC time-series in the subarctic region, Kuroshio Extension, and subtropical region from 1993 to 2008 were calculated in the surface (0–100 m), subsurface (100–500 m), middle (500–1000 m), and deeper (1000–2000 m) layers. The DIC
time-series in the surface layer (Fig. 12a) varies more than those in the deeper layers (Fig. 12b–d). In the subsurface layer (Fig. 12b), DIC values decrease in the subarctic region, but increase in the Kuroshio Extension and subtropical region. In the layers
deeper than 500 m (Fig. 12c–d), a widespread decrease in DIC is observed. The DIC
changes gradually with time and its fluctuations become smaller, except in surface waters(Fig. 12a). This DIC time-series from 1993 to 2018 was used to define the spin-up process
for the model output prior to 1999.
Figures 13 and 14

13

29

Page 16 of 20

Fig. 12  Time-series of annualmodeled DIC variations from1993 to 2018 in each layer inthe subtropical region (black),
Kuroshio Extension (green), andsubtropical region (red). (a–d)
The time-series in the surface(0–100 m depth), subsurface(100–500 m depth), middle(500–1000 m depth), and deeper
layers (1000–2000 m depth),respectively

13

Climatic Change(2021) 168:29

Climatic Change(2021) 168:29

Page 17 of 20

29

Fig. 13  Annual mean values and linear regression trends in (a–c) modeled DIC and (d–f) modeled ­pHin situ(d–f) in the subtropical region, Kuroshio Extension, and subarctic region, respectively (Table 1). Red,
green, blue, purple, and black colors are data for 0, 100, 200, 500, and 1000 m, respectively. Regression
coefficients and p-values are given in parentheses on the right side of the figure

13

29

Page 18 of 20

Climatic Change(2021) 168:29

Fig. 14  Same as for Appendix Fig. 13, but for modeled ­pH25 and Ωarg
Acknowledgements We acknowledge support from the Sasakawa Peace Foundation of the Ocean Policy
Research Institute (OPRI-SPF). This study is part of the Japan Coastal Ocean Predictability Experiment(JCOPE) supported by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as longas you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article
are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to thematerial. If material is not included in the article’s Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directlyfrom the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

References
Bates NR (2007) Interannual variability of the ocean C­ O2 sink in the subtropical gyre of the North Atlantic
Ocean over the last 2 decades. J Geophys Res 112:C09013. https://​doi.​org/​10.​1029/​2006J​C0037​59
Bates NR, Astor YM, Church MJ, Currie K, Dore JE, Gonzalez-Davila M, Lorenzoni L, Muller-Karger F,
Olafsson J, Santana-Casiano JM (2014) A time-series view of changing surface ocean chemistry due
to ocean uptake of anthoropogenic ­CO2 and ocean acidification. Oceanography 27(1):126–141. https://​doi.​org/​10.​5670/​ocean​og.​2014.​16.​201413

Climatic Change(2021) 168:29

Page 19 of 20

29

Biondi F, Gershunov A, Cayan RD (2001) North Pacific decadal climate variability since 1661. J Clim
14:5–10Carter BR, Feely AR, Mecking S, Cross JN, Macdonald MA, Siedlecki AS, Talley DL, Sabine LC, Millero
JF, Swift HJ, Dickson GA, Rodgers BK (2017) Two decades of Pacific anthropogenic carbon storageand ocean acidification along global ocean ship-based hydrographic investigations program sections
P16 and P02. Glob Biogeochem Cycles 31:306–327. https://​doi.​org/​10.​1002/​2016G​B0054​85
Chang YL, Miyazawa Y, Guo X (2015) Effects of the STCC eddies on the Kuroshio based on the 20-year
JCOPE2 reanalysis results. Prog Oceanogr 135:64–76. https://​doi.​org/​10.​1016/j.​pocean.​2015.​04.​006
Chen CTA, Lui H-K, Hsieh C-H, Yanagi T, Kosugi N, Ishii M, Gong G-C (2017) Deep oceans may acidify faster than anticipated due to global warming. Nat Commun:890–895.https://​doi.​org/​10.​1038/​
s41558-​017-​0003-y
Cook ER, Meko MD, Stockton WC (1997) A new assessment of possible solar and lunar forcing of bidecadal drought rhythm in the western United States. J Clim 10:1343–1356
Dore JE, Lukus R, Sadler DW, Church MJ, Karl DM (2009) Physical and biogeochemical modulation ofocean acidification in the central North Pacific. Proc Natl Acad Sci 106:12235–12240
Fujii Y, Nakano T, Usui N, Matsumoto S, Tsujino H, Kamachi M (2013) Pathways of the North PacificIntermediate Water identified through the tangent linear and adjoint models of an ocean general circulation model. J Geophys Res Oceans 118:2035–2051. https://​doi.​org/​10.​1002/​jgrc.​20094
Gonzalez-Davila M, Santana-Casiano JM, Gonzalez-Davila EF (2007) Interannual variability of the upper
ocean carbon cycle in the northeast Atlantic Ocean. Geophys Res Lett 34:L07608. https://​doi.​org/​10.​1029/​2006G​L0281​45
Intergovernmental Panel on Climate Change (IPCC) (2013) Climate Change 2013: The Physical Science
Basis Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panelon Climate Change. In: Stocker TF, Qin D, G-K Plattner, Tignor MMB, Allen SK, Boschung J, Nauels A, Zia Y, Bex V, Midgley PM (eds) Cambridge University Press, Cambridge and New York, pp
1−1535Ishii M, Rodgers BK, Inoue YH, Toyama K, Sasano D, Kosugi N, Enyo K, Nakano T, Iudicone D, Blanke
B, Aumount O, Feely AR (2020) Ocean acidification from below in the tropical Pacific. Biogeochem
Cycles 34:e2019GB006368. https://​doi.​org/​10.​1029/​2019G​B0063​68
Ishizu M, Miyazawa Y, Tsunoda T, Guo X (2019) Development of a biogeochemical and carbon model
related to ocean acidification indices with an operational ocean model product in the North Western
Pacific. Sustainability 11:2677. https://​doi.​org/​10.​3390/​su110​92677
Ishizu M, Miyazawa Y, Guo X, Tsunoda T (2020) Seasonal variability in the inorganic ocean carbon cyclein the Northwest Pacific evaluated using a biogeochemical and carbon model coupled with an operational ocean model. Clim Chang:1–26. https://​doi.​org/​10.​1007/​s10584-​020-​02779-2
Masuzawa J (1969) Subtropical Mode Water. Deep-Sea Res 16:463–472
McKinley A, Galen Pilcher J. D, Fay RA, Lindsay K, Long CM, Lovenduski SN (2016) Timescales for
detection of trends in the ocean carbon sink. Nature 530:469–472. https://​doi.​org/​10.​1038/​natur​e16958
Mellor GL (2001) An equation of state for numerical models of oceans and estuaries. J Atmos Ocean Technol 8:609–611
Midorikawa T, Ishii M, Saito S, Sasano D, Kosugi N, Motoi T, Kamiya H, Nakadate A, Nemoto K, Inoue H(2010) Decreasing pH trend estimated from 25-yr time series of carbonate parameters in the western
North Pacific. Tellus 62B:649–659
Miyama T, Minobe S, Goto H (2021) Marine heatwave of sea surface temperature of the Oyashio region in
the summer in 2010–2016. Front Mar Sci 7:576240. https://​doi.​org/​10.​3389/​fmars.​2020.​576240
Miyazawa Y, Zhang R, Guo X, Tamura H, Ambe D, Lee JS, Yoshinari H, Setou T (2009) Water mass variability in the western north Pacific detected in a 15-year eddy resolving ocean reanalysis. J Oceanogr
65:737–756
Miyazawa Y, Varlamov MS, Miyama T, Guo X, Hihara T, Kiyomatsu K, Kachi M, Kurihara Y, Murakami
H (2017) Assimilation of high-resolution sea surface temperature data into an operational nowcast/
forecast system around Japan using a multi-scale three-dimensional variational scheme. Ocean Dyn
67:713–728
Miyazawa Y, Kuwano-Yoshida A, Doi T, Nishikawa H, Narazaki T, Fukuoka T, Sato K (2019) Temperature
profiling measurements by sea turtles improve ocean state estimation in the Kuroshio-Oyashio Confluence region. Ocean Dyn 69:267–282
Olafsson J, Olafsdottir SR, Benoit-Cattin A, Danielsen M, Arnarson TS, Takahashi T (2009) Rate of Iceland
Sea acidification from time series measurements. Biogeosciences 6(2):661–2, 668. https://​doi.​org/​10.​
5194/​bg-6-​2661-​2009

13

29

Page 20 of 20

Climatic Change(2021) 168:29

Ono H, Kosugi N, Toyama K, Tsujino H, Kojima A, Enyo K, IidaNakano YT, Ishii M (2019) Accelerationof ocean acidification in the western North Pacific. Geophys Res Lett 46:13161–13169. https://​doi.​org/​
10.​1029/​2019G​L0851​21
Orr CJ, Epitalon MJ (2015) Improved routines to model the ocean carbonate system: mocsy 2.0. Geosci
Model Dev 8:485–499
Orr JC, Najjar R, Sabine CL, Joos F (1999) Abiodic-HOWTO; Internal OCMIP Report; LSCE/CEA
Scalay: Gif-sur-Yvette, France; 25p
Saeyanto E, Guo X, Ono J, Miyazawa Y (2014) Interannual variations of Kuroshio transport in the East
China Sea and its relation to the Pacific decadal oscillation and mesoscale eddies. J Geophys Res
Oceans 119:3595–3616. https://​doi.​org/​10.​1002/​2013J​C0095​29
Sauzede R, Lavigne H, Claustre H, Schmechtig C, Ortenzio DF, Guinet C, Pesant S (2015) Vertical
distribution of chlorophyll a concentration and phytoplankton community composition from in situ
flurorescence profiles: a first database for the global ocean. Earth Syst Sci Data 2015(7):261–273.
https://​doi.​org/​10.​5194/​essd-7-​261-​2015
Schlunegger S, Rodgers BK, Sarmiento LJ, Frolicher LT, Dunne PJ, Ishii M, Slater R (2019) Emergenceof anthropogenic signals in the ocean carbon cycle. Nat Clim Chang 9:719–725. https://​doi.​org/​10.​1038/​s41558-​019-​0553-2
Suga T, Takai Y, Hanawa K (1997) The subtropical mode water circulation in the North Pacific. J Phys
Oceanogr 25:958–970
Sverdrup H, Johnson WM, Flemming HR (1942) The Oceans. Their Physics. Chemistry and Biology.
Prentice-Hall, 1087pp
Takahashi T, Sutherland CS, Sweeney C, Poisson A, Metzl N, Tilbrook B, Bates N, Wanninkhof R, Feely
AR, Sabin C, Olafsson J, Nojiri Y (2002) Global sea-air ­CO2 flux based on climatological surface
ocean ­pCO2, and seasonal biological and temperature effects. Deep-Sea Res II 49(49):1601–1622
Takahashi T, Sutherland CS, Wanninkhof R, Sweeney C, Feely AR, Chipman WD, Hales B, Friedrich
G, Chaves F, Sbaine C, Watson A, Bakker CED, Schuster U, Metzl N, Yoshikawa-Inoue H, Ishii
M, Midorikawa T, Nojiri Y, Kortzinger A, Steinhoff T, Hoppema M, Olafsson J, Arnarson ST, Tilbrook B, Johannessen T, Olsen A, Bellerby R, Wong CS, Delille B, Bates NR, de Baar JWH (2009)
Climatological mean and decadal change in surface ocean ­pCO2 and net sea-air ­CO2 flux over the
global oceans. Deep-Sea Res II 56(56):554–577
Talley LD (1993) Distribution and formation of North Pacific Intermediate Water. J Phys Oceanogr
23:517–537. https://​doi.​org/​10.​1175/​1520-​0485
Toyama K, Rodger BK, Blanke B, Iudicone D, Ishii M, Aumont O, Sarmiento LJ (2017) Large Reemergence of anthropogenic carbon into the ocean’s surface mixed layer sustained by the ocean’s overturning circulation. J Clim 30:8615–8631. https://​doi.​org/​10.​1175/​JCLI-D-​16-​0725.1
Valsala V, Maksyutov S (2010) Simulation and assimilation of global ocean ­pCO2 and air-sea ­CO2 fluxes
using ship observation of surface ocean ­pCO2 in a simplified biogeochemical offline model. Tellus
Ser B Chem Phys Meteorol 62(5):821–840. https://​doi.​org/​10.​1111/j.​1600-​0889.​2010.​00495.x
Wakita M, Watanabe S, Honda M, Nagano A, Kimoto K, Matsumoto K, Kitamura M, Sasaki K,
Kawakami H, Fujiki T, Sasaoka K, Nakano Y, Murata A (2013) Ocean acidification from 1997 to2011 in the subarctic western North Pacific Ocean. Biogeosciences 10:7817–7827
Yasuda I (1997) The origin of the North Pacific Intermediate Water. J Geophys Res 102(C1):893–909
Yasuda I (2004) North Pacific Intermediate Water: progress in SAGE (SubArctic Gyre Experiment) andrelated projects. J Oceanogr 60:385–395. https://​doi.​org/​10.​1023/B:​JOCE.​00000​38344.​25081.​42
Yasunaka S, Nojiri Y, Nakaoka S, Ono T, Mukai H, Usui N (2013) Monthly maps of sea surface dissolved inorganic carbon in the North Pacific: basin-wide distribution and seasonal variation. J Geophys Res Oceans 118:3843–3850. https://​doi.​org/​10.​1002/​jgrc.​20279
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps andinstitutional affiliations.

13