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

Abstract

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 pH_insitu 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.

1. 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.

The levels of ocean acidification are evaluated using parameters such as the pH and/or the aragonite saturation (Ω_arg). However, there are limited observed data for pH and Ω_arg, so it is difficult to assess the current status of ocean acidification both locally and globally [3,9]. This situation is unlikely to change even if new automatic monitoring systems, such as biological profiling floats [10], are established and operational monitoring, using volunteer observing ship (VOS) pCO₂ transects [11,12], increases.

Numerical model simulations are important tools that support the interpretation of observed data at the required spatial and temporal scales in a dynamically consistent manner, and are used to characterize the underlying physical–biological processes [13,14,15,16,17]. The term ‘modeling’ describes various processes including those associated with nitrate, phytoplankton, zooplankton, and detritus (NPZD processes), which are connected to the carbon cycle throughout photosynthesis, remineralization, and calcination [18]. These NPZD processes influence ocean acidification indices, so marine carbon and biogeochemical models are useful for understanding past, present, and future ocean acidification.

Biogeochemical models focused on NPZD processes in the North Pacific have been developed by Kawamiya et al. [19,20,21], Kishi et al. [22], Fujii et al. [23], Yoshie et al. [24] and Sasai et al. [25]. Onitsuka and Yanagi [26], Guo and Yanagi [27] and Yoshie et al. [24] have developed biogeochemical models for marginal seas including the Japan Sea, the East China Sea and the Seto Inland Sea, respectively. These models included primary productivity among the NPZD processes, but did not include the carbon cycle. Zhao and Guo [28] established an inorganic carbon system module for the Yellow and East China Seas, and Luo et al. [29] constructed a carbon biogeochemical model using the NORWeigian Ecological Model (NORWECOM) and a hydrodynamic module [28].

Earth System Models (ESMs) are climate models that may include physical processes and biogeochemical cycles, and so can facilitate representation of anthropogenic effects on the carbon cycle. In modeling future scenarios, marine ecosystems in ESMs are adapted to particular environmental conditions based on social and economic scenarios (Coupled Model Intercomparison Project Phase 5; [30,31,32,33,34]). Several modeling studies have used coupled climate–marine biogeochemical models to investigate how environmental stressors might evolve under climate change scenarios [30,34]. While ESMs provide past/present/future distributions of several biogeochemical components, and approximately reproduce the large-scale patterns of the subpolar and subtropical gyres, their outputs are not necessarily adequate for reproducing smaller-scale patterns, including shelf processes and/or mesoscale phenomena.

The Japan Coastal Ocean Predictability Experiment (JCOPE: http://www.jamstec.go.jp/jcope/) is an operational eddy-resolving ocean physical model. The latest version of the model (JCOPE2M) [35,36] assimilates satellite sea surface height anomalies and sea surface temperature, and temperature and salinity data obtained from floating profilers and/or ship measurements. The JCOPE2M covers the marginal seas of Japan (including the Okhotsk, Japan and East China Seas) in the geographic range defined by 10.5°–62° N and 108°–180° E (Figure 1). It reproduces the Kuroshio, Oyashio, mesoscale eddies, and upwelling and mixing, which play crucial roles in biological processes related to primary productivity in the upper ocean. The model products provide reanalysis data based on assimilations of all the available observed data. The model outputs include data on velocities, temperature, salinity, and sea elevation that are corrected using data assimilation techniques, and the model more realistically reproduces the physical conditions than non-operational physical models.

To improve our insights into the drivers of ocean acidification, we aimed to develop the JCOPE further to incorporate a new marine ecosystem and carbon cycle model (hereafter JCOPE_EC), driven by the physical processes represented by the JCOPE2M model outputs [37]. In this contribution, we provide a description of the JCOPE_EC, and report on the quality of the model outputs by comparing them with available observed data. The present version of the JCOPE_EC includes damping terms for the climatological states of dissolved inorganic nitrate (DIN), dissolved inorganic phosphate (DIP), dissolved inorganic carbon (DIC) and total alkalinity (ALK) and so can be used to determine the distributions of ocean acidification indices in the North Western Pacific, with a constraint to the climatological information. Sensitivity analysis was performed using fixed values for temperature, salinity, DIC, and ALK, to assess which variables were important in determining the ocean acidification indices, pH and aragonite saturation (Ω_arg).

This is the first time the development of an operational carbon cycle model with eddy-resolving horizontal resolution focused on the North Western Pacific has been reported. Some fundamental data required for our study (e.g., climatology data) had to be acquired because only a few of the necessary data for research in this challenging scientific field were available.

Details of the model formulation are provided in Section 2. The data and model parameters used in the study are described in Section 3. Section 4 describes how the model was validated with observed data. In Section 5, the influence of some environmental variables on ocean acidification indices (pH and Ω_arg) are discussed. The main findings of the study are overviewed, and future issues are discussed, in Section 6.

2. 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.

2.1. 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₃), dissolved inorganic phosphate (DIP, PO₄), phytoplankton (P), zooplankton (Z) and detritus (D) [26,27]. The carbon cycle model comprised three components: DIC, ALK and calcium carbonate (CaCO₃), and was coupled with the NPZD model. Detritus corresponded to the concentrations of particulate organic nitrogen (PON). The pH (pH_insitu, pH₂₅) 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:

