Seasonal variability in the inorganic ocean carbon cycle in the Northwest Pacific evaluated using a biogeochemical and carbon model coupled with an operational ocean model

Here, we investigate the seasonal variability in the dissolved inorganic carbon (DIC) cycle in the Northwest Pacific using a high-resolution biogeochemical and carbon model coupled with an operational ocean model. Results show that the contribution to DIC from air–sea CO2 exchange is generally offset by vertical mixing at the surface at all latitudes, with some seasonal variation. Biological processes in subarctic regions are evident at the surface, whereas in the subtropical region these processes take place within the euphotic layer and then DIC consumption deepens southward with latitude. Such latitudinal differences in biological processes lead to marked horizontal and vertical contrasts in the distribution of DIC, with modulation by horizontal and vertical advection–diffusion processes.


Introduction
The atmospheric partial pressure of CO 2 (pCO 2 ) has been increasing at a rate of~1.8 ppm by volume (ppmv) per year in recent decades as a result of human activities such as fossil-fuel burning, deforestation, and cement production (Takahashi et al. 2009;IPCC 2013). In the preindustrial era, the ocean was generally a net source of CO 2 emissions to the atmosphere because of the mineralization of land-derived organic matter in addition to that produced by in situ production, and CaCO 3 precipitation (Mackenzie et al. 2004). Rising atmospheric CO 2 concentrations caused by fossil-fuel combustion and land-use changes (Mackenzie et al. 2004;Bauer et al. 2013;IPCC 2013;IGBP, IOC, SCOR 2013) reversed the direction of the air-sea CO 2 flux, leading the global ocean to become a net sink of anthropogenic CO 2 . The present thickness of the upper thermocline, where large amounts of anthropogenic CO 2 emissions are stored, is estimated to be of the order of a few hundred meters (Mackenzie et al. 2004). The oceanic coastal zone changed from being a source to a sink during the industrial era (Mackenzie et al. 2004;Bauer et al. 2013). Several estimates of CO 2 sinks and sources in ocean provinces (Cai et al. 2006) and/or spatially explicit typology (Laurelle et al. 2010) showed that marginal seas in the tropics are sources of CO 2 , whereas those in temperate regions and at high latitudes act as sinks (Cai et al. 2006;Laurelle et al. 2010).
Data-based estimates of variability and trends in ocean CO 2 uptake are limited by the short record of observations. Although high-quality measurements of CO 2 in surface waters and air commenced in the early 1960s, the amount of available information is still limited (Wanninkhof et al. 2013). The principal observational approaches for estimating sea-air fluxes of CO 2 are to measure ΔpCO 2 from ships (Takahashi et al. 2009;Nakaoka et al. 2013) and moorings (Sutton et al. 2017), and apply a parameterization using a function of wind speed (Wanninkhof et al. 2013). Other approaches rely on simulations made by ocean biochemistry general circulation models (OBGCMs) with parameterization of biogeochemical processes and total dissolved inorganic carbon (DIC) measurements in the ocean interior, and/or atmospheric data. However, gaps remain in the understanding of ocean CO 2 uptake, especially the spatiotemporal variability of the seasonal inorganic/organic carbon cycle, because CO 2 concentrations and other related oceanic variables are difficult to observe simultaneously, frequently, and widely. The seasonal variability in pCO 2 shows differences at a local scale (Takahashi et al. 2009;Sutton et al. 2017). Model estimates of temporal trend detection (Keller et al. 2014;Lovenduski et al. 2015) show the influence of both decadal/interannual and seasonal variabilities and suggest that the time of emergence of a trend signal is basically around 10 years in the surface but the tropical area needs more time.
The detailed processes that generate variation in the DIC of the ocean interior are still uncertain. Several studies have proposed possible mechanisms for the oceanic annual carbon cycle (Palmer and Totterdell 2001;Takahashi et al. 2002;Xiu and Chai 2013). Palmer and Totterdell (2001) discussed physical and biological mechanisms that contribute to the global annual mean carbon cycle using an ecosystem model without considering the contribution of air-sea CO 2 exchange. They reported that the effects of vertical mixing were largely offset by biological processes in the latitudinal range of 25−60°N over ocean surfaces and that the effects of advection were mostly offset by biological processes at latitudes of < 20°N. Takahashi et al. (2009) focused on the relative importance of temperature and biological effects to the global seasonal cycle of air-sea CO 2 exchange by evaluating monthly climatological maps of air-sea CO 2 flux and pCO 2 . Xiu and Chai (2013) investigated the seasonal and decadal variability of the upper-ocean carbon cycle in the North Pacific using a physical-biogeochemical model. Their results showed that the seasonal variability in pCO 2 and CO 2 flux in the North Pacific followed the change in sea surface temperature closely, with high and low values in summer and winter, respectively, and that surface pCO 2 variations at the modeled sites correspond to well-known observational monitoring points controlled primarily by anthropogenic CO 2 and modulated by decadal variations.
The Japan Coastal Ocean Predictability Experiment (JCOPE; http://www.jamstec.go. jp/jcope/) is an operational eddy-resolving physical ocean model for the Northwest Pacific, the Japan Sea, the Okhotsk Sea, and the East China Sea (Miyazawa et al. 2009(Miyazawa et al. , 2014. Ishizu et al. (2019) recently developed a biogeochemical and carbon model coupled with the JCOPE (JCOPE_EC). This model generally reproduces the observed seasonal variability of chlorophyll-a (Chl-a), dissolved inorganic nitrogen (DIN), phosphorus (DIP), and inorganic carbon (DIC), and total alkalinity (ALK), but includes monthly climatological damping for DIN, DIP, DIC, and ALK (Ishizu et al. 2019). The damping forcibly constrains the calculated biological parameters (DIN, DIP, DIC, and ALK) around the monthly climatological values with a timescale of 30 days (Ishizu et al. 2019), meaning that those authors were unable to discuss the mechanism of the inorganic carbon cycle. In this study, we therefore performed a simulation using JCOPE_EC without any climatological damping and examined the physical and biological mechanisms represented by the model dynamics, focusing on the critical roles of horizontal and vertical advection-diffusion processes in generating seasonal variation in DIC.
We present the results from JCOPE_EC (without climatological damping) and discuss the mechanisms responsible for the simulated seasonal inorganic carbon cycle for the Northwest Pacific. Details of the model configuration are given in Section 2, model accuracy is described in Section 3, the processes involved in the inorganic carbon cycle in the model are discussed in Section 4, and the conclusions of the study are provided in Section 5.

Model configuration
The JCOPE_EC (Ishizu et al. 2019) 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) based on the Princeton Ocean Model with a generalized sigma coordinate (Mellor 2001). The model is a three-dimensional high-resolution regional model covering the Northwest Pacific (108-180°E, 10.5−62°N) with a horizontal resolution of 1/12°( 4.4-9.1 km) and 46 vertical active levels. The model structure is the same as that described by Ishizu et al. (2019), but our model differs in that the governing equations for DIN, DIP, DIC, and ALK (equations 4, 5, 17, and 18 of Ishizu et al. 2019) remove climatological damping.
We determined the biogeochemical model parameters using multi-optimized operations (Ishizu et al. 2019) separately for subarctic and subtropical regions (Table 1). These biogeochemical parameters are the maximum growth rate from photosynthesis (V max ), the phytoplankton mortality rate at 0°C (M P ), the phytoplankton respiration rate (R), the maximum grazing rate of zooplankton (G z ), the zooplankton mortality rate (M) z , the decomposition rate (V PN ), and the optimum light intensity for phytoplankton (I opt ). In addition to the parameters described in Ishizu et al. (2019), we introduced latitudinal changes in V max , M P , R, M z , V PN , and I opt according to the results of several sensitivity experiments as follows: where V max , M P , R, M z , V PN , and I opt change latitudinally from 0.28 to 0.97 day −1 , 0.054 to 0.155 (mmol Nm −1 ) m −3 day −1 , 0.0011 to 0.00256 day −1 , 0.044 to 0.12 (°C) −1 , 0.0954 to 0.47 day −1 , and from 20 to 100 W m −2 , respectively; V max 0 , R 0 , M P 0 , M z 0 , V PN 0 , I opt 0 , V max 1 , R 1 , M P 1 , M z 1 , V PN 1 , and I opt 1 are the tunable parameters (Table 1); Lat bnd.vmax , Lat bnd , and Lat slp are the coefficients representing the values at latitudinal boundaries and the latitudinal slopes for these parameters, respectively (Table 1 and Eqs. 1, 3, 4, 5, and 6). The first version of JCOPE_EC had a model bias resulting in a large decrease in ALK in summer due to anomalously high CaCO 3 production with a large increase in Chl-a during summer in subarctic regions (Ishizu et al. 2019). To suppress this large decrease in ALK in this model, we set the CaCO 3 to non-photosynthetic POC production ratio to a much smaller value (0.00035) in the version of the model used here (Table 1). This improvement allows the model to well represent the seasonal variability of ALK in our target region (Section 3).
The model was driven by forcing from daily oceanic (JCOPE2M) and six-hourly atmospheric (NCEP/NCAR) reanalysis data for a 1-year period (2015). The initial concentrations of phytoplankton were set to 0.1 and 0.0 mmol N m −3 for depths above and below 150 m, respectively. The initial zooplankton concentrations were set to 10% of the phytoplankton concentration. The initial detritus concentration was set to 0.0 mmol N m −3 . The variables DIN, DIP, and DIC were initialized using the climatology for January, and ALK was initialized using the annual climatology (Ishizu et al. 2019).

