Development of a Biogeochemical and Carbon Model Related to Ocean Acidification Indices with an Operational Ocean Model Product in the North Western Pacific

We developed a biogeochemical and carbon model (JCOPE_EC) coupled with an operational ocean model for the North Western Pacific. JCOPE_EC represents ocean acidification indices on the background of the risks due to ocean acidification and our model experiences. It is an off-line tracer model driven by a high-resolution regional ocean general circulation model (JCOPE2M). The results showed that the model adequately reproduced the general patterns in the observed data, including the seasonal variability of chlorophyll-a, dissolved inorganic nitrogen/phosphorus, dissolved inorganic carbon, and total alkalinity. We provide an overview of this system and the results of the model validation based on the available observed data. Sensitivity analysis using fixed values for temperature, salinity, dissolved inorganic carbon and total alkalinity helped us identify which variables contributed most to seasonal variations in the ocean acidification indices, pH and Ωarg. The seasonal variation in the pHinsitu was governed mainly by balances of the change in temperature and dissolved inorganic carbon. The seasonal increase in Ωarg from winter to summer was governed mainly by dissolved inorganic carbon levels.


Introduction
Ocean acidification poses a serious risk to marine organisms and ecosystems, including finfish and coral reefs in subtropical regions, and species or groups of organisms in polar regions [1][2][3][4][5]. The effects of ocean acidity have previously been reported and the effects of acidification are predicted to increase, with great risks to marine organisms [2,[4][5][6][7]. The global economic loss of organisms from ocean acidification has been estimated at $24 billion, $0.7 billion, $37 billion, $65 billion, and $30-375 billion for molluscs, echinoderms, crustaceans, finfish, and corals, respectively. Some organisms including molluscs and echinoderms may become locally extinct, and corals will be damaged by the combined effects of global warming and ocean acidification, which will reduce calcification, increase bio-erosion, and have synergistic effects [8]. All organisms except echinoderms could be seriously affected, and substantial economic losses are likely [8]. Therefore, these issues can no longer be ignored, and urgent action is necessary.  (Tables 2 and 3).

Model Formation
The JCOPE_EC was constructed on an offline tracer model driven by physical processes simulated by the model outputs of the JCOPE2M, and an operational eddy-resolving ocean general circulation model based on the Princeton Ocean Model, with a generalized sigma coordinate [38]. The JCOPE_EC is a three-dimensional high-resolution regional model that covers the western North Pacific (10.5°-62°N, 108°-180°E) at a horizontal resolution of 1/12° (9.1-4.4 km) and has 46 vertical active levels.

Equations Governing the Marine Ecosystem Model
A simple nitrogen-based plankton ecosystem model (NPZD) was used. It comprised five components including dissolved inorganic nitrogen (DIN, NO3), dissolved inorganic phosphate (DIP, PO4), phytoplankton (P), zooplankton (Z) and detritus (D) [26,27]. The carbon cycle model comprised three components: DIC, ALK and calcium carbonate (CaCO3), and was coupled with the NPZD model. Detritus corresponded to the concentrations of particulate organic nitrogen (PON). The pH (pHinsitu, pH25) and Ωarg values were calculated from temperature, salinity, DIC and ALK using a Fortran program based on the CO2sys MATLAB software [39,40]. The development of the nitrogen-based plankton ecosystem model followed the outlines described by Guo and Yanagi [27], Onitsuka and The colored red points indicate the JMA observation points. The red points in the shaded blue, orange, pink and green regions were used in the calculation of optimized operations and of correlation coefficients between observation and model outputs for the subarctic region, the Kuroshio extension, the subtropical region and the Japan Sea, respectively (Tables 2 and 3).

Model Formation
The JCOPE_EC was constructed on an offline tracer model driven by physical processes simulated by the model outputs of the JCOPE2M, and an operational eddy-resolving ocean general circulation model based on the Princeton Ocean Model, with a generalized sigma coordinate [38]. The JCOPE_EC is a three-dimensional high-resolution regional model that covers the western North Pacific (10.5 • -62 • N, 108 • -180 • E) at a horizontal resolution of 1/12 • (9.1-4.4 km) and has 46 vertical active levels.

Equations Governing the Marine Ecosystem Model
A simple nitrogen-based plankton ecosystem model (NPZD) was used. It comprised five components including dissolved inorganic nitrogen (DIN, NO 3 ), dissolved inorganic phosphate (DIP, PO 4 ), phytoplankton (P), zooplankton (Z) and detritus (D) [26,27]. The carbon cycle model comprised three components: DIC, ALK and calcium carbonate (CaCO 3 ), and was coupled with the NPZD model. Detritus corresponded to the concentrations of particulate organic nitrogen (PON). The pH (pH insitu , pH 25 ) and Ω arg values were calculated from temperature, salinity, DIC and ALK using a Fortran program based on the CO2sys MATLAB software [39,40]. The development of the nitrogen-based plankton ecosystem model followed the outlines described by Guo and Yanagi [27], Onitsuka and Yanagi [26] and Sasai et al. [25] (Figure 2a). The carbon cycle model is referred from Orr et al. [41], Ewen et al. [42] and Schmitnner et al. [18] (Figure 2b).
The source-sink terms for each individual biological tracer were defined by the following equations: dP dt = Phtosynthesis (DIN, DIP, P) − Respiration (P) − Mortality (P) − Grazing (P, Z), where W D is the sinking velocity of detritus, R P:N is the Redfield ratio of P (phosphorus):N (nitrogen) ( Table 1). The biogeochemical processes represented in Equations (1) where T is the seawater temperature, V max is the maximum growth rate from photosynthesis, and K DIN and K DIP are half saturation constants for DIN and DIP, respectively. We adopted optimal uptake kinetics for K DIN relating to V max , according to the theory of physiological acclimation of phytoplankton [25,43], defined as: V max · R P:N DIP/A (µmol · l −1 ) , where A is an affinity coefficient of basic cellular physiology representing the efficiency of nutrient encounters at the phytoplankton cell surface, R P:N is the stoichiometry of N to P [44,45], C P T is a temperature-dependent coefficient for photosynthesis, and I opt is the optimum light intensity. I(P;z) is the light intensity at depth z, and is determined by solar radiation at the sea surface (I s ) and the light extinction coefficient, k(P). These relationships are given by where k(P) = c dom + 0.054 P 0.6667 + 0.088 P, c dom is the light dissipation coefficient of seawater. k(P) depends on c dom and a shelf shading effect influenced by phytoplankton [21,27]. The respiration (extracellular excretion) of phytoplankton, and the mortality of phytoplankton and zooplankton were defined as follows: where C RP T , C MP T , and C MZ T are temperature-dependent coefficients for the respiration of phytoplankton, and the mortality of phytoplankton and zooplankton, respectively. The respiration of phytoplankton was assumed to be proportional to its concentrations and rate of photosynthesis. The mortality rates for phytoplankton and zooplankton were assumed to be proportional to the square of their concentrations (biomass). The grazing of phytoplankton by zooplankton was defined as a function of temperature and the concentrations of phytoplankton and zooplankton, expressed as where C GZ T is the temperature dependent-coefficient of grazing by zooplankton, G z is the maximum rate of grazing of phytoplankton at 0 • C, P* is the threshold phytoplankton concentration for possible grazing by zooplankton, and λ is the Ivlev constant. Zooplankton were assumed to graze only on phytoplankton.
Excretion and egestion by zooplankton were assumed to be proportional to their grazing and predating, as defined by where α z and β z are assimilation and growth efficiency constants for zooplankton, respectively. The decomposition of detritus is dependent on the water temperature and was expressed as follows: where V D is the decomposition rate at 0 • C, and C λ D T is the temperature coefficient of decomposition. Climatology damping conditions were adapted for the DIN and DIP source-sink terms, and were defined as follows: where γ is the inverse of the time scale (= 30 days), and DIN clim and DIP clim are the monthly climatology data for DIN and DIP, respectively. Some details of climatology damping are described in Section 3.1.

Equations Governing the Carbon Cycle Model
The modeling of the air-sea gas exchange and carbon chemistry followed the protocols from the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) [41,42]. Biological uptake and release occur in fixed elemental ratios for C, P and/or N. The production of DIC and ALK was controlled by changes in the inorganic nutrients (DIN, DIP) and calcium carbonate (CaCO 3 ) in molar relationships as follows: where R C:P and R ALK:N are the molar elemental ratios of C to P, and of ALK to N, respectively. Climatology damping conditions were also adapted for the DIC and ALK source-sink terms, as follows: where DIC clim and ALK clim are monthly climatology data for DIC and ALK, respectively (see Section 3.1). Calcium carbonate (CaCO 3 ) immediately sinks after its production. We did not treat the CaCO 3 concentration as a predictable variable in this model. The equations for approximating CaCO 3 were as follows: Pr (CaCO 3 ) = (Egsetion(Z) + Mortality(P) + Mortality(Z))R CaCO 3 /POC R C:P ), (22) where D CaCO3 is the CaCO 3 remineralization e-folding scale depth, and R CaCO3/POC is the ratio of CaCO 3 to the non-photosynthetic production of particulate organic carbon (POC). Pr (CaCO 3 ) is parameterized as a fixed ratio (R CaCO3/POC ) of the production of non-diazotrophic detritus [18].

