Long-term variations in ocean acidification indices in the Northwest Pacific from 1993 to 2018

Long-term variations in ocean acidification indices in the Northwest Pacific were examined using 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 and Kuroshio Extension regions and relative alkalinization (i.e., reduced acidification) in the subarctic 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 biological processes 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 of ocean acidification, along with the strength of the ocean carbon sink, which are key aspects of global climate change.


Introduction
The atmospheric partial pressure of CO 2 (pCO 2 , 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 of human activities such as fossil fuel burning, deforestation, and cement production IPCC 2013). The contemporary global ocean now acts as a net sink of anthropogenic CO 2 , but this net flux reflects an imbalance between regions of CO 2 outgassing and uptake. In general, the subtropics act as a source for CO 2 and the subarctic is a sink ), but there is strong regional and temporal heterogeneity in the air-sea CO 2 flux, because the sea surface CO 2 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 CO 2 concentration can become higher than that of the atmospheric CO 2 concentration, and vice versa depending on the regional internal oceanic processes. Interannual and decadal variations in the climate system can also alter the patterns of mean-state CO 2 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 CO 2 flux.
The uptake of anthropogenic CO 2 by seawater changes the natural chemical equilibrium of the seawater CO 2 -carbonate system, resulting in ocean acidification (Bates et al. 2014). In the Northwest Pacific, Midorikawa et al. (2010) estimated surface pH 25 trends from a 25-year time-series along the 137°E line south of 35°N using observational data from the Japan Meteorological Agency (JMA). The results revealed that pH 25 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) and found that pH in 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 pH in 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 of 26.9σ θ kg m -3 , and is greater than any previously reported pH in situ decrease in the surface layer in the open North Pacific.
Monitoring for biogeochemical variables including ocean acidification is important for 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 pH in situ measurements (pH 25 ) 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 (Miyazawa et al. 2017 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(Miyazawa et al. , 2017Saeyanto et al. 2014;Chang et al. 2015;Miyama et al. 2021). Using the JCOPE reanalysis model products as the physical background, Ishizu et al. (2019Ishizu et al. ( , 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 showed that the model reliably reproduces the distributions of the biogeochemical state variables (DIN, DIC, and ALK) on a seasonal scale (Ishizu et al. 2019(Ishizu et al. , 2020. In this study, we examined the recent trends of long-term (i.e., multi-decadal) variations in the inorganic marine carbon cycle, using JCOPE_EC. In Sect. 2, details of the model and 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 CO 2 fluxes and potential mechanisms of the multidecadal variations in ocean acidification are discussed. Finally, the conclusions are presented in Sect. 5.

Model and data
JCOPE_EC (Ishizu et al. 2019(Ishizu et al. , 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. 2017Miyazawa et al. , 2019, which is based on the Princeton Ocean Model with 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 (Orr et al. 1999). The JCOPE_EC model calculates pCO 2 based on temperature, salinity, DIC, and ALK (Ishizu et al. 2019(Ishizu et al. , 2020. The atmospheric CO 2 values are obtained from an empirical formula based on observational data (Ishizu et al. 2019). The reproducibility of the 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 initialized using the annual climatology at the beginning of the model run (Ishizu et al. 2019(Ishizu et al. , 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(Ishizu et al. , 2020. JCOPE_EC is described in more detail in Ishizu et al. (2019Ishizu et al. ( , 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 that changes in pH in situ (pH 25 and Ω arg ) are primarily linked to changes in DIC (not shown). To validate the reproducibility of the long-term variations in ocean acidification from the model, we compared observations of DIC from 1999 to 2018 with the model outputs along the same section (line 165°E) in Sect. 3.2. There are no available observations of ALK that encompass at least a 10-year interval within the target period from 1999 to 2018.
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 of change in state variables C (x, y, z, and t) from 1999 to 2018 in Sect. 3.1: where < > denotes the annual mean of the state variable for the specified year. We also calculated linear trends for DIC, pH in situ , pH 25 , 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 the calculation method of Eq. (1).

Long-term variations in ocean acidification indices obtained from JCOPE_EC
The rates of change in the ocean acidification index (pH in situ ) obtained from JCOPE_EC are shown in Fig. 1. The surface pH in 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 of ocean acidification (rate of change of pH in 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 pH in 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). pH in situ increases (i.e., alkalinization) are also evident in some regions ( Fig. 1b-d), although the areas showing a pH in situ increase at the surface are very limited in extent (Fig. 1a). Large increases in pH in 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 and Kuril Islands (Fig. 1c-d).
Vertical sections of pH in 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. Alkalinization is occurring in the area north of 25°N at depths below 100 m, especially in 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 the regions from 15°-25°N to 40°-45°N, respectively (Fig. 1e). Although not shown here, the other ocean acidification indices (pH 25 and Ω arg ) also exhibit similar patterns of acidification and alkalinization as the pH in situ values (Fig. 1).

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 shown in 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 temporal variations 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 with increasing 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 with a DIC increase being in the subarctic region (Fig. 3a, b, d, e); (4) the DIC increase in the subarctic region being located near the surface layer or close to the Kuroshio Extension (35°-40°N; Fig. 3a-f). 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 than in 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). Midorikawa et al. (2010) found that ocean acidification is occurring in the surface waters of the subtropical region, and estimated pH 25 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 pH 25 changes at a rate of 0.0019 yr -1 at 30°N and 137°E, which is consistent with the estimate of Midorikawa et al. (2010). In the northern part of our study area, Wakita et al.

Comparison with observational data
(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  (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 time required 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 that obtain 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.

Air-sea CO 2 flux
The horizontal distribution of annual air-sea CO 2 fluxes (Fig. 4) shows that most of the Pacific plays a role in outgassing CO 2 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 CO 2 . Several historical studies have noted the poor agreement between measured and modeled pCO 2 values which is calculated from DIC and ALK (Orr and Epitalon 2015;Valsala and Maksyutov 2010;Ishizu et al. 2019Ishizu et al. , 2020. In our study region, the calculated oceanic pCO 2 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 than that suggested by the observational data (Fig. 2). The annual air-sea CO 2 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 CO 2 release in the subarctic region, Kuroshio Extension, and subtropical region decreases with increasing atmospheric CO 2 concentrations. Only the p-values for the Kuroshio Extension and subtropical region are > 0.05, indicating that the statistical regression in the subarctic region is not reliable. Large variability in the air-sea CO 2 fluxes in the Kuroshio Extension (Fig. 5) is likely driven by strong eddy activity.

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): 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 CO 2 exchange, respectively (positive air-sea values indicate a transfer of CO 2 from air to sea). We refer to these as the "DIC variation terms", and the DIC variations induced by air-sea CO 2 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 the right-hand side of the equation. The averaged terms for the period from 1999-2018 were analyzed. The temporally averaged total variation term [DIC] t (left-hand side of Eq. (2)) is thus 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 (2) Fig. 5 Time-series of annual air-sea CO 2 fluxes of the model outputs from 1999 to 2018 in the subarctic region (blue), Kuroshio Extension (red), and subtropical region (gray). The annual air-sea CO 2 flux is divided into 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). Positive and negative values represent CO 2 uptake and release, respectively 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 and 9a-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 to 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 CO 2 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 CO 2 exchange was only considered in surface waters 1 3 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 the subtropical and subarctic regions due to biological processes (Fig. 7d, i) can be linked to latitudinal 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 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 the subarctic region. Results from the simulation suggest that vertical mixing balances the emission of DIC at the surface, despite CO 2 release to the atmosphere. However, in the subtropical region, the positive DIC values due to remineralization also occur at the surface ( 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 100 m is less than that of the upper layers (Figs. 6,7,. This process is responsible for the DIC decrease in the subarctic region , and likely leads to the alkalinization below the subsurface layers (Fig. 1b-e). DIC  (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) 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 1997Yasuda , 2004Fujii 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 the mechanisms of acidification and alkalinization in a vertical section in the Northwest Pacific.  Fig. 11 A schematic of the acidification 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 colored area shows the phytoplankton zone, where modeled DIC is consumed in the surface and subsurface layers. The blue arrows indicate the movement of DIC and carbon input due to air-sea CO 2 exchange

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 increase in DIC and ocean acidification over most of the surface layer. The Northwest Pacific became less of a CO 2 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 Pacific show 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 to ocean circulation and those due to increased sequestration of atmospheric CO 2 as a consequence of increasing atmospheric CO 2 , our model does allow us to examine the mechanisms responsible for the stronger ocean acidification in the subtropical gyre south of the Kuroshio Extension relative to the subarctic region. CO 2 may be sequestered by the Kuroshio Extension in winter by subduction processes (not shown), and the injected CO 2 in the subsurface and deeper layers may gradually reach the surface, with possibly more rapid transport 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  (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