Model validation
To validate phytoplankton concentrations in the model, MODIS-Aqua Ocean Color Data for Chl-a from 2015 were used, as downloaded from the website of the Physical Oceanography Distributed Active Archive Center (PODAAC). To validate the model results for DIN, DIP, and DIC, we used monthly climatological DIN, DIP, and DIC data (World Ocean Atlas 2013 (WOA13); Yasunaka et al. 2013) and Japan Meteorological Agency (JMA) observational data for 2015, as in Ishizu et al. (2019). There are no applicable monthly climatological datasets for ALK in our target region (Goyet et al. 2000;Key et al. 2004;Takatani et al. 2014, Takahashi et al. 2014. We therefore used ALK observational data obtained by JMA in 2015 for comparison with model results (Section 3).

Accuracy of modeled Chl-a, DIN, DIC, and ALK
Model results presented here are slightly less accurate than those of the model with the climatological conditions described by Ishizu et al. (2019), except for ALK. However, the Growth efficiency of zooplankton 0.  (Table 2), but the reproducibility of ALK is much improved, especially in the subarctic region, the Kuroshio Extension, and the Japan Sea. The correlation between observed ALK and simulated ALK is low and negative (R = − 0.24) in the Japan Sea because the simulated ALK values there have near-uniform values with depth (not shown).

Ocean acidification indices pH in situ , pH 25 , and Ω arg
The ocean acidification indices pH in situ , pH 25 , and aragonite saturation (Ω arg ) were calculated from model results for temperature, salinity, DIC, and ALK. The pH in situ values change throughout the year (Figs. 6 and 7). The summer pH in situ values in the subarctic region (150-175°E, 50−60°N) are slightly higher (8.05-8.10) than those shown in Fig. 9    The seasonal pattern of Ω arg (Fig. 6(i-l)) is similar to that of pH 25 (Fig. 6(e-h)). The Ω arg values in the same area of pH 25 north of 35°N become 0.5-1.0 larger in summer and in autumn. The summer increase in Ω arg is less evident in the south.
Correlation coefficients for modeled pH in situ , pH 25 , and Ω arg are higher than those reported by Ishizu et al. (2019) because the accuracy of the modeled ALK concentrations is much improved compared to them.

Processes affecting the seasonal inorganic carbon cycle
To identify processes affecting the inorganic carbon cycle, we examined the physical and biological mechanisms underlying the seasonal inorganic carbon cycle in the Northwest Pacific using the model results. We separately evaluate each process included in the governing equation as follows: Model outputs of surface chlorophyll-a data in 165E Line 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 (positive values indicate a transfer of CO 2 from air to sea), respectively. Note that we refer to these as "DIC variation terms" here, and the DIC variations induced by air-sea CO 2 exchange are included only for the top (surface) level and not for the levels below it. The total DIC variation term on the left side of Eq. (7) is a summation of all the terms on the right-hand side of the equation. Air-sea CO 2 exchange shows negative values (CO 2 emission to the atmosphere) in winter north of 35°N ( Fig. 8(a)) and positive values (CO 2 absorption from the atmosphere) south of 35°N, and vice versa in summer for these regions. The subarctic region releases CO 2 to the atmosphere in winter north of 35°N ( Fig.  8(a)), and absorbs CO 2 from the atmosphere in summer north of 35°N (Fig. 8(c)). In contrast, the subtropical region south of 35°N intensely absorbs CO 2 south of the Kuroshio Extension in winter (Fig. 8(a)). A weak release of CO 2 from the ocean to the atmosphere occurs south of 40°N in summer. The areas of strongest release are located in the northeast of the Kuril Islands in the Oyashio region and in the Okhotsk Sea in autumn-winter ( Fig. 8(a, d)). (Further details on air-sea CO 2 exchange represented in the model are described in the Appendix.) Horizontal distributions and monthly mean balances of the DIC variation terms in the surface layer induced by all subprocesses (Figs. 8 and 9(a−d)) indicate that the air-sea CO 2 exchange shown in Fig. 8(a−d) is generally offset by vertical mixing. Biological processes, however, made a subordinate contribution to the total DIC variation (Figs. 8 and 9(c, d, g, h, j, k, m)). Negative biological process values indicate consumption of DIC through photosynthesis. Biological process contributions vary spatially and with depth (Figs. 8 and 9; Appendix Figs. 16-18). At the surface level, the DIC consumption induced by biological processes is high in the subarctic region around 50°N throughout the year (Figs. 8(m−p) and 9(a−d)). The peak consumption of DIC moves deeper beneath the surface southward ( Fig. 9(d, e-g, j, m); Appendix Figs. 17 and 18). The highest DIC consumption occurs at 50-100 m depth at 30-40°N ( Fig. 9(g, j, k); Appendix Figs. 16(m-p) and 17(m-p)) and at 200 m depth in the subtropical region south of 30°N ( Fig. 9(m); Appendix Fig.  18(m-p)). The contributions from biological processes at the surface and at 200 m depth are of opposite sign (Fig. 8(m-p); Appendix Fig. 18(m-p)). Instead of the contributions of the vertical mixing or biological processes below the surface, we see that advection processes are relatively contributed to the total DIC variation ( Fig. 9(ik, n, m)).
The relative contribution of biological processes compared with the other terms at 165°E is calculated as as shown in Fig. 10. The highest DIC consumption occurs above 100 m depth north of 40°N throughout the year, gradually deepening in the range 30-40°N, and spreading vertically at 100-350 m depth in the subtropical region at latitudes south of 30°N. These patterns of DIC consumption/production may be caused by the latitudinal difference of the Chl-a maximum depth (Ishizu et al. 2019;Sauzede et al. 2015). Vertical distributions of Chl-a along 165°E from Ishizu et al. (2019) and global ocean Chl-a data from Sauzede et al. (2015) indicate that the Chl-a maximum deepens southward; the Chl-a maximum is located in the surface layer in subarctic regions and at~150 m depth in subtropical regions. The highest Chl-a consumption in the subtropical region from our results (Fig. 10) occurs at greater depths (200-500 m), but the magnitude is relatively low and makes a negligible contribution to the total variation in DIC ( Fig. 9(n−p)).

Mechanisms of the seasonal inorganic carbon cycle
The relative contributions of terms in the governing equation of DIC (Figs. 8 and 9; Appendix Figs. 16-18) suggest possible mechanisms for the seasonal carbon cycle, as follows. CO 2 (DIC) is absorbed from the atmosphere to the ocean during winter south of 40°N. The corresponding absorbed volume of DIC is conveyed from the surface to the subsurface by vertical mixing (Figs. 8(i,j) and 9(a−c)). In summer and autumn, CO 2 (DIC) is released to the atmosphere by air-sea interactions, and vertical mixing generally offsets carbon emissions ig. 8 Surface horizontal distributions of monthly mean DIC variation terms generated by air-sea CO 2 exchange to the ocean (a-d), advection (e-h), vertical mixing (i-l), biological processes (m-p), and total DIC time variation terms (q-t) at the surface (0 m depth) for January, April, July, and October, respectively (Figs. 8(c,d) and 9(a−c)). In the subarctic region, this interaction is opposite that in the subtropical region (Figs. 8(a-k) and 9(a−d)), with a negative contribution from biological processes (i.e., a sink) throughout the year (Figs. 8(m−p) and 9(a−d)). Below 50 m depth, the advection term becomes more prominent in the subtropical region and in the Kuroshio Fig. 9 Monthly mean balances of DIC variation terms (mmol C m −3 day −1 ) along 165°E for 20°N, 30°N, 40°N , and 50°N at depths of 0, 50, 100, and 200 m. The monthly mean values were calculated within a range of five grid cells (22.0-45.5 km) from the target location. Colored lines indicate monthly means of the processes of DIC variation terms (ΔDIC) induced by advection, horizontal mixing, vertical mixing, biological processes, and airsea exchange process, and total DIC time variation (from Eq. 7). Positive and negative values indicate an increase and decrease in each DIC time variation term, respectively. The DIC variation term influenced by air-sea CO 2 exchange is considered only for the surface (a-d), not below it (e-p), where the green line representing zero is given as a reference for the other terms  15  20  25  30  35  40  45  50 15  20  25  30  35  40  45  50 15  20  25  30  35  40  45  50 15  20  25  30  35  40  45  50 Latitude (N) Latitude (N) Latitude (N) Latitude (N) Fig. 10 Vertical sections showing the relative contribution of biological processes to total DIC time variation along 165°E (from Eq. 7) for January, April, July, and October. Positive and negative percentages indicate the production and consumption of DIC induced by biological processes, respectively Extension region, connecting with total DIC time variation (Fig. 9(i-k, n, m)). Working against the small negative contribution of biological processes below 200 m depth in our model ( Fig. 9(n−p)) is advection, which can produce an increase or decrease in DIC in the subtropical region ( Fig. 9(n)). The modeled carbon cycle also includes latitudinal differences in total DIC time variation, as shown in Figs. 12 and 13(a). Differences in the contribution from biological processes are related to the total DIC time variation both horizontally and vertically (Figs. 12 and 13(a)). A distinctive area of DIC decrease exists near 0-200 m depth, and deepens southward from the subarctic region to the subtropical region (Figs. 11(a−d) and 12b). A comparison between the annual mean DIC time variation and the Chl-a maximum depth (Ishizu et al. 2019; Fig. 12) shows some similarity above 200 m depth. Advection dominates the trend in annual mean DIC time variation in the subtropical region south of 30°N below 200 m depth ( Fig. 9(n)). These results suggest that the positive and negative patterns of annual mean DIC time variation are caused by ocean currents (Fig. 9(n)).
The density range of this distinctive zone of contrast is comparable with that of the North Pacific Ocean Central Mode Water (CMW; σ θ = 26.0-26.5 kg m −3 ; Oka and Suga 2005; Oka and Qui 2012; Fig. 12b). Those waters are formed by winter ventilation around thermocline fronts, including the Kuroshio Extension front, the Kurhoshio Bifurcation front, and the subarctic front, and are then spread by advection (Oka and Suga 2005;Oka and Qui 2012). We therefore suggest that the positive and negative contrast in the temporal DIC variations depends on uptake in the ventilation areas north of 40°N and could be transported by advection (Fig. 12). A sensitivity experiment was performed in which the maximum depth of Chl-a was decreased by adjusting the biological parameters related to photosynthesis (I opt ), to check whether advection affects the contrast in DIC time variation below the surface layer.

Comparison with previous work on the marginal seas
The balance of the monthly mean DIC time variation terms in the East China and Japan seas ( Fig. 13(a, b)) shows that air-sea CO 2 exchange is generally offset by vertical mixing. These balances differ from the simulation results for the Yellow and East China Seas of Luo et al. (2015), which indicates that the contributions of vertical mixing, advection, and biological processes to the inorganic carbon cycle largely offset each other on the continental shelves and vary seasonally. The difference between the results of Luo et al. (2015) and those presented here may be due to missing processes in our model that may be required to represent local variability in these regions (e.g., the DIC input from river discharge), although the reproducibilities for the Japan and East China seas are high, with high correlation coefficients for several variables, including Chl-a, DIN, DIC, pH 25 , and Ω arg ( Table 2). The reproducibility for the Okhotsk Sea could not be evaluated because of a lack of observational data. Variations in DIC in Okhotsk Sea simulated by our model (Fig. 13(c)) are similar to those for the subarctic region of the Northwest Pacific, where the contributions of air-sea CO 2 exchange and vertical mixing are generally offset by each other (Fig. 13(c)) and biological processes make a subordinate contribution.
The relatively small contribution of biological processes in the DIC cycle in marginal seas can be explained by the latitude-dependent functions of biological parameters (Section 2), as we adjusted the biological parameters to focus mainly on the Northwest Pacific in the subtropical region and on subarctic regions (Ishizu et al. 2019;Eqs. 1-6). Tittensor et al. (2010) showed that biological communities in the marginal seas are dissimilar to those of the Pacific at the same latitude, suggesting that there is a limitation for the ability of the parameters optimizing method applied in the present model (Section 2.1). Yasunaka et al. (2013) discussed the inorganic carbon cycle using surface DIC climatological data, omitting physical processes such as advection and vertical mixing. Our results agree with those of Yasunaka et al. (2013) in terms of the importance of biological processes from spring to summer, especially in the subarctic region. However, in our model, biological processes are a subordinate contributor to variation in DIC compared with both vertical mixing and air-sea exchange throughout the year. The subsurface vertical structure of the biological contribution varies with latitude. Average daily net community production (NCP) values were estimated by Yasunaka et al. (2013) from March to July as > 14 mmol m −2 C day −1 in the Kuroshio Extension region (140-170°E, 30-40°N) and 2-6 mmol C m −2 day −1 in the subarctic region. The corresponding values of NCP in our model can be calculated as ∫ MLD 0 ∂ DIC ½ ∂t dz, where MLD indicates the mixed layer depth, which is defined as the depth at which the density is 0.125 kg m −3 greater than the density at the surface (Ohishi et al. 2019;Suga et al. 2004;Ohno et al. 2004). The modeled NCP (Fig. 14) shows positive values in the subtropical region of less than 2 mmol C m −2 day −1 and negative values in the subarctic region of less than − 4 mmol C m −2 day −1 , respectively (Fig. 14), which are smaller than their respective estimated values ). In addition, our model predicts that the distinctive positive/ negative (maximum/minimum) values spread sparsely and weakly only at the boundary between the Kuroshio Extension and the subarctic region. Such differences between our model and that of Yasunaka et al. (2013) can be explained by the horizontal and vertical advection/diffusion processes represented in our model, which should be considered for better understanding of the inorganic carbon cycle in our target regions. Fig. 13 The same as in Fig. 9 but for the monthly mean balance of each DIC term in the surface layer (at 0 m depth) in the East China Sea (30°N, 125°E), the Japan Sea (40°N, 135°E), and the Okhotsk Sea (55°N, 150°E )