Inorganic Carbon Cycling in Water
The chemical reactions in water involve the carbon-water-borate system [47].
where BOR can be approximated using a known parameter of salinity, S, and X = 1/[H + ] (note that pH = −log10[H + ]). These parameters were solved iteratively, as [ALK] approximates the carbonate alkalinity according to As γ c = [ALK]/[DIC], Equation (23) can be reduced to an iterative solution of the following two equations for two unknowns, X and γ c [47,48]: γ c is approximately 1.1, and this value can be used as a first guess. If X is known, the various ion concentrations can be defined [47] as follows: The apparent dissolution constants k 0 , k 1 , k 2 , k b , k w , and k arg are empirical functions of ambient temperature and salinity (see Mucci [49], Feely et al. [3] and other sources, cosys.m as summarized by Zeebe and Wolf-Gladrow [50]). The algorithm was checked using typical representative values provided by Machkenzie et al. [51], except for Ω arg .
The ALK flux at the sea surface was set to zero, while the DIC flux in the form of CO 2 gas was assumed to be proportional to the difference in the partial pressures of CO 2 in surface waters (PCO 2 ) and the atmosphere (P atm CO 2 ) [47,50], as described by where ρ w is the seawater density, and k 0c is the solubility of CO 2 in seawater in mol kg −1 atm −1 [52][53][54], and is expressed as follows: where T is in units of Kelvin, and 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 [53] as follows: S c is a function of the ambient temperature T (in • C), and, for CO 2, is expressed as follows:

Climatology and Forcing Data
We created monthly DIN and DIP climatological data on the model grid by interpolating the original data from the World Ocean Atlas 2013 (WOA13). Yasunaka et al. [55] created a monthly climatology for surface DIN and DIP in the North Pacific (10 • -60 • N, 120 • E-90 • W). We replaced the surface values for DIN and DIP in the WOA13 with the data from Yasunaka's study, in an area defined by 10.5 • -60 • N and 120 • -180 • E.
A climatological DIC database was created by combining an annual mean climatology product estimated using the method of Goyet et al. [56] with temperature and salinity data from the WOA climatology and the annual mean climatological datasets created by Key et al. [57]. There were no monthly climatology datasets available, except for the monthly surface climatology datasets of Yasunaka et al. [58]. We therefore replaced the surface values in the climatology database with Yasunaka's climatology data, for an area defined by 10.5 • -60 • N and 120 • -180 • E.
There were no sophisticated climatological datasets for ALK for the study region [46,56,57,59]. Consequently, we first combined two climatology datasets estimated by Goyet et al. [27] with temperature and salinity data from the WOA13 climatology and the annual mean climatology of Key et al. [57], and then replaced the surface data for depths above 200 m with the data estimated using Takatani's method [46].
We found that the first version of the created ALK data was significantly different from the ALK observed data obtained for the Japan Sea for the period 1997-2017, downloaded from the Japan Meteorological Agency (JMA) website, so we replaced the ALK data for the Japan Sea with the JMA observed data, on the assumption that the JMA data had the characteristics of ALK profiles in the Japan Sea.
After all the processes for creating the monthly climatological data for DIN, DIP, DIC and ALK were completed, we applied Gaussian smoothing with a spatial decorrelation scale (e-folding scale) of 100 km to remove any horizontal gaps in the data that may have formed during the repeated operations that combined the different data sources.
The model required forcing data to drive photosynthesis and CO 2 gas exchange in the surface/subsurface layer. Short wave radiation (I s ) was derived from the 6-hourly NCEP/NCAR (National Centers for Environmental Projection/National Center for Atmospheric Research) reanalysis data. Monthly mean values for the partial pressure of carbon dioxide (P atm CO 2 ) are available at the JMA website. Using data from the Ayasato Observatory in Japan, we created an empirical formula for CO 2 (in ppm) with the addition of a second-order polynomial (CO 2 poly ) and seasonal sinusoidal curve equations (CO 2 curv ), as follows: where t yr and t hr are the Julian year and the time in hours from 1 January each year, respectively. The left-hand terms in Equations (1)- (5), (17) and (18) denote transport-diffusion equations for tracers. The daily mean velocities, sea level, and turbulent diffusion coefficient from the JCOPE2M were included in the transport-diffusion equations. The daily mean temperature and salinity from the JCOPE2M were also used to evaluate the source-sink terms in Equations (1)-(5), (17), (18) and (21). Table 1 shows the biogeochemical parameters used in our model, and the relevant reference values used by Onitsuka and Yanagi [26]. We adjusted the biogeochemical parameters by comparing the model outputs and the observed data in specific target regions ( Figure 1, Table 2). The observed data were downloaded from the JMA website ( Figure 1). The phytoplankton concentrations were also validated with chlorophyll-a data downloaded from an File Transfer Protocol (FTP) site for the Moderate Resolution Imaging Spectroradiometer (MODIS)-Aqua Ocean Color Data.

Biogeochemical Parameters and Observational Data for Validation
In previous studies, where biogeochemical models were applied to the North Pacific and the marginal seas of Japan, the target area was assumed to be occupied by the same type of water mass and/or the same water circulation gyres in the subarctic area [19][20][21][22][23]25], the subtropical area [21,25], and/or the marginal seas [27,60] (Figure 1). In contrast, our model included both subarctic and subtropical areas and adjoining marginal seas. As biogeochemical situations differ depending on oceanic conditions [61][62][63], in our model we applied latitudinal differences for I opt and c dom , which were defined as follows: where I opt and c dom change latitudinally from 20 to 120 W/m −2 and from 0.015 to 0.045 m −1 , respectively; I 0 opt , I 1 opt , c 0 dom , and c 1 dom are parameters for adjusting from subtropical to subarctic regions (Table 1 and Figure 3); and Lat bnd and Lat slp (Lat slp_cdom ) are the coefficients for latitudinal boundaries and the latitudinal slopes for these parameters, respectively (Table 1 and Figure 3). The effects of latitudinal differences on I opt and c dom are described in Appendix A.
Using the Green's function method [36], we optimized the biogeochemical parameters Vmax, R, Mp, Gz, Mz, VPN, and WD (Equations (6)-(16) and Table 1). Chlorophyll-a data for 2015 were used in these optimization operations, and we adopted the optimized operations for the subarctic and subtropical regions separately. All the estimated biogeochemical parameters were found to be similar to those for the subarctic and subtropical regions ( Table 2), but we manually adjusted Vmax to 0.0492 day −1 , after adopting the latitudinal differences for cdom and Iopt (Appendix A).