$$\frac{dP}{dt}=Phtosynthesis (DIN,DIP,P)-\mathrm{Re}spiration (P)-Mortality (P)-Grazing (P,Z),$$

(1)

$$\frac{dZ}{dt}=Grazing (P,Z)-Mortality (Z)-Egestion (Z)-Excertion (Z),$$

(2)

$$\frac{dD}{dt}=Egestion (Z)+Mortality (P)+Mortality (Z)-Decomposition (D)-\frac{\partial }{\partial z}\left(W_{D}D\right)\hspace{1em},$$

(3)

$$\begin{array}{ll}\frac{d(DIN)}{dt}=& -Phtosynthesis (DIN,DIP,P)+Respiration (P)+Excertion (Z)\\ & +Decomposition (D)-\mathit{Clim}\_Damp(DIN)\end{array}$$

(4)

$$\frac{d(DIP)}{dt}=\frac{1}{R_{P:N}}\frac{d(DIN)}{dt}-\mathit{Clim}\_Damp(DIP),$$

(5)

where W _D is the sinking velocity of detritus, R_P:N is the Redfield ratio of P (phosphorus):N (nitrogen) (

Table 1. Biogeochemical parameters used in the JCOPE_EC. The biogeochemical values used by Onitsuka and Yanagi [26] are also provided as a reference.

Symbol	Definition	Value	Units	Reference Values Onitsuka and Yanagi [26]
Ecosystem Model
For phytoplankton
Vmax	Growth rate for phytoplankton	0.0492	day−1	0.6
A	Affinity coefficient of basic cellular physiology	6.75	mmolN−1m−1day−1	−
KDIN	Half saturation constant for DIN	1.50	mmolN m−1	1.5
KDIP	Half saturation constant for DIP	0.0940	mmolN m−1	0.09
MP	Phytoplankton mortality rate at 0 °C	0.0400	(mmolNm−1) m−3 day−1	0.07
Pmin	Threshold of phytoplankton mortality	0.0587	(mmolNm−1) m−3	−
R	Phytoplankton respiration rate at 0 °C	0.0317	day−1	0.3
CPT	Temperature coefficient for photosynthesis	0.0392	°C−1	0.0693
CRPT	Temperature coefficient for phytoplankton respiration	0.0519	°C−1	0.0519
CMPT	Temperature coefficient for phytoplankton mortality	0.0693	°C−1	0.0693
Iopt	Optimum light intensity for phytoplankton	* 20–120	W m−2	70
Iopt0		100.0		−
Iopt1		70.0		−
latbnd	Boundary for latitudinal differences	45.0		−
latslp	Slope for latitudinal differences	4.0		−
cdom	Light dissipation coefficient of seawater	* 0.015–0.045	m−1	0.04
cdom0		0.0275	m−1	−
cdom1		0.0250	m−1	−
Latslp_dom	Slope for latitudinal differences for cdom	1.5		−
For zooplankton
GZ	Maximum grazing rate of zooplankton at 0 °C	0.05	day−1	0.3
Λ	Ivlev constant	1.4	(mmolN m−3)−1	1.4
MZ	Zooplankton mortality rate at 0 °C	0.05	°C−1	0.07
βz	Growth efficiency of zooplankton	0.3		0.3
αz	Assimilation efficiency of zooplankton	0.7		0.7
P*	Zooplankton threshold value for grazing on phytoplankton	0.0430	(mmolNm−1) m−3	0.043
CGZT	Temperature coefficient for zooplankton grazing	0.0390	°C−1	0.0693
CMPT	Temperature coefficient for zooplankton mortality	0.0693	°C−1	0.0693
For diatoms
WD	Singing velocity of detritus	6.7	m day−1	10
VPN	Decomposition rate at 0 °C (DET→DIN)	0.050	day−1	0.05
CλDT	Temperature coefficient for decomposition	0.0693	°C−1	0.0693
Carbon Cycle Model
RP:N	Stoichiometry of nitrogen to phosphorus [46]	16.0		−
RC:P	Molar elemental ratios	112.0		−
RCaCO3/POC	CaCO3 over nonphotosynthetical POC production ratio	0.035		−
RALK:N	Alkalinity over non-photosynthetic N production ratio	0.001		−
DCaCO3	CaCO3 remineralization e-folding depth	3500.0	m	−

* indicates that latitudinal differences are the biogeochemical parameters adopted latitudinal differences. The underlines indicate the model types used the setting parameters.

). The biogeochemical processes represented in Equations (1)–(5) are as follows:

$$\begin{gathered} Photosynthesis (DIN,DIP,P)= \\ {V}_{\max }\exp (C_T^{P}T)\times \min \left(\frac{DIN}{DIN+K_{DIN}},\frac{DIP}{DIP+K_{DIP}}\right)\times \frac{I(P;z)}{I_{opt}}\exp \left(1-\frac{I(P;z)}{I_{opt}}\right)\times P, \end{gathered}$$

(6)

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: $\begin{array}{cc}{K}_{DIN}=V_{\max }/A+2\sqrt{V_{\max }\cdot DIN/A}& (\mu \mathrm{mol}\cdot {\mathrm{l}}^{-1})\end{array}$ and $\begin{array}{cc}{K}_{DIP}=V_{\max }/A+2\sqrt{V_{\max }\cdot R_{P:N}DIP/A}& (\mu \mathrm{mol}\cdot {\mathrm{l}}^{-1})\end{array}$, 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_T^P$ 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, ĸ(P). These relationships are given by

$$I(P;z)=I_s\exp \left(-{\displaystyle \underset{0}{\overset{z}{\int }}\kappa (P)dz}\right),$$

(7)

where ĸ(P) = c_dom + 0.054 P^0.6667 + 0.088 P, c_dom is the light dissipation coefficient of seawater. ĸ(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:

$$\mathit{Respiration} (P)=R\times \exp (C_T^{RP}T)\times P,$$

(8)

$$\mathit{Mortality}(P)=M_P\times \exp (C_T^{MP}T)\times P^2,$$

(9)

$$\mathit{Mortality}(Z)=M_Z\times \exp (C_T^{MZ}T)\times Z^2,$$

(10)

where $C_T^{RP}$, $C_T^{MP}$, and $C_T^{MZ}$ 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

$$\mathit{Grazing}(P,Z)=G_Z\times \max \{0, 1-\exp (\lambda (P_{}^{\ast }-P))\}\times \exp (C_T^{GZ}T)\times Z,$$

(11)

where $C_T^{GZ}$ 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

$$\mathit{Egestion} (Z)=(1-{\alpha }_Z)\times Grazing(Z),$$

(12)

$$\mathit{Excretion} (Z)=({\alpha }_Z-{\beta }_Z)\times Grazing(Z),$$

(13)

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:

$$\mathit{Decomposition} (D)=V_D\times \exp (C_T^{{\lambda }_D}T)\times D,$$

(14)

where V_D is the decomposition rate at 0 °C, and $C_T^{{\lambda }_D}$ is the temperature coefficient of decomposition.

Climatology damping conditions were adapted for the DIN and DIP source–sink terms, and were defined as follows:

$$\mathit{Clim}\_Damp(DIN)=\gamma (DIN-DI{N}_{\mathit{clim}}),$$

(15)

$$\mathit{Clim}\_Damp(DIP)=\gamma (DIP-DI{P}_{\mathit{clim}}),$$

(16)

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.

2.2. 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₃) in molar relationships as follows:

$$\frac{d\left(DIC\right)}{dt}=\frac{d\left(DIP\right)}{dt}\cdot R_{C:P}-\frac{d\left(CaC{O}_3\right)}{dt}-\mathit{Clim}\_Damp (DIC),$$

(17)

$$\frac{d\left(ALK\right)}{dt}=-\frac{d\left(DIN\right)}{dt}\cdot R_{ALK:N}-2\cdot \frac{d\left(CaC{O}_3\right)}{dt}-\mathit{Clim}\_Damp (ALK),$$

(18)

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:

$$\mathit{Clim}\_Damp(DIC)=\gamma (DIC-DI{C}_{\mathit{clim}}),$$