Conclusions
To understand the seasonal inorganic carbon cycle in the Northwest Pacific, we performed simulations using a biogeochemical and carbon model coupled with an operational ocean model. The reproducibility of the model was sufficient to evaluate processes related to the seasonal inorganic carbon cycle. We found that contributions to the inorganic carbon cycle from air-sea CO 2 exchange generally offset those of vertical mixing at the surface in the Northwest Pacific. Biological processes are a subordinate contributor to variation in DIC and show a latitudinal dependence in the euphotic layer. Advection actively contributes to variation in DIC below the layer where biological processes contribute. A schematic representation of the main inorganic carbon cycle in the Northwest Pacific is shown in Fig. 15. DIC is absorbed from the atmosphere during winter south of 40°N and is released to the atmosphere north of 40°N (Fig. 15a). The DIC introduced in the subtropical region is conveyed from the surface to the subsurface by vertical mixing, and the DIC released in the subarctic region is compensated for by vertical mixing from subsurface layers (Fig. 15). In summer, the opposite pattern occurs because vertical mixing is weaker and the corresponding influence of mixing on variation in DIC decreases (Fig. 15b). Unlike the contribution from seasonal air-sea CO 2 exchange, the negative contribution induced by biological processes (i.e., a sink) below the surface occurs throughout the year (Fig. 15) and becomes stronger to the subarctic region. Below the surface, the contribution of advection is prominent in the subtropical and Kuroshio Extension regions. This contribution of advection connects with the increase/decrease in DIC in the subtropical region (Fig. 15).
Ocean circulation (Tally et al. 1993;Suga and Hanawa 1995a;Yasuda et al. 1997;Qui and Chen 2006) transports DIC-rich and DIC-poor water and redistributes them over time. The horizontal pattern of long-term DIC trends varies with depth. The impacts of long-term trends in ocean acidification indices (pH and aragonite saturation) also show horizontal and vertical dependencies. Future studies will use modeling experiments to evaluate decadal variations in the NPZD and carbon cycle. Results of these experiments are expected to improve our understanding of the variability in carbon, biological processes, and ocean acidification. . Arrows indicate the movements of inorganic carbon due to air-sea CO 2 exchange, to vertical mixing, to consumption of phytoplankton, and to advection Fig. 17 Horizontal distributions of the monthly DIC variation term induced by air-sea CO 2 exchange (a-d), advection (e-h), vertical mixing (i-l), biological processes (m-p), and the total DIC time variation term (q-t) at 100 m depth for January, April, July, and October, respectively Fig. 18 Horizontal distributions of the monthly DIC variation term induced by air-sea CO 2 exchange (a-d), advection (e-h), vertical mixing (i-l), biological processes (m-p), and the total DIC time variation term (q-t) at 200 m depth for January, April, July, and October, respectively We calculated air-sea CO 2 fluxes ( Fig. 8(a-d)) as ρ w V p k 0C P CO2 −P atm CO 2 (Kantha, 2004;Ishizu et al. 2019), where ρ w is the density of seawater (kg m −3 ); k 0C is the solubility of CO 2 in seawater (mol kg −1 atm −1 ; Weiss 1974); V p is the piston velocity, which depends on wind speed (U 10 ) and the Schmidt number (S c ) and is expressed in units of m s -1 (V p ¼ 8:61 Â 10 −7 U 2 10 S c =660 ð Þ −1=2 ; Wanninkhof, 1992); and P CO 2 and P atm CO2 are the partial pressures of CO 2 in surface waters and the atmosphere (μatm), respectively.
Compared with the global air-sea CO 2 fluxes evaluated on the basis of monthly climatological data (Takahashi et al. 2002(Takahashi et al. , 2009Yasunaka et al. 2013), our model provides a finer spatiotemporal resolution of variability, with both strong and weak contrasts in air-sea CO 2 fluxes in the target region (not shown). The seasonal changes in air-sea CO 2 flux simulated by our model (Fig. 8(a−d); note that the values in Fig. 8(a−d) do not represent the actual air-sea CO 2 flux, but the contribution of air-sea CO 2 exchange to the total DIC balance) are similar to the seasonal climatology (Takahashi et al. 2002;Yasunaka et al. 2013), although the annual mean air-sea CO 2 flux (Fig. 19) indicates a smaller sink in the Kuroshio Extension region between 30°N and 40°N, and the transition between sink and source areas is less distinct than that in the climatology (Takahashi et al. 2002;Yasunaka et al. 2013). For example, the climatology  shows an extensive sink area between 30°N and 40°N of > 6 mmol m −2 C day −1 , but our model predicts a mixture of sink and source areas in the range from − 5 to 5 mmol m −2 C day −1 . One reason for the relatively weak sink areas predicted by our model may be uncertainty in the calculation of pCO 2 , as revealed by the MATLAB program CO 2 sys.m (Orr and Epitalon 2015;Mackenzie et al. 2004). The calculation of pCO 2 performed here gives higher pCO 2 values in the ocean compared with the pCO 2 climatology . While et al. (2012) also referred to large estimation in the calculation of pCO 2 . 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 long as 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 the material. 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 directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.