Initial and Boundary Conditions
We used a 1-year period, 2015, for all the experiments, which were driven by the forcing data for that year. In a preliminary experiment, the phytoplankton were represented by chlorophyll-a concentrations from the World Ocean Atlas 2001. There was undesirable noise in the model outputs for phytoplankton in the preliminary experiment, so, for the initial conditions, i.e., the beginning of the first month in the numerical experiments, constant values for phytoplankton of 0.05 and 0.0 μg/L were set for depths above 150 m and deeper depths, respectively, and the model was then run continuously. The initial zooplankton concentrations were set at 10% of the phytoplankton concentration. The initial detritus concentration was set at 0.0 mmol/kg.
As a first step towards developing reliable operational systems for biogeochemical parameters related to ocean acidification indices, we adopted a diagnostic approach to obtain realistic  Using the Green's function method [36], we optimized the biogeochemical parameters V max , R, M p , G z , M z , V PN , and W D (Equations (6)-(16) and Table 1). Chlorophyll-a data for 2015 were used in these optimization operations, and we adopted the optimized operations for the subarctic and subtropical regions separately. All the estimated biogeochemical parameters were found to be similar to those for the subarctic and subtropical regions ( Table 2), but we manually adjusted V max to 0.0492 day −1 , after adopting the latitudinal differences for c dom and I opt (Appendix A).

Initial and Boundary Conditions
We used a 1-year period, 2015, for all the experiments, which were driven by the forcing data for that year. In a preliminary experiment, the phytoplankton were represented by chlorophyll-a concentrations from the World Ocean Atlas 2001. There was undesirable noise in the model outputs for phytoplankton in the preliminary experiment, so, for the initial conditions, i.e., the beginning of the first month in the numerical experiments, constant values for phytoplankton of 0.05 and 0.0 µg/L were set for depths above 150 m and deeper depths, respectively, and the model was then run continuously. The initial zooplankton concentrations were set at 10% of the phytoplankton concentration. The initial detritus concentration was set at 0.0 mmol/kg.
As a first step towards developing reliable operational systems for biogeochemical parameters related to ocean acidification indices, we adopted a diagnostic approach to obtain realistic distributions for biogeochemical variables. The variables DIN, DIP and DIC were initialized using the climatology data for these parameters at the beginning of the 1-year simulation (January), and ALK was initialized using the annual climatology at the beginning. The climatology data for DIN, DIP, DIC and ALK were applied at all open lateral boundaries, and were included in the damping terms in Equations (15), (16), (19) and (20), with a 30-day time scale. The DIN, DIP and DIC climatology data used for the open lateral boundaries and the damping conditions (Equations (15), (16), (19) and (20)) were replaced by monthly data at the beginning of each month. Preliminary results showed that the damping condition was effective in representing smooth and realistic variations in the biogeochemical variables (not shown). Using this approach, we were able to predict biogeochemical variables over the short-term (within one month), driven by the forcing including the ocean and atmosphere reanalysis data with constraints on the climatology of DIN, DIP, DIC and ALK. All the experiments in this study were driven by the forcing created from the ocean (JCOPE2M) and atmospheric (NCEP/NCAR) reanalysis data.