(19)

$$\mathit{Clim}\_Damp(ALK)=\gamma (ALK-AL{K}_{\mathit{clim}}),$$

(20)

where DIC_clim and ALK_clim are monthly climatology data for DIC and ALK, respectively (see Section 3.1).

Calcium carbonate (CaCO₃) immediately sinks after its production. We did not treat the CaCO₃ concentration as a predictable variable in this model. The equations for approximating CaCO₃ were as follows:

$$\frac{d\left(CaC{O}_3\right)}{dt}=\mathrm{Pr} \left(CaC{O}_3\right)-{\displaystyle \int \mathrm{Pr} \left(CaC{O}_3\right)dz\cdot \frac{d}{dz}\left(e^{-z/D_{CaC{O}_3}}\right)},$$

(21)

$$\mathrm{Pr} (CaC{O}_3)=(Egsetion(Z)+Mortality(P)+Mortality(Z))R_{CaC{O}_3/POC}{R}_{C:P}),$$

(22)

where D_CaCO3 is the CaCO₃ remineralization e-folding scale depth, and R_CaCO3/POC is the ratio of CaCO₃ to the non-photosynthetic production of particulate organic carbon (POC). Pr (CaCO₃) is parameterized as a fixed ratio (R_CaCO3/POC) of the production of non-diazotrophic detritus [18].

2.3. Inorganic Carbon Cycling in Water

The chemical reactions in water involve the carbon–water–borate system [47]. This system describes the relationships among [DIC], [ALK], [H⁺], and borate (BOR), whereby

$$[ALK]=\left(\frac{2k_1^{}{k}_2^{}{X}^2}{1+k_1^{}X+k_1^{}{k}_2^{}{X}^2}\right)[DIC]+\left(\frac{k_b^{}}{1+k_b^{}X}\right)[BOR]+k_w^{}X-X^{-1},$$

(23)

$$[BOR]=[B{(OH)}_4^{-}]+[H_{3}B{O}_3]\cong (0.00232/10.81)(S/1.80655),$$

(24)

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

$${[ALK]}_c=2[C{O}_3^{2-}]+[HC{O}_3^{-}]=\left(\frac{2k_1^{}{k}_2^{}{X}^2}{1+k_1^{}X+k_1^{}{k}_2^{}{X}^2}\right)[DIC].$$

(25)

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]:

$${\gamma }_c=\frac{[ALK]}{[DIC]}-\left(\frac{k_b^{}{X}^{}}{k_b^{}+1}\right)\frac{[BOR]}{[DIC]}-[k_w^{}X-X^{-1}]\frac{1}{[DIC]},$$

(26)

$$X=\frac{({\gamma }_c-1)k_1+{[{({\gamma }_c-1)}^2{k}_1^2+4{\gamma }_c(2-{\gamma }_c)k_1{k}_2]}^{1/2}}{2(2-{\gamma }_c)k_1{k}_2}.$$

(27)

γ_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:

$$[C{O}_2]=\frac{[DIC]}{1+k_{1}X+k_1{k}_2{X}^2},$$

(28)

$$P_{C{O}_2}=\frac{[DIC]/k_0}{1+k_{1}X+k_1{k}_2{X}^2},$$

(29)

$$C{O}_3^{2-}=\frac{k_1{X}^2[DIC]}{1+k_{1}X+k_1{k}_2{X}^2},$$

(30)

$${\Omega }_{\arg }=\frac{\frac{k_1{X}^2[DIC]293.86S}{k_{\arg }}}{1+k_{1}X+k_1{k}_2{X}^2}.$$

(31)

The apparent dissolution constants k₀, k₁, k₂, 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₂ gas was assumed to be proportional to the difference in the partial pressures of CO₂ in surface waters (PCO₂) and the atmosphere ($P_{C{O}_2}^{atm}$) [47,50], as described by

$$F_{DIC}={\rho }_w{V}_p{k}_{0C}(P_{C{O}_2}-P_{C{O}_2}^{atm}),$$

(32)

$$F_{ALK}=0,$$

(33)

where ρ_w is the seawater density, and k_0c is the solubility of CO₂ in seawater in mol kg⁻¹ atm⁻¹ [52,53,54], and is expressed as follows:

$$\begin{array}{l}{k}_{0C}=\exp (9345.17\left(\frac{1}{T}\right)-60.2409+23.3585\log \left(\frac{T}{100}\right)\\ +S\left(0.023517-0.023659\left(\frac{T}{100}\right)+0.0047036{\left(\frac{T}{100}\right)}^2\right)),\end{array}$$

(34)

where T is in units of Kelvin, and V_p is the piston velocity, which depends on wind speed U₁₀ and the Schmidt number (S_c), and is expressed in units of m s⁻¹ [53] as follows:

$$V_p=8.61\times {10}^{-7}{U}_{10}^2{(S_c/660)}^{-1/2}.$$

(35)

S_c is a function of the ambient temperature T (in °C), and, for CO_2, is expressed as follows:

$$S_c=2073.1-125.62 T+3.627 T^2-0.043219 T^3.$$

(36)

3. Data and Model Parameters

3.1. 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₂ 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_{C{O}_2}^{atm}$) are available at the JMA website. Using data from the Ayasato Observatory in Japan, we created an empirical formula for CO₂ (in ppm) with the addition of a second-order polynomial (CO₂^poly) and seasonal sinusoidal curve equations (CO₂^curv), as follows:

$$C{O}_2^{poly}=0.000633\times t_{yr}^2-0.650417\times t_{yr}-859.0625,$$

(37)

$$\begin{array}{ll}C{O}_2^{curv}=& -0.00064+6.09254\cos \left(\frac{2\pi }{24\times 365}{t}_{hr}-55.20975\frac{2\pi }{360}\right)\\ & +0.39431\cos \left(\frac{2\pi }{24\times 91.25}{t}_{hr}-120.3614\frac{2\pi }{360}\right),\end{array}$$

(38)

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).

3.2. Biogeochemical Parameters and Observational Data for Validation

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. Biogeochemical parameters obtained from optimized parameterization using chlorophyll-a, for subarctic and subtropical regions.

	Chlorophyll-a
	Subarctic Region	Subtropical Region
Vmax	0.692	0.692
R	0.02	0.04
MP	0.045	0.3
GZ	0.05	0.05
MZ	0.05	0.05
VPN	0.05	0.05
WD	6.79	6.78

). 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.

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:

$$I_{opt}=0.5I_{opt}^0(\tanh ((Lat-La{t}_{bnd})/La{t}_{slp})+1)+I_{opt}^1,$$

(39)

$$c_{dom}=0.5c_{dom}^0(\tanh ((Lat-La{t}_{bnd})/La{t}_{slp\_cdom})+1)+c_{dom}^1,$$

(40)

where I_opt and c_dom change latitudinally from 20 to 120 W/m⁻² and from 0.015 to 0.045 m⁻¹, respectively; I⁰_opt, I¹_opt, c⁰_dom, and c¹_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 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⁻¹, after adopting the latitudinal differences for c_dom and I_opt (Appendix A).

3.3. 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.

4. Results

4.1. 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 Figure 4. The satellite data show contrasting larger and smaller values in the subarctic and subtropical regions, respectively, throughout the year (Figure 4a–d and Figure 5a). The chlorophyll-a gradually increased to 2.0–3.0 milli mol/m³ at 50° N from March to May, and to 1.0–2.0 milli mol/m³ at 40° N, and remained high above 50° N during summer (June to July), with double peaks in spring and autumn around this latitude. There were also double peaks at 40° N in spring (March to May) and autumn (September to October) (Figure 5a). The model outputs basically reflected the observed chlorophyll-a distribution. The chlorophyll-a values were relatively large for areas above the region of 40° N, and were lower for areas at lower latitudes (Figure 4e–h). The model gave reasonable simulations of the chlorophyll-a variability, and the ranges were similar to the observed values (Figure 4e–h), but gave poor simulations of the seasonal patterns, particularly in the subarctic region (Figure 5a,b). The model reproduced only one notable peak in summer from May to August at 50° N (Figure 5b), and this peak lagged behind the observed data by 1–2 months (Figure 5a,b). Small and large double peaks were reproduced at 30° N and at 40° N, respectively. Consistent with the observed data (Figure 5a), the model reproduced small peaks for April and October at 20° N (Figure 5b).