Comparison of the Observed Data with the Simulated Outputs
The monthly mean distributions of satellite and modeled chlorophyll-a values at the surface are shown in showed a relatively larger decrease in DIC in summer (from July to September; Figure 5f) compared with the observed data for 40° N (Figure 5e). Due to a lack of observed data, we could not compare the seasonal variations in ALK [13,26,52,53]. Instead of using the ALK climatology, we plotted the JMA observed data (Figure 8a-c). These data showed that the ALK concentrations were relatively low (2200-2250 μmol/kg) in the subarctic region above 40°N and in the Okhotsk Sea, and were higher (2250-2300 μmol/kg) below 40°N from spring to autumn (Figure 8a   The model gave adequate simulations of DIN and reproduced the basic features and seasonal variability, when the monthly climatology data were dampened (Figures 5c,d and 6e-h). The climatology indicated large differences in the DIN concentrations between the subarctic and subtropical regions throughout the year (Figure 6a-d). The DIN values ranged from more than 4 µmol/kg above 45 • N in the subarctic region to almost 0 µmol/kg in subtropical regions at various times of the year (Figure 6a-d).
The modelled concentrations of DIN in the subtropical and subarctic regions were reasonably consistent with observed data. The seasonal variability in the climatology for the subarctic region above 40 • N (Figures 5c and 6a-d) was reproduced by the model (Figures 5d and 6e-h). The modeled DIP distributions also reflected the observed data (not shown).  The concentrations of DIC (Figure 7a-d) were high (>2120 µmol/kg) in the Okhotsk Sea and the northwestern North Pacific along the Kamchatka Peninsula and the Kuril Islands in winter (January) and spring (April) (Figure 7a,b). The concentrations were less than 2040 µmol/kg in the Kuroshio Extension region and in the subtropical region below 35 • N (Figure 7a-d). There was noticeable seasonal variation in the concentrations (Figure 7a-d). For the Okhotsk Sea and northwestern North Pacific area mentioned above, DIN was not detected in the summer (July) or autumn (October) (Figure 7a-d).
The model simulations were similar (Figure 7e-h), with a minor difference in the continental area of the Okhotsk Sea in summer and autumn (Figure 7c,d,g,h). The model reproduced the observed temporal variability in the surface DIC values for the latitudes in subtropical and subarctic regions along longitude 165 • E (Figure 5e,f). The time-series in simulated DIC for 40 • N showed a relatively larger decrease in DIC in summer (from July to September; Figure 5f) compared with the observed data for 40 • N (Figure 5e).
Due to a lack of observed data, we could not compare the seasonal variations in ALK [13,26,52,53]. Instead of using the ALK climatology, we plotted the JMA observed data (Figure 8a-c). These data showed that the ALK concentrations were relatively low (2200-2250 µmol/kg) in the subarctic region above 40 • N and in the Okhotsk Sea, and were higher (2250-2300 µmol/kg) below 40 • N from spring to autumn (Figure 8a [55], and from the model outputs, respectively.  [55], and from the model outputs, respectively.   Figure 9 shows the surface pHinsitu (pH25) and Ωarg values calculated from the model outputs for temperature, salinity, DIC and ALK (Figure 9a-d). The pHinsitu ranged from 7.85 to 8.10 in the northwestern Pacific area (Figure 9a-d), except in winter and spring (January and April) for the East China Sea near the coast of China (Figure 9a,b). The values in the side of Pacific basically ranged from 8.00 to 8.05. The values were lower (7.85-7.90) in the Okhotsk Sea and the Japan Sea throughout the year, and on the Pacific Ocean side along the Kamchatka Peninsula and Kuril Islands from winter to   Figure 9 shows the surface pHinsitu (pH25) and Ωarg values calculated from the model outputs for temperature, salinity, DIC and ALK (Figure 9a-d). The pHinsitu ranged from 7.85 to 8.10 in the northwestern Pacific area (Figure 9a-d), except in winter and spring (January and April) for the East China Sea near the coast of China (Figure 9a,b). The values in the side of Pacific basically ranged from 8.00 to 8.05. The values were lower (7.85-7.90) in the Okhotsk Sea and the Japan Sea throughout the year, and on the Pacific Ocean side along the Kamchatka Peninsula and Kuril Islands from winter to  Figure 9 shows the surface pH insitu (pH 25 ) and Ω arg values calculated from the model outputs for temperature, salinity, DIC and ALK (Figure 9a-d). The pH insitu ranged from 7.85 to 8.10 in the northwestern Pacific area (Figure 9a-d), except in winter and spring (January and April) for the East China Sea near the coast of China (Figure 9a,b). The values in the side of Pacific basically ranged from 8.00 to 8.05. The values were lower (7.85-7.90) in the Okhotsk Sea and the Japan Sea throughout the year, and on the Pacific Ocean side along the Kamchatka Peninsula and Kuril Islands from winter to spring (Figure 9a N (Figure 10a).

Ocean Acidification Indices pH and Ω arg
The pH 25 values were lower in the northern region and higher in the southern region (Figure 9e-h). The lowest values (7.5-7.6) occurred in the Okhotsk Sea and along the Kamchatka Peninsula on the Pacific Ocean side in winter (January) and spring (April) (Figure 9e,f). The highest values (8.10-8.15) were found in the subtropical region at latitudes below 20 • N. The pH 25 values increased in summer (Figures 9e-h and 10b). The amplitudes were larger in the northern region than in the southern region. The lowest and highest values at each latitude occurred in April and August, respectively.
The surface Ω arg values (>3.5) were high in the subtropical region, and low (2.0-3.0) in the subarctic region (Figure 9i-l). The values were lowest (1.0-1.5) in the Okhotsk Sea and along the Kamchatka Peninsula and the Kuril Islands on the Pacific Ocean side in winter (January) and spring (April). The values increased in summer (July), with similar spatial variability in autumn (October). The time-series at longitude 165 • E showed that the Ω arg values were highest in August-September and lowest in April (Figure 10c), and the seasonal variation was similar to that for pH 25 (Figure 10b).
Takahashi et al. [59] and Jiang et al. [64] created a global map of the surface climatology for pH insitu and Ω arg in spring/summer and autumn/winter, but not for the Okhotsk, Japan and East China Seas. Their climatology showed that the pH insitu and the Ω arg varied seasonally from 0.06 to 0.08 and from 0.50 to 1.00, respectively, in the subtropical, Kuroshio Extension, and were similar to the modeled values (Figure 9a

Reproducibility of the Biogeochemical Variables
The correlation coefficients between the JMA observed data and the model outputs for each region for 2015 were calculated for chlorophyll-a, DIN, DIC, pH 25 and ALK ( Table 3). The observed data and model outputs could be compared as they overlapped in time and location. The natural logarithm values of chlorophyll-a were calculated, as the plankton concentrations were not normally distributed but followed a lognormal distribution [66].
The  (Table 2). Scatter plots (Figures 11-13) showed that the parameters were highly correlated with each other (not shown for DIP, DIC and ALK).

Sensitivity Experiment for pH and Ωarg Based on Constant Temperature, Salinity, DIC and ALK
To assess which variables were important in determining the ocean acidification indices pH and Ωarg, we performed a series of sensitivity experiments assuming a constant set of conditions for temperature (= 21 C), salinity (= 34.44), DIC (= 1995.2 μmol/kg) and ALK (= 2260.9 μmol/kg). In the sensitivity experiments, pH and Ωarg were calculated with the assumed fixed variables of temperature, salinity, DIC and ALK, individually, in place of these original target variables of the model outputs. The other variables, except the assumed constant variables, were the same as those

Sensitivity Experiment for pH and Ω arg Based on Constant Temperature, Salinity, DIC and ALK
To assess which variables were important in determining the ocean acidification indices pH and Ω arg , we performed a series of sensitivity experiments assuming a constant set of conditions for temperature (= 21 • C), salinity (= 34.44), DIC (= 1995.2 µmol/kg) and ALK (= 2260.9 µmol/kg). In the sensitivity experiments, pH and Ω arg were calculated with the assumed fixed variables of temperature, salinity, DIC and ALK, individually, in place of these original target variables of the model outputs. The other variables, except the assumed constant variables, were the same as those of the model outputs (Figures 4e-h, 6e-h, 7e-h and 8; Section 4). The assumed constant variables were determined by the values corresponding to median values calculated from the model outputs in the region of 10.5 • -62 • N and 108 • -180 • E in the surface in January (e.g., Figures 6e and 7e). Figure 14 shows the horizontal distributions of the surface pH difference (∆pH), which was calculated by subtracting the model output for each sensitivity experiment using each fixed temperature, salinity, DIC and ALK value from the basic model outputs shown in Figure 9a-d. Measuring the differences enabled the individual effects of each target variable to be assessed. The temperature and DIC had most influence on the pH insitu in winter, spring and autumn, while ALK had a considerable influence in summer.
The variables that had most influence on Ω arg throughout the year were DIC and ALK ( Figure 15). The sensitivity experiments using the fixed values for temperature, DIC and ALK ( Figure 15) showed positive effects in the subarctic region and negative effects in the subtropical region. The effects of temperature and salinity, while detectable, were small compared with the effects of DIC and ALK (Figure 15e-h).

Effect of Seasonal Processes on pH insitu
The sensitivity analysis of the effect of temperature on pH showed a clear positive and negative contrast between the northern and southern regions, respectively (Figure 14a-d). This was related to the dissociation equilibrium: H 2 O↔H + +OH − . A relatively low (high) temperature in the northern (southern) region resulted in an increase (decrease) in the pH by shifting the equilibrium in the direction: H 2 O←H + +OH − (H 2 O→H + +OH − ). With the transition from winter to summer, the pH decreased throughout the entire region as a result of increasing temperature (Figure 14a-c), while the pH increased with cooling in autumn in some parts of the northern region (Figure 14d).
The horizontal distribution of surface DIC (Figure 7e-h) led to a significant negative and positive contrast between the northern and southern regions, respectively (Figure 14i-l). This was associated with the CO 2 equilibrium in sea water, defined by: CO 2 +H 2 O↔H + +HCO 3 − . Relatively high (low) DIC concentrations in the northern (southern) region tended to result in lower (higher) pH values by shifting the equilibrium in the direction: in DIC from winter to summer (Figure 7e-g) caused the pH to increase throughout the entire region (Figure 14j,k), by shifting the equilibrium in the direction: The relatively high concentration of ALK in the subtropical gyre (Figure 8a-d) resulted in a relatively high pH (Figure 14m-p). As the ALK decreased from winter to summer (Figure 8a-c) the pH decreased throughout the entire region (Figure 14m-o). This trend was more evident in the northern region than in the southern region (Figure 14o), reflecting the relatively enhanced decrease in ALK in the northern region (Figure 8c).
The sensitivity of salinity to pH was not significant compared with the other variables, and the seasonal variability in pH caused by the variation in salinity was also minor (Figure 14e-h). Lower salinity levels led to higher pH values, and was locally evident around the mouth of the Changjiang River.
From the sensitivity analyses, we inferred that overall seasonal variations in the pH insitu (Figure 9a-d and Figure 10a) were largely governed by the balance between temperature and DIC, and to a lesser extent by ALK. In the southern region, the relationship with temperature was relatively important, and thus the pH insitu decreased from winter to summer. In contrast, the pH insitu increased slightly from winter to summer in the northern region, where the relationship with DIC and ALK was slightly stronger. The pH 25 variation mainly reflected the effects of DIC and ALK variations (Figure 9e-h and Figure 10b). However, the model probably overestimated a decrease in ALK in summer, as shown in Figure 8b,f, and the variations in the pH insitu may actually be mainly governed by temperature and DIC.

Seasonal Processes Affecting Ωarg Variations
The sensitivity of the target variables to Ωarg are shown in Figure 15. The temperature dependence of Ωarg can be explained by the characteristics of the aragonite solubility product, karg (arg = ([Ca 2+ ][CO3 2− ])/karg). The value of karg decreases as the water temperature increases [34]. Latitudinal differences in the sensitivity of temperature resulted in higher and lower Ωarg values in the southern and northern regions (Figure 15a-d), respectively. Seasonal warming/cooling affected the seasonal variation in Ωarg throughout the entire region.
The reaction for the dissolution of marine carbonate [3,6], CO2+H2O+CO3 − ↔2HCO3 − , was related to latitudinal differences in the sensitivity of DIC (Figure 15i-l). As the DIC concentration was relatively high in the northern region (Figure 7e-h), the concentration of carbonate ions (CO3 − ) was relatively low, and consequently the Ωarg value was relatively low in this region (Figure 15i-l). The seasonal decrease in DIC from winter to summer caused the Ωarg values to increase throughout the entire region (Figure 15k). The sensitivity of ALK to Ωarg showed that higher ALK levels led to higher Ωarg values because the higher ALK levels were related to higher carbonate ion concentrations. With

Seasonal Processes Affecting Ω arg Variations
The sensitivity of the target variables to Ω arg are shown in Figure 15. The temperature dependence of Ω arg can be explained by the characteristics of the aragonite solubility product, k arg (Ω arg = ([Ca 2+ ][CO 3 2− ])/k arg ). The value of k arg decreases as the water temperature increases [34]. Latitudinal differences in the sensitivity of temperature resulted in higher and lower Ω arg values in the southern and northern regions (Figure 15a-d), respectively. Seasonal warming/cooling affected the seasonal variation in Ω arg throughout the entire region. The reaction for the dissolution of marine carbonate [3,6], CO 2 +H 2 O+CO 3 − ↔2HCO 3 − , was related to latitudinal differences in the sensitivity of DIC (Figure 15i-l). As the DIC concentration was relatively high in the northern region (Figure 7e-h), the concentration of carbonate ions (CO 3 − ) was relatively low, and consequently the Ω arg value was relatively low in this region (Figure 15i-l). The seasonal decrease in DIC from winter to summer caused the Ω arg values to increase throughout the entire region ( Figure 15k). The sensitivity of ALK to Ω arg showed that higher ALK levels led to higher Ω arg values because the higher ALK levels were related to higher carbonate ion concentrations. With a decrease in ALK in the northern region from winter to summer (Figure 8a-c), the Ω arg values tended to decrease (Figure 15m-o). The sensitivity of salinity to Ω arg (Figure 15e-h) was low, and similar to that of the pH insitu (Figure 14e-h).
From the sensitivity analysis, we inferred that DIC and ALK had most impact on the seasonal variability in Ω arg across the entire region throughout the year, with less of a temperature effect. However, the summer increase in Ω arg (Figure 10c a decrease in ALK in the northern region from winter to summer (Figure 8a-c), the Ωarg values tended to decrease (Figure 15m-o). The sensitivity of salinity to Ωarg (Figure 15e-h) was low, and similar to that of the pHinsitu (Figure 14e-h).
From the sensitivity analysis, we inferred that DIC and ALK had most impact on the seasonal variability in Ωarg across the entire region throughout the year, with less of a temperature effect. However, the summer increase in Ωarg (Figure 10c

Conclusions
We developed a new biogeochemical model (JCOPE_EC) for carbon processes, based on the physical background of a numerical model product, JCOPE2M, which is a three-dimensional operational eddy-resolving model product. Comparison with observed data for 2015 showed that the JCOPE_EC model adequately reproduced the basic features of chlorophyll-a, DIN (DIP), DIC and

Conclusions
We developed a new biogeochemical model (JCOPE_EC) for carbon processes, based on the physical background of a numerical model product, JCOPE2M, which is a three-dimensional operational eddy-resolving model product. Comparison with observed data for 2015 showed that the JCOPE_EC model adequately reproduced the basic features of chlorophyll-a, DIN (DIP), DIC and ALK. The seasonal variability in these biogeochemical variables was similar to the observed (climatological) variability, although the modelled seasonal variability in chlorophyll-a and ALK deviated somewhat from the observed.
The overestimation of the increase in chlorophyll-a in the subarctic region (Figures 4 and 5b) may be related to iron cycling, which was not included in JCOPE_EC [67]. Iron cycling may also explain the overestimates of biological production of CaCO 3 and also ALK reduction (Figures 5g and  8f). The present model will need to be improved to account for these processes.
The JCOPE_EC also represented the ocean acidification indices pH and Ω arg , based on the model outputs for temperature, salinity, DIC and ALK. The simulated values were consistent with observed data, although the pH insitu and Ω arg values in some parts in the Okhotsk Sea and the Japan Sea from winter to spring were lower than expected (Figure 9a-d,i-l). This is probably a consequence of the uncertainty in the climatology data for these parameters, and should be investigated more fully.
Sensitivity experiments, aimed at assessing which variables had most effect on the ocean acidification indices, were performed using fixed values for temperature, salinity, DIC and ALK. The results showed that seasonal variations in pH insitu were largely governed by the corresponding balances between temperature and DIC, although ALK had some influence in summer in model. DIC was the main driver of the increases in the Ω arg from winter to summer, with some lesser effects of ALK and temperature. Salinity had little effect on pH and Ω arg , except around the Changjiang River mouth.
At present, it is difficult to operationally acquire nowcast/forecast information on biogeochemical variables at the basin scale. However, the JCOPE_EC was designed to use operational physical data from the JCOPE2M. Thus, with further improvements in its design, the JCOPE_EC could become an effective tool for providing realistic information about biogeochemical variables in the studied marine area.
The daily snapshots of the distributions of chlorophyll-a, DIN, pH, DIC and ALK in spring 2015 in Figures 16 and 17 show that the JCOPE_EC was able to represent the distribution of biogeochemical factors associated with the meandering of the Kuroshio Current and other eddy phenomena. The Kuroshio Current transports water that is low in nutrients, DIC and ALK to downstream regions. Prominent areas of enhanced chlorophyll-a concentrations were evident, and their occurrence depended on ocean currents and the distribution of nutrient-rich waters ( Figure 16). These mesoscale phenomena involve complex factors affecting DIC and ALK, and so affect the associated pH and Ω arg ( Figure 17). Therefore, the JCOPE_EC has the potential to provide current realistic information on ocean acidification indices, which vary in space and time because of various oceanic phenomena, although validation and improvements in the model will be necessary.
By using this operational model (JCOPE_EC), we have already launched a trial to perform a nowcast/forecast experiment targeting ocean acidification indices (pH insitu and Ω arg ), which are reproduced with the nowcast/forecast operational ocean model products [35,36,38]. A website (https: //www.marinecrisiswatch.jp/mcwatch/prediction/jcope/index.html) experimentally shows up-to-date results of the nowcast/forecast system to demonstrate our studies effectively.
The current version of the JCOPE_EC system includes climatology damping conditions for some biogeochemical variables, which highlights that the present model is diagnostic. To investigate the mechanisms that drive the modeled ecosystem and the carbon cycle, we will develop a more prognostic model without climatology damping in future studies, based on the present model. We also note that there is an extreme lack of publicly available observed biogeochemical data for the Okhotsk Sea and the East China Sea. We hope that the availability of observed data will improve, as, with more data, we can reduce the uncertainty in future model outputs for these regions.    Monthly climatology data for DIN, DIP, and DIC from Yasunaka et al. [24,60] were downloaded from: http://soop. jp/index.html. In-situ observed data obtained from the Japan Meteorological Agency (JMA) were downloaded from https://www.data.jma.go.jp/kaiyou/db/vessel_obs/data-report/html/ship/ship.php. Partial pressure data for carbon dioxide obtained at the Ayasato observatory were downloaded from https://ds.data.jma.go.jp/ghg/kanshi/ obs/co2_monthave_ryo.html. Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua Ocean Color data for 2015 were downloaded from the Physical Oceanography Distributed Active Archive Center (PODAAC) at ftp://podaac-ftp.jpl.nasa.gov/allData/modis/L3/aqua/4um/v2014.0/4km/daily/2015/.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Preliminary results showed that the chlorophyll-a maximum depth ( Figure A1b) reproduced using constant coefficients for c dom and I opt were shallower than the observed depths ( Figure A1a). The chlorophyll-a reached a maximum at 0-50 m deep in the subarctic region of the northwestern Pacific, and this gradually deepened with distance southwards, and reached a maximum at 100-150 m deep in the subtropical region in summer ( Figure A1a). The depths estimated by the preliminary model at which chlorophyll-a reached a maximum in the subtropical region (10 • -20 • N) were much shallower (above 50 m) for all latitudes ( Figure A1b).
To reduce the bias in the preliminary model results for the chlorophyll-a maximum depth, we adopted values for c dom and I opt that varied by latitude (Equations (39) and (40)) (Figure 3a,b). Using this approach, the model reproduced the vertical distribution and maximum depth of chlorophyll-a very effectively for the subtropical region, although the chlorophyll-a values ( Figure A1c-e) were higher than the observation values ( Figure A1a,e). As the parameter dependences of c dom and I opt were nonlinear, we manually adjusted the parameter values after comparing the modeled and observational results.
Previous NPZD models of the marginal seas of Japan have used various constant I opt values. For example, Kawamiya et al. [19][20][21] used an I opt value of 48.8 W/m 2 for station Papa (50.1 • N, 144.9 • W) [19,20] and the North Pacific [21], Yamanaka et al. [48] used an I opt value of 104.7 W/m 2 for the North Western Pacific, Sasai et al. [25] used an I opt value of 100 W/m 2 for stations S1 and K2, Kishi et al. [22] used an I opt value of 104.7 W/m 2 for the subarctic region of the North Pacific, Guo and Yanagi et al. [27] used an I opt value of 150.0 W/m 2 for the East China Sea, and Onitsuka and Yanagi [26] used an I opt value of 70.0 W/m 2 .
Edwards et al. [61] characterized growth-irradiance relationships from information about 308 growth-irradiance experiments performed on 119 species of marine phytoplankton sampled from all major functional groups, and reported that the I opt values did not vary among taxonomic groups. However, all traits varied considerably within most groups, and optic-ocean isolates tended to have lower I opt values than coastal isolates. This implies that phytoplankton may change their physiological characteristics depending on the light environment. Their study appears to be consistent with the latitudinal difference we found for I opt .
Various studies have used constant values ranging from 0.035 to 0.04 m −1 for c dom in the North Pacific [19][20][21][22]25,48]. These c dom values represent light dissipation coefficients, which are linked to the distribution of chromophoric dissolved organic matter. Nelson and Sigel [63] showed that the absorption of surface chromophoric dissolved organic matter plus detrital particle at 443 nm was 0.005 m −1 in the subtropical region of the North Pacific but gradually increased to 0.05 m −1 in the subarctic region. This indicates that the distribution of chromophoric dissolved organic matter make more invisible the subsurface layer in the subarctic region than in the subtropical region. The latitudinal difference for c dom used in the present study represents a change from 0.015 m −1 in the subtropical region to 0.045 m −1 in the subarctic region (Figure 3a), which is consistent with observed distributions of chromophoric dissolved organic matter [63].