The model gave adequate simulations of DIN and reproduced the basic features and seasonal variability, when the monthly climatology data were dampened (Figure 5c,d and Figure 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 (Figure 5c and Figure 6a–d) was reproduced by the model (Figure 5d and Figure 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–c). The model simulations of the patterns of ALK in the subarctic region above 45° N in winter, spring and autumn (Figure 8d–f), and in the subtropical region in summer (Figure 8f), were similar. However, unlike the observed data, there was a marked decrease (50 μmol/kg) for summer in the subarctic region and the Okhotsk Sea (Figure 8f). The low values in summer were followed by higher values in autumn (Figure 8g). The time series for ALK (Figure 5g) also showed a decrease in summer around 40° N and 50° N. The summer decrease arose from the modeled variations in CaCO₃ (not shown), and indicates that the summer decrease in ALK may have been overestimated.

4.2. Ocean Acidification Indices pH and Ω_arg

Figure 9 shows the surface pH_insitu (pH₂₅) 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,b). In the Kuroshio Extension region (140°–180° E, 30° N), the pH values were low during July (summer). The values were highest (8.10–8.15) in winter (January) and spring (April) (Figure 9a) between the subtropical and subarctic regions in the area corresponding to the Kuroshio extension. The values varied with latitude along the 165° E line (Figure 10a). The highest and lowest values occurred in February or March and June or August, respectively, in lower latitude regions around 20° N and 30° N, but the opposite occurred at around 40° N and 50° N (Figure 10a).

The pH₂₅ 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₂₅ values increased in summer (Figure 9e–h, Figure 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₂₅ (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–d,i–l). Yara et al. [65] studied the Ω_arg distributions in a marginal area around Japan (24°–48° N, 118°–157° E). Their isothermal contours (Ω_arg = 2.3 and 3.0), which reproduced the conditions in the 2010s (their Figure 7b), were similar to our model outputs (Figure 9i–l).

4.3. 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₂₅ and ALK (

Table 3. Correlation coefficients for chlorophyll-a (Chl-a), DIN, DIC, ALK, pH₂₅ and Ω_arg between JMA observed data and the model outputs for each region in 2015 (Figure 1). The observed data and model outputs that matched in time and place were compared. Values in brackets indicate p-values the significance level of 0.05 and the significant level of p < 0.05. Only for Chl-a, log₁₀ (Chl-a) values were evaluated.

Parameter	Subtropical Region	Subarctic Region	Kuroshio Extension	Japan Sea
Chl-a	0.57 (0.07)	0.88 (0.07)	0.86 (0.12)	0.87 (0.15)
DIN	0.72 (0.04)	0.88 (0.12)	0.87 (0.10)	0.95 (0.11)
DIC	0.74 (0.09)	0.87 (0.20)	0.78 (0.17)	0.92 (0.26)
ALK	0.48 (0.09)	0.69 (0.20)	0.48 (0.17)	0.35 (0.26)
pH25	0.78 (0.09)	0.86 (0.20)	0.81 (0.17)	0.90 (0.26)
Ωarg	0.99 (0.09)	0.98 (0.20)	0.98 (0.17)	0.99 (0.26)

). 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 correlation values were 0.57, 0.88, 0.86 and 0.87 for chlorophyll-a; 0.72, 0.88, 0.87 and 0.95 for DIN; 0.74, 0.87, 0.78 and 0.92 for DIC; 0.48, 0.69, 0.48 and 0.35 for ALK; 0.74, 0.86, 0.81 and 0.90 for pH₂₅; and 0.99, 0.98, 0.98 and 0.99 for Ω_arg, in the subtropical region, the subarctic region, the Kuroshio Extension and the Japan Sea, respectively (Table 2). Scatter plots (Figure 11, Figure 12 and Figure 13) showed that the parameters were highly correlated with each other (not shown for DIP, DIC and ALK).

5. Discussion

5.1. 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 (Figure 4e–h, Figure 6e–h, Figure 7e–h and Figure 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., Figure 6e and Figure 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).

5.2. 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₂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₂O←H⁺+OH⁻ (H₂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₂ equilibrium in sea water, defined by: CO₂+H₂O↔H⁺+HCO₃⁻. 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: CO₂+H₂O→H⁺+HCO₃⁻ (CO₂+H₂O←H⁺+HCO₃⁻). The decrease 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: CO₂+H₂O←H⁺+HCO₃⁻.

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₂₅ 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.

5.3. 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²⁺][CO₃²⁻])/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₂+H₂O+CO₃⁻↔2HCO₃⁻, 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₃⁻) 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) was mainly the result of the decrease in DIC in summer (Figure 9i–l and Figure 10c), because the model probably overestimated a decrease in ALK in summer and the decrease in ALK and the increase in temperature in summer together caused a decrease in Ω_arg.

6. 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 (Figure 4 and Figure 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₃ and also ALK reduction (Figure 5g and Figure 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 Figure 16 and Figure 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.