sustainability
ArticleDevelopment of a Biogeochemical and Carbon
Model Related to Ocean Acidification Indiceswith an Operational Ocean Model Product in the
North Western PacificMiho Ishizu 1, *, Yasumasa Miyazawa 1 , Tomohiko Tsunoda 2 and Xinyu Guo 1,3
1

2
3

*

Application Laboratory, Japan Agency for Marine-Earth Science and Technology, 3173-25, Showa-machi,
Kanazawa-ku, Yokohama-city, Kanagawa 236-0001, Japan; miyazawa@jamstec.go.jp (Y.M.);
guo.xinyu.mz@ehime-u.ac.jp (X.G.)
Ocean Policy Research Institute, Sasakawa Peace Foundation, 1-15-16 Toranomon, Minato-ku,
Tokyo 105-0001, Japan; t-tsunoda@spf.or.jp
Center for Marine Environmental Studies, Ehime University, 2-5, Bunkyo-ku, Matsuyama,
Ehime 790-8577, Japan
Correspondence: mishizu@jamstec.go.jp; Tel.: +81-45-778-5875

Received: 28 February 2019; Accepted: 6 May 2019; Published: 10 May 2019


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 resultsof 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 temperatureand dissolved inorganic carbon. The seasonal increase in Ωarg from winter to summer was governed
mainly by dissolved inorganic carbon levels.
Keywords: biogeochemical model; carbon; NPZDC; North Western Pacific; ocean acidification; JCOPE;
pH; aragonite saturation

1. Introduction
Ocean acidification poses a serious risk to marine organisms and ecosystems, including finfish andcoral reefs in subtropical regions, and species or groups of organisms in polar regions [1–5]. The effectsof ocean acidity have previously been reported and the effects of acidification are predicted to increase,with great risks to marine organisms [2,4–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, andsubstantial economic losses are likely [8]. Therefore, these issues can no longer be ignored, and urgent
action is necessary.
Sustainability 2019, 11, 2677; doi:10.3390/su11092677

www.mdpi.com/journal/sustainability

Sustainability 2019, 11, 2677

2 of 28

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 tochange even if new automatic monitoring systems, such as biological profiling floats [10], are establishedand operational monitoring, using volunteer observing ship (VOS) pCO2 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–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–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 thecarbon cycle. Zhao and Guo [28] established an inorganic carbon system module for the Yellow andEast 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 andbiogeochemical 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–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 anoperational 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 andsalinity data obtained from floating profilers and/or ship measurements. The JCOPE2M covers themarginal 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, mesoscaleeddies, 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, andsea 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 forthe 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 theclimatological information. Sensitivity analysis was performed using fixed values for temperature,

Sustainability2019,2019,11,
11,2677
2677
Sustainability

33 of 28
30

climatological information. Sensitivity analysis was performed using fixed values for temperature,salinity,
andandALK,
to assesswhich which
variables
were importantin determining
the ocean acidification
salinity,DIC,
DIC,
ALK,
to assessvariables
were importantin determining
the ocean
indices,pH
andaragonitesaturation(Ω).
acidification indices, pH and aragonite arg
saturation (Ωarg).
This
Thisisisthethefirst
firsttime
timethethe development
developmentofof an
an operational
operational carbon
carbon cycle
cycle model
model with
with eddy-resolving
eddy-resolving
horizontal
resolution
focused
on
the
North
Western
Pacifichas
beenreported.
Somefundamental
data
horizontal resolution focused on the North Western Pacific has been reported. Somefundamental
required
for
our
study
(e.g.,
climatology
data)
had
to
be
acquiredbecause
only
a
few
of
the
necessarydata required for our study (e.g., climatology data) had to be acquired because only a few of thedata
for researchinresearch
this challenging
scientific field
were available.
necessarydata forin this challenging
scientificfield were available.
Details
of
the
model
formulation
areprovided
Section
2. Themodel
parameters
Details of the model formulation are provided inin
Section
2. Thedatadata
and andmodel
parameters
used
used
in study
the study
described
in Section
Sectiondescribeshow
howthethe model
model was
was validated
in theare aredescribed
in Section
3. 3.
Section4 4describes
validated withwithobserved
data.
In
Section5,
theinfluenceofsomeenvironmental
variableson
ocean
acidification
observed data. In Section 5, the influence of some environmental variables on ocean acidification
indices) are discussed. The main findings of the study are overviewed, and future issues
indices(pH(pHandandΩΩargarg) are discussed. The main findings of the study are overviewed, and future issues
arearediscussed,
discussed,in
inSection
Section6.6.

Figure 1.
1. ModelMeteorological
Agency(JMA)
observationsitessites
in 2015.
The
FigureModeltarget
targetarea
areaandandJapanJapan
Meteorological
Agency(JMA)
observationin 2015.
coloredred points
indicate
the JMAobservation
points.
The redpointsthe shaded
blue, orange,Thecoloredred points
indicate
the JMAobservation
points.
The
red inpointsin the shaded
blue,pink andgreen
regions
were wereusedusedin thecalculation
of of
optimized
correlation
orange,pink
and green
regionsin thecalculation
optimizedoperations
operationsand
and ofof correlation
coefficients
betweenobservation
and
model
outputsforthesubarctic
region,theKuroshioextension,coefficients between observation and model outputs for the subarctic region, the Kuroshio extension,
thesubtropicalsubtropicalregion
regionandandthetheJapanJapanSea,
Sea,respectively
respectively(Tables(Tables22andand3).
3).
the2.
2. ModelModel Formation
FormationThe
constructed
on antracer model
physical
TheJCOPE_EC
JCOPE_ECwas
wasconstructed
onoffline
an offline
tracer driven
model bydrivenby processes
physical simulated
processesbythe modelof theJCOPE2M,
an operational
eddy-resolving
ocean generalcirculation
simulatedby outputs
the model
outputsof the and
JCOPE2M,
and an operational
eddy-resolving
ocean
general
model
based
on
the
Princeton
Ocean
Model,
withageneralized
sigma
coordinate
[38].
The
JCOPE_ECiscirculation model based on the Princeton Ocean Model, with a generalized sigma coordinate [38].
◦ –62◦ N,
aThethree-dimensional
high-resolution
regional
model
thatcoversthewestern
North
Pacific(10.5
JCOPE_EC is a three-dimensional high-resolution regional model that covers the western North◦ –180◦ E) at a horizontal resolution of 1/12◦ (9.1–4.4 km) and has 46 vertical active levels.
108
Pacific(10.5°–62°N, 108°–180°E) at a horizontal resolution of 1/12° (9.1–4.4 km) and has 46 verticalactive levels.
2.1. Equations Governing the Marine Ecosystem Model
2.1. Equations
the Marine
Ecosystem
Model model (NPZD) was used. It comprised five
A simpleGoverning
nitrogen-basedplanktonecosystemcomponents
includingdissolved
inorganic
nitrogen model(DIN, NO
inorganic
phosphate3 ), dissolved
A simple
nitrogen-basedplanktonecosystem(NPZD)
was used.
It comprised
five
(DIP,
PO
),
phytoplankton(P),zooplankton(Z)anddetritus(D)[26,27].
The
carbon
cycle
model
4
components
including dissolved inorganic nitrogen (DIN, NO3), dissolved inorganic phosphate(DIP,
comprised
three components:
DIC, ALK(CaCO
was coupledwith the3 ), andcyclePO4), phytoplankton(P), zooplankton(Z)andandcalcium
detrituscarbonate(D) [26,27].
The carbon
model comprised
NPZD
model.
Detrituscorresponded
totheconcentrationsofparticulateorganic
nitrogen
three components: DIC, ALK and calcium carbonate (CaCO3), and was coupled with the NPZD (PON).
model.
The
pH (pHandtemperature,salinity,DIC
arg values were
insitu , pH25 )toDetrituscorresponded
theΩconcentrations
ofcalculatedparticulatefrom
organic
nitrogen (PON).
The
pHand
(pHALK
insitu,usingaFortranprogram
based
on
the
CO2sys
MATLABsoftware[39,40].
The
developmentofthepH25) and Ωarg values were calculated from temperature, salinity, DIC and ALK using a Fortran
nitrogen-basedplanktonecosystem
model
followed
theoutlines
described
by
Guo
and
Yanagi[27],
program based on the CO2sys MATLAB software [39,40]. The development of the nitrogen-basedplankton ecosystem model followed the outlines described by Guo and Yanagi [27], Onitsuka andSustainability 2019, 11, 2677

4 of 28

Onitsuka and Yanagi [26] and Sasai et al. [25] (Figure 2a). The carbon cycle model is referred fromOrr 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= Phtosynthesis (DIN, DIP, P) − Respiration (P) − Mortality (P) − Grazing (P, Z),
dt

(1)

dZ
= Grazing (P, Z) − Mortality (Z) − Egestion (Z) − Excertion (Z),
dt

(2)

dD
∂
= Egestion (Z) + Mortality (P) + Mortality (Z) − Decomposition (D) − (WD D),
dt
∂z

(3)

d(DIN )=
dt

−Phtosynthesis (DIN, DIP, P) + Respiration (P) + Excertion (Z)+Decomposition (D) − Clim_Damp(DIN )d(DIP)
1 d(DIN )=− Clim_Damp(DIP),
dt
RP:N
dt

(4)
(5)

where W D is the sinking velocity of detritus, RP:N is the Redfield ratio of P (phosphorus):N (nitrogen)(Table 1). The biogeochemical processes represented in Equations (1)–(5) are as follows:
Photosynthesis (DIN, DIP, P) = I(P;z)I (P;z)× Iopt exp 1 − Iopt × P,

DIP
Vmax exp(CPT T ) × min DINDIN+KDIN , DIP+KDIP(6)

where T is the seawater temperature, V max is the maximum growth rate from photosynthesis,
and KDIN and KDIP are half saturation constants for DIN and DIP, respectively. We adopted
optimal uptake kinetics for KDIN relating to V max , according to the theory of physiological√acclimation of phytoplankton [25,43], defined as: KDIN = Vmax /A + 2 Vmax · DIN/A (µmol · l−1 )√and KDIP = Vmax /A + 2 Vmax · RP: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,
RP:N is the stoichiometry of N to P [44,45], CPT is a temperature-dependent coefficient for photosynthesis,
and Iopt 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 (Is ) and the light extinction coefficient, k(P). These relationships are given by
 z

 Z



I (P; z) = Is exp− κ(P)dz,



(7)

0

where k(P) = cdom + 0.054 P0.6667 + 0.088 P, cdom is the light dissipation coefficient of seawater.
k(P) depends on cdom 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:
Respiration (P) = R × exp(CRP
T T ) × P,

(8)2
Mortality(P) = MP × exp(CMP
T T) × P ,(9)2
Mortality(Z) = MZ × exp(CMZ
T T) × Z ,(10)

where CRP
, CMP
, and CMZ
are temperature-dependent coefficients for the respiration of phytoplankton,
TT
T
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 theirconcentrations (biomass).

Sustainability 2019, 11, 2677

5 of 28

The grazing of phytoplankton by zooplankton was defined as a function of temperature and theconcentrations of phytoplankton and zooplankton, expressed asGrazing(P, Z) = GZ × max 0, 1 − exp(λ(P∗ − P)) × exp(CGZ
T T ) × Z,

(11)

where CGZ
is the temperature dependent-coefficient of grazing by zooplankton, Gz is the maximumTrate 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 andpredating, as defined by
Egestion (Z) = (1 − αZ ) × Grazing(Z),(12)
Excretion (Z) = (αZ − βZ ) × 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:λ

Decomposition (D) = VD × exp (CTD T ) × D,

(14)

λ

where VD is the decomposition rate at 0 ◦ C, and CTD is the temperature coefficient of decomposition.
Climatology damping conditions were adapted for the DIN and DIP source–sink terms, and were
defined as follows:
Clim_Damp(DIN ) = γ(DIN − DINclim ),(15)
Clim_Damp(DIP) = γ(DIP − DIPclim ),(16)

where γ is the inverse of the time scale (= 30 days), and DINclim and DIPclim are the monthly climatology
data forDIN and
DIP,
respectively. Some details of climatology damping are described in Section
Sustainability
2019,11, 2677
6 of 30 3.1.

Figure 2. Schematic view of the basic components of the current biogeochemical model: (a) model of
Figure 2. Schematic view of the basic components of the current biogeochemical model: (a) model ofbiogeochemical and carbon cycling; and (b) the ocean carbon cycle model.
biogeochemical and carbon cycling; and (b) the ocean carbon cycle model.

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

Value

Units

Reference Values
Onitsuka andSustainability 2019, 11, 2677

6 of 28

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

Value

Units

Reference
Values Onitsuka
and Yanagi [26]Growth rate for phytoplankton
Affinity coefficient of basic cellular physiology
Half saturation constant for DIN
Half saturation constant for DIP
Phytoplankton mortality rate at 0 ◦ C
Threshold of phytoplankton mortality
Phytoplankton respiration rate at 0 ◦ C
Temperature coefficient for photosynthesis
Temperature coefficient for phytoplankton
respiration
Temperature coefficient for phytoplanktonmortality
Optimum light intensity for phytoplankton

0.04926.75
1.50
0.09400.04000.05870.03170.0392

day−1
mmolN−1 m−1 day−1
mmolN m−1
mmolN m−1(mmolNm−1 ) m−3 day−1(mmolNm−1 ) m−3
day−1
◦ C−1

0.6
−
1.5
0.09
0.07
−
0.3
0.06930.0519◦ C−1

0.05190.0693◦ C−1

0.0693* 20–120
100.0
70.0
45.0
4.0
* 0.015–0.045
0.02750.0250
1.5

W m−2

70
−
−
−
−
0.04
−
−
−

0.05
1.4
0.05
0.3
0.7

day−1
(mmolN m−3 )−1
◦ C−1

0.3
1.4
0.07
0.3
0.7

0.0430(mmolNm−1 ) m−3

0.043

0.0390◦ C−1

0.06930.0693◦ C−1

0.0693Singing velocity of detritus
Decomposition rate at 0 ◦ C (DET→DIN)
Temperature coefficient for decomposition

6.7
0.050
0.0693m day−1
day−1
◦ C−1

10
0.05
0.0693

Stoichiometry of nitrogen to phosphorus [46]Molar elemental ratios
CaCO3 over nonphotosynthetical POC
production ratio
Alkalinity over non-photosynthetic N
production ratio
CaCO3 remineralization e-folding depth

16.0
112.0

−
−

0.035

−

Ecosystem Model
For phytoplankton
Vmax
A
KDIN
KDIP
MP
Pmin
R
CP T
CRP T
CMP T
Iopt
Iopt 0Iopt 1
latbndlatslpcdom
cdom 0cdom 1
Latslp_dom
For zooplankton
GZ
Λ
MZ
βz
αz
P*
CGZ T
CMP T
For diatoms
WD
VPN
CλD T
Carbon Cycle Model
RP:N
RC:P
RCaCO3/POC
RALK:N
DCaCO3

Boundary for latitudinal differences
Slope for latitudinal differences
Light dissipation coefficient of seawater

Slope for latitudinal differences for cdom
Maximum grazing rate of zooplankton at 0 ◦ C
Ivlev constant
Zooplankton mortality rate at 0 ◦ C
Growth efficiency of zooplankton
Assimilation efficiency of zooplankton
Zooplankton threshold value for grazing on
phytoplankton
Temperature coefficient for zooplankton
grazing
Temperature coefficient for zooplanktonmortality

m−1
m−1
m−1

−

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

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 releaseoccur in fixed elemental ratios for C, P and/or N. The production of DIC and ALK was controlled bychanges in the inorganic nutrients (DIN, DIP) and calcium carbonate (CaCO3 ) in molar relationships
as follows:d(DIC)d(DIP)d(CaCO3 )=· RC:P −− Clim_Damp (DIC),
(17)
dt
dt
dt
d(ALK)d(DIN )d(CaCO3 )=−· RALK:N − 2 ·− Clim_Damp (ALK),dt
dt
dt
where RC:P and RALK:N are the molar elemental ratios of C to P, and of ALK to N, respectively.

(18)

Sustainability 2019, 11, 2677

7 of 28

Climatology damping conditions were also adapted for the DIC and ALK source–sink terms,
as follows:
Clim_Damp(DIC) = γ(DIC − DICclim ),(19)
Clim_Damp(ALK) = γ(ALK − ALKclim ),(20)

where DICclim and ALKclim are monthly climatology data for DIC and ALK, respectively (see Section 3.1).
Calcium carbonate (CaCO3 ) immediately sinks after its production. We did not treat the CaCO3
concentration as a predictable variable in this model. The equations for approximating CaCO3 wereas follows:Zd(CaCO3 )d  −z/DCaCO 
3 ,
e
(21)
= Pr (CaCO3 ) − Pr (CaCO3 )dz ·dtdzPr (CaCO3 ) = (Egsetion(Z) + Mortality(P) + Mortality(Z))RCaCO3 /POC RC:P ),(22)

where DCaCO3 is the CaCO3 remineralization e-folding scale depth, and RCaCO3/POC is the ratio of CaCO3
to the non-photosynthetic production of particulate organic carbon (POC). Pr (CaCO3 ) is parameterized
as a fixed ratio (RCaCO3/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 systemdescribes the relationships among [DIC], [ALK], [H+ ], and borate (BOR), whereby!
kb
2k1 k2 X2[DIC] +[ALK] = [BOR] + kw X − X−1 ,1 + kb X
1 + k1 X + k1 k2 X2 

(23)[BOR] = [B(OH )−4 ] + [H3 BO3 ]  (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[CO2−3 ] + [HCO3 ] = 

2k1 k2 X2[DIC].
1 + k X + k k X2 
1

(25)

1 2

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 =

X=

!
kb X [BOR]
[ALK]
1
−
− [kw X − X−1 ],kb + 1 [DIC]
[DIC]
[DIC]

(γc − 1)k1 + [(γc − 1)2 k12 + 4γc (2 − γc )k1 k2 ]2(2 − γc )k1 k2

(26)

1/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:[CO2 ] =[DIC],
1 + k1 X + k1 k2 X2(28)

PCO2 =[DIC]/k0,
1 + k1 X + k1 k2 X2(29)

CO32− =

k1 X2 [DIC],
1 + k1 X + k1 k2 X2(30)

Sustainability 2019, 11, 2677

8 of 28Ωarg =

k1 X2 [DIC]293.86S
karg

1 + k1 X + k1 k2 X 2

.

(31)

The apparent dissolution constants k0 , k1 , k2 , kb , kw , and karg 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 CO2 gas was
assumed to be proportional to the difference in the partial pressures of CO2 in surface waters (PCO2 )and the atmosphere (Patm) [47,50], as described by
CO
2

FDIC = ρw Vp k0C (PCO2 − Patm
CO ),

(32)

FALK = 0,(33)2where ρw is the seawater density, and k0c is the solubility of CO2 in seawater in mol kg−1 atm−1 [52–54],
and is expressed as follows:

 
 
T
k0C = exp 9345.17 T1 − 60.2409 + 23.3585 log 100  2 
T
T+S 0.023517 − 0.023659 100+ 0.0047036 100,(34)

where T is in units of Kelvin, and Vp is the piston velocity, which depends on wind speed U10 and theSchmidt number (Sc ), and is expressed in units of m s−1 [53] as follows:2
Vp = 8.61 × 10−7 U10(Sc /660)−1/2 .

(35)

Sc is a function of the ambient temperature T (in ◦ C), and, for CO2, is expressed as follows:
Sc = 2073.1 − 125.62 T + 3.627 T2 − 0.043219 T3 .

(36)

3. Data and Model Parameters3.1. Climatology and Forcing Data
We created monthly DIN and DIP climatological data on the model grid by interpolating theoriginal 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 thesurface values for DIN and DIP in the WOA13 with the data from Yasunaka’s study, in an area definedby 10.5◦ –60◦ N and 120◦ –180◦ E.
A climatological DIC database was created by combining an annual mean climatology productestimated 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] withtemperature and salinity data from the WOA13 climatology and the annual mean climatology ofKey et al. [57], and then replaced the surface data for depths above 200 m with the data estimatedusing 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

Sustainability 2019, 11, 2677

9 of 28

observed data, on the assumption that the JMA data had the characteristics of ALK profiles in theJapan 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 CO2 gas exchange in thesurface/subsurface layer. Short wave radiation (Is ) 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 (Patm) are available at the JMA
CO2
website. Using data from the Ayasato Observatory in Japan, we created an empirical formula for
CO2 (in ppm) with the addition of a second-order polynomial (CO2 poly ) and seasonal sinusoidal curve
equations (CO2 curv ), as follows:
poly

CO2

COcurv=
2= 0.000633 × t2yr − 0.650417 × t yr − 859.0625,

(37)



2π
2π
thr − 55.20975 360−0.00064 + 6.09254 cos 24×365


2π
2π
+0.39431 cos 24×91.25
thr − 120.3614 360,(38)

where tyr and thr 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 fortracers. 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). The observed
data were downloaded from the JMA website (Figure 1). The phytoplankton concentrations werealso validated with chlorophyll-a data downloaded from an File Transfer Protocol (FTP) site for theModerate Resolution Imaging Spectroradiometer (MODIS)-Aqua Ocean Color Data.
In previous studies, where biogeochemical models were applied to the North Pacific and themarginal 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–23,25], the subtropical area [21,25],
and/or the marginal seas [27,60] (Figure 1). In contrast, our model included both subarctic andsubtropical areas and adjoining marginal seas. As biogeochemical situations differ depending on
oceanic conditions [61–63], in our model we applied latitudinal differences for Iopt and cdom , whichwere defined as follows:
01
Iopt = 0.5Iopt(tanh((Lat − Latbnd )/Latslp ) + 1) + Iopt,

(39)

cdom = 0.5c0dom (tanh((Lat − Latbnd )/Latslp_cdom ) + 1) + c1dom ,(40)

where Iopt and cdom change latitudinally from 20 to 120 W/m−2 and from 0.015 to 0.045 m−1 , respectively;
I0 opt , I1 opt , c0 dom , and c1 dom are parameters for adjusting from subtropical to subarctic regions (Table 1
and Figure 3); and Latbnd and Latslp (Latslp_cdom ) are the coefficients for latitudinal boundaries and thelatitudinal slopes for these parameters, respectively (Table 1 and Figure 3). The effects of latitudinaldifferences on Iopt and cdom 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 inthese optimization operations, and we adopted the optimized operations for the subarctic andsubtropical regions separately. All the estimated biogeochemical parameters were found to be similar
to those forthe11,subarctic
and subtropical regions (Table 2), but we manually adjusted Vmax to10
0.0492
Sustainability
2019,
2677
of 28
−1
day , after adopting the latitudinal differences for cdom and Iopt (Appendix A).

(a)

(b)

Figure3.3.Latitudinal
Latitudinaldifferences
differencesinin
biogeochemicalparameters,dissipation coefficient
coefficient
Figurethethebiogeochemicalparameters,(a) the light dissipationcdom)) andand(b)
thetheoptimum
lightintensity(Iopt(I).opt ).
(c(domoptimum
lightintensity
Table2.2.Biogeochemical
Biogeochemicalparametersparametersobtained
obtainedfrom
fromoptimized
optimizedparameterization
parameterizationusing
usingchlorophyll-a,chlorophyll-a,
Table
forsubarctic
subarcticand
andsubtropicalsubtropicalregions.
regions.
for

Vmax
R
MP
GZ
MZ
VPN
WD

Vmax
R
MP
GZ
MZ
VPN
WD

Chlorophyll-a
Chlorophyll-a
Subarctic region Subtropical region
Subarctic Region
Subtropical Region
0.692
0.692
0.692
0.692
0.02
0.04
0.02
0.04
0.045
0.3
0.045
0.3
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
0.05
6.79
6.78
6.79
6.78

3.3. Initial and Boundary Conditions
Using the Green’s function method [36], we optimized the biogeochemical parameters Vmax , R, Mp ,used a 1-year period, 2015, for all the experiments, which were driven by the forcing data
Gz , MWez , VPN , and WD (Equations (6)–(16) and Table 1). Chlorophyll-a data for 2015 were used in these
for that year.operations,
In a preliminary
phytoplankton
wererepresented
optimization
and we experiment,
adopted the theoptimized
operations
for thesubarctic by
andchlorophyll-a
subtropicalconcentrations
from
the
World
Ocean
Atlas
2001.
There
was
undesirable
noise
in
the
model
outputs
regions separately. All the estimated biogeochemical parameters were found to be similar
to thoseforfor
phytoplankton
inthepreliminary
experiment,
so,
for
the
initialconditions,
i.e.,
the
beginningof−1
the subarctic and subtropical regions (Table 2), but we manually adjusted Vmax to 0.0492 day , after
the first monthin the numericalexperiments,
values for
adopting
the latitudinaldifferences
for cdom andconstantIopt (Appendix
A).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
3.3.
Initial and Boundary
Conditions
concentration. The initial detritus concentration was set at 0.0 mmol/kg.
We used a 1-year period, 2015, for all the experiments, which were driven by the forcing data
As a first step towards developing reliable operational systems for biogeochemical parameters
for that year. In a preliminary experiment, the phytoplankton were represented by chlorophyll-a
related to ocean acidification indices, we adopted a diagnostic approach to obtain realistic
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 thefirst month in the numerical experiments, constant values for phytoplankton of 0.05 and 0.0 µg/L wereset 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Sustainability 2019, 11, 2677

11 of 28

using the annual climatology at the beginning. The climatology data for DIN, DIP, DIC and ALK wereapplied 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 lateralboundaries 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 theforcing created from the ocean (JCOPE2M) and atmospheric (NCEP/NCAR) reanalysis data.
4. Results
4.1. Comparison
the Observed Data with the Simulated Outputs
Sustainability
2019, 11,of
2677

13 of 30

The monthly mean distributions of satellite and modeled chlorophyll-a values at the surface

showed a relatively larger decrease in DIC in summer (from July to September; Figure 5f) compared
are shown in Figure 4. The satellite data show contrasting larger and smaller values in the subarctic
with the observed data for 40° N (Figure 5e).
and subtropical regions, respectively, throughout the year (Figures 4a–d and 5a). The chlorophyll-a
Due toincreaseda lack ofobserved
could
not Marchcomparetheandseasonal
variations
in3 ALK
3 at 50◦ N from
gradually
to 2.0–3.0
millidata,
mol/mweto May,to 1.0–2.0
milli mol/m
at[13,26,52,53].
Instead
of
using
the
ALK
climatology,
weplottedtheJMAobserved
data(Figure8a–c).
◦◦
40 N, and remained high above 50 N during summer (June to July), with double peaks in spring
These
showedthat
ALK concentrations
relatively
μmol/kg)in theand data
autumnaroundthisthelatitude.
There were alsowere
doublepeaks at low
40◦ N(2200–2250
in spring (March
to May)subarctic
region
above
40°N
and
in
the
Okhotsk
Sea,
and
were
higher(2250–2300μmol/kg)
below
and autumn (September to October) (Figure 5a). The model outputs basically reflected the observed
40°N
from spring
to autumn
8a–c). values
The model
simulations
patternsof ALK
chlorophyll-a
distribution.
The(Figure
chlorophyll-awere relatively
largeofforthe
areas
above theregioninofthe◦
subarctic
region45°areas
N inspring(Figure
and autumn(Figureand in thesubtropical
40 N, andwereabove
lower forat winter,lower latitudes4e–h). The
model8d–f),gave reasonable
simulations
regioninchlorophyll-asummer (Figure8f), weresimilar.
However,
unlike theobserved
data,
there(Figure
was a 4e–h),markedof thevariability,
andthe ranges
were similar
to theobserved
valuesbut gavesimulationsof the seasonalpatterns, region
particularlyin thesubarctic(Figure5a,b).
decrease(50poorμmol/kg)
for summer
in the subarctic
and theOkhotsk
Searegion(Figure8f). Thelow
The model
reproduced
one notablepeak values
in summerfrom May
to August50◦ time
N (Figure5b),valuesin summer
were only
followedby higher
in autumn(Figure8g).atTheseriesforandALK
this peakbehindatheobserved
by 1–2months(FigureSmalllarge double
peaks(Figure5g) lagged
also showeddecreasein datasummeraround40°
N and5a,b).
50° N.
Theandsummerdecrease
arose
◦
◦
were
at 30 N and
40 3N,
respectively.
Consistentwiththeobserved(Figurein5a),fromthereproduced
modeled variations
in at
CaCO
(not
shown), andindicatesthatthesummerdata
decrease
ALK◦ N (Figure 5b).
the
model
reproduced
small
peaks
for
April
and
Octoberat
20
may have been overestimated.

Figure4. 4.
Surface
distributions
concentrationsfrom
from
MODIS_Aqua
FigureSurface
distributionsofofthe
themonthly
monthlymean
mean chlorophyll-achlorophyll-a concentrations
MODIS_Aqua
data,
and
from
model
outputsforJanuary(a,e),
April(b,f),
July
(c,g)
and
October(d,h).
data, and from model outputs for January
April (b,f), July (c,g) and October (d,h).

Sustainability 2019, 11, 2677

12 of 28

The model gave adequate simulations of DIN and reproduced the basic features and seasonalvariability, 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
nability 2019, 11,2677
distributions
also reflected the observed data (not shown).

Figure 5. Time series for surface chlorophyll-a, dissolved inorganic nitrate (DIN), dissolved inorganic

14

Figure 5. Time series
for surface chlorophyll-a, dissolved inorganic nitrate (DIN), dissolved inorganiccarbon (DIC) and total alkalinity (ALK) data for longitude 165◦ E from observed data and model outputs.
◦ N and 50◦ N, respectively.
Thetotal
purple, alkalinity
blue, green and(ALK)
red colorsdata
show thefor 20◦ N, 30◦ 165
N, 40°
carbon (DIC) andfordatalongitude
E from observed data and mode
outputs. The purple, blue, green and red colors show the data for 20° N, 30° N, 40° N and 50° N
respectively.

Sustainability 2019, 11, 2677

13 of 28

The concentrations of DIC (Figure 7a–d) were high (>2120 µmol/kg) in the Okhotsk Sea and thenorthwestern 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 seasonalvariation in the concentrations (Figure 7a–d). For the Okhotsk Sea and northwestern North Pacificarea 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 areaof 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 toautumn (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), weresimilar. However, unlike the observed data, there was a marked decrease (50 µmol/kg) for summer in
theFiguresubarctic
regionandOkhotskSea (Figure8f). Thelow values
summerwere followed
5. Timeseriesforthesurfacechlorophyll-a,
dissolved
inorganic
nitratein(DIN),
dissolved
inorganic byhighervaluesautumn(Figure 8g).
Thedata
timeforseriesfor ALK5g)observed
also showeda decreasecarbon(DIC)inand
total alkalinity(ALK)
longitude
165(Figure° E from
data andmodel in
◦ N and 50◦ N. The summer decrease arose from the modeled variations in CaCO
summeraround40
outputs. The purple, blue, green and red colors show the data for 20° N, 30° N, 40° N and 50° N, 3
(notshown), and indicates that the summer decrease in ALK may have been overestimated.
respectively.

Figure6. Surface
monthly
distributions(DIN:µmol/kg)
μmol/kg)ininJanuary
January(a,e),(a,e),
Figure6. Surface
monthly
distributionsofofdissolved
dissolvedinorganic
inorganic nitrate
nitrate (DIN:
March(b,f),
July(c,g),
andOctober
(d,h)
in
the
climatology
created
by
World
Ocean
Atlas
2013
March (b,f), July (c,g), and October (d,h) in the climatology created by World Ocean Atlas 2013 (WOA13)(WOA13)and Yasunaka
al. [55],
from themodelrespectively.
outputs, respectively.
and Yasunaka
et al. [55],etandfrom andthe model
outputs,

Sustainability 2019, 11, 2677
Sustainability
2019,
2677
Sustainability
2019,11,11,
2677

14 of 28
15 15
of of
30 30

Figure7. Surface
monthly
distributions
ofof
dissolved
inorganiccarbon(DIC:μmol/kg)in January(a,e),
Figure7. Surface
monthly
distributionsof dissolved
dissolved
inorganiccarbon(DIC:μmol/kg)in in
January(a,e),
Figure7.
Surfacemonthly
distributions
inorganiccarbon(DIC:µmol/kg)January(a,e),
April(b,f),
July(c,g),
andOctober(d,h)
fromclimatology(combined
datasets
from
Goyet
et
al.
[56],
April
(b,f),
July(c,g),
andOctober(d,h)
fromclimatology(combined
datasets
from
Goyet
et
al.
[56],
April (b,f), July (c,g), and October (d,h) from climatology (combined datasets from Goyet et al.
[56],
Key
etal.
al.
[57]
and
Yasunaka
et
[58];
see
Section
3.13.1
forfordetails)andmodel
outputs,
respectively.
Key
et al.
[57]
andYasunaka
et
[58];
see
Section3.1
fordetails)andmodel
outputs,
respectively.
Key
et
[57]
and
Yasunakaetal.
al.al.
[58];
see
Sectiondetails)andmodel
outputs,
respectively.

Figure 8. Plots of the JMA observed values for total alkalinity averaged above a depth of 20 m in spring
Figure
8. 8.
Plots
of of
thetheJMA
observed
valuesforfortotal
alkalinity
averagedabove
a depth
of of
20 20
mm
in in
FigurePlots
JMAobserved
valuestotal
alkalinity
averagedabove
a depth(a; March to May), summer (b; June to August) and autumn (c; September to November), and surface
spring
(a;
March
to
May),
summer(b;
JunetoAugust)andautumn(c;
SeptembertoNovember),andspring (a; March to May), summer (b; June to August) and autumn (c; September to November), andmonthly
distributionsof totalof
alkalinity(ALK: (ALK:µmol/kg)in January(d), April(e), July(f) andOctobersurfacemonthly
distributions
total
alkalinity
μmol/kg)in January(d),(d),April(e),(e),July(f) (f)andandsurfacemonthly
distributionsof
total
alkalinity(ALK:
μmol/kg)in January
April
July
(g)
from
model
outputs.
The
color
bar
indicates
thealkalinity
concentration.
October (g) from model outputs. The color bar indicates the alkalinity concentration.
October (g) from model outputs. The color bar indicates the alkalinity concentration.

4.2. Ocean Acidification Indices pH and Ωarg

4.2.
Ocean
Acidification
Indices
pHpH
andand
Ωarg
4.2.
Ocean
Acidification
IndicesΩarg

Figure 9 shows the surface pHinsitu (pH25 ) and Ωarg values calculated from the model outputs
Figure9 shows
thethesurface
pHpH
insitu(pH25) 25andΩargvaluescalculated
fromthethemodel
outputs
forforFigure9 shows
surface
insitu
(pH
) and
Ωarg
valuescalculated
from
model
outputs
for temperature, salinity, DIC and ALK (Figure 9a–d). The pHinsitu ranged from 7.85 to 8.10 in thetemperature,
DIC
and
ALK
(Figure9a–d).
The
pHpH
insituranged
from
7.85
to to
8.10
in in
thethetemperature,salinity,
salinity,
DIC
and
ALK
(Figure9a–d).
The
insituranged
from
7.85
8.10
northwestern
Pacificareaarea(Figure(Figure9a–d),
exceptin winterand spring(January
andforApril)for thenorthwestern
Pacific9a–d),
exceptin in
winterandspring(January
andApril)
thetheEastnorthwestern
Pacific area
(Figure9a–d),
exceptwinterandspring(January
andApril)
for
East
East
China
Sea
near
the
coast
of
China
(Figure9a,b).
The
values
in
the
side
of
Pacific
basically
ranged
China
Sea
near
the
coast
of
China
(Figure9a,b).
The
values
in
the
side
of
Pacific
basically
rangedfrom
China Sea near the coast of China (Figure 9a,b). The values in the side of Pacific basically ranged from
from
8.00
to
8.05.
The
valueswerelower(7.85–7.90)
inOkhotsk
the Okhotsk
Sea
and
the
Japan
Sea throughout
8.00
to to
8.05.
The
valueswerelower(7.85–7.90)
in inthetheOkhotsk
SeaSeaandthetheJapan
Seathroughout
thethe
8.00
8.05.
The
valueswerelower(7.85–7.90)
and
Japan
Sea
throughout
the
year,
and
on
the
Pacific
Ocean
side
along
the
Kamchatka
Peninsulaand
Kuril
Islands
fromwinter
year,
and
onon
thethePacific
Ocean
side
along
thetheKamchatka
Peninsulaand
Kuril
Islands
fromwintertoyear,
andPacific
Ocean
side
along
Kamchatka
Peninsulaand
Kuril
Islands
fromwintertoto 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)Sustainability 2019, 11, 2677

15 of 28

(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 pH25 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 pH25 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 thesubarctic region (Figure 9i–l). The values were lowest (1.0–1.5) in the Okhotsk Sea and along theKamchatka 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 pH25 (Figure 10b).
Takahashi et al. [59] and Jiang et al. [64] created a global map of the surface climatology for pHinsitu
and Ωarg in spring/summer and autumn/winter, but not for the Okhotsk, Japan and East China Seas.
Their climatology showed that the pHinsitu and the Ωarg varied seasonally from 0.06 to 0.08 and from0.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 theconditions in the 2010s (their Figure 7b), were similar to our model outputs (Figure 9i–l).
Sustainability 2019, 11, 2677

17 of 30

Figure9. 9.Horizontal
surfacedistributions
distributions
ofinsitu
pH(a–d),(a–d),
pHand(e–h),andsaturation
aragonite(Ωsaturation
insitu pH
FigureHorizontal surfaceof pH
25 (e–h), 25
aragonite
arg;
(Ωarg
in January,
July,and Octobertheoutputs.
model outputs.
i–l); i–l)in January,
April,April,
July, and
Octoberfrom thefrom
modelFigure9. Horizontal surface distributions of pHinsitu (a–d), pH25 (e–h), and aragonite saturation(Ωarg;
Sustainability 2019, 11, 2677
16 of 28i–l) in January, April, July, and October from the model outputs.

FigureTimeofseriesof surface
pH(a),(a),25pH(b) andat longitude165°
165 E
E fromarg (c)insitupH
Figure 10.
Time10.series
surface
pHinsitu(b)25 andΩargΩ(c)at longitude
fromthethe model
model outputs. The purple, blue, green and red colors show the data for 20◦ N, 30◦ N, 40◦ N andoutputs. The
purple, blue, green and red colors show the data for 20° N, 30° N, 40° N and 50° N,
50◦ N, respectively.
respectively.
◦

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, pH25 and ALK (Table 3). The observed data andmodel 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
pH25 ; 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 (Figures 11–13) showed that theparameters were highly correlated with each other (not shown for DIP, DIC and ALK).
Table 3. Correlation coefficients for chlorophyll-a (Chl-a), DIN, DIC, ALK, pH25 and Ωarg between
JMA observed data and the model outputs for each region in 2015 (Figure 1). The observed data andmodel 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, log10 (Chl-a) values
were evaluated.
Parameter

Subtropical
RegionSubarctic
RegionKuroshio
Extension

Japan Sea

Chl-a
DIN
DIC
ALK
pH25
Ωarg

0.57 (0.07)
0.72 (0.04)
0.74 (0.09)
0.48 (0.09)
0.78 (0.09)
0.99 (0.09)

0.88 (0.07)
0.88 (0.12)
0.87 (0.20)
0.69 (0.20)
0.86 (0.20)
0.98 (0.20)

0.86 (0.12)
0.87 (0.10)
0.78 (0.17)
0.48 (0.17)
0.81 (0.17)
0.98 (0.17)

0.87 (0.15)
0.95 (0.11)
0.92 (0.26)
0.35 (0.26)
0.90 (0.26)
0.99 (0.26)

significance level of 0.05 and the significant level of p < 0.05. Only for Chl-a, log10 (Chl-a) values were
evaluated.
Parameter
Chl-a
DIN
DIC
Sustainability 2019, 11, 2677 ALK
pH25
Ωarg

Subtropical Region
0.57 (0.07)
0.72 (0.04)
0.74 (0.09)
0.48 (0.09)
0.78 (0.09)
0.99 (0.09)

Subarctic Region
0.88 (0.07)
0.88 (0.12)
0.87 (0.20)
0.69 (0.20)
0.86 (0.20)
0.98 (0.20)

Kuroshio Extension
0.86 (0.12)
0.87 (0.10)
0.78 (0.17)
0.48 (0.17)
0.81 (0.17)
0.98 (0.17)

Japan Sea
0.87 (0.15)
0.95 (0.11)
0.92 (0.26)
0.35 (0.26)
0.90 (0.26)
0.99 (0.26)

17 of 28

Figure 11. Scatter plots showing the chlorophyll-a concentrations for the subtropical (a) and subarctic

Figure 11. Scatterplots showing the chlorophyll-a concentrations for the subtropical (a) and subarctic(b) regions, the Kuroshio-extension (c), and the Japan Sea (d). The horizontal and vertical axes denote(b) regions, theKuroshio-extension(c),values,
and respectively.
the JapanTheSea(d).indicate
The thehorizontal
and vertical axes
the simulated and JMA observedcolorsobserved and/or
simulated
depths.
denote the simulated and JMA observed values, respectively. The colors indicate the observed and/or
simulatedSustainability
depths. 2019, 11, 2677
19 of 30
Sustainability 2019, 11, 2677

19 of 30

Figure 12. Scatter plots showing the pH25 values for the subtropical (a) and subarctic (b) regions, theFigure 12. Scatterplots showing the pH25 (d).
valuesfor the subtropical (a) and subarctic (b) regions,
Kuroshio-extension (c) and the Japan Sea
The horizontal and vertical axes denote simulated andFigure 12. Scatter plots showing the pH25 values for the subtropical (a) and subarctic (b) regions, thethe Kuroshio-extension
(c)
and
the
Japan
Sea
(d).
The
horizontal
andsimulated
verticaldepths.
axes denote simulated
JMA observed values, respectively. The colors indicate theobserved and/or
Kuroshio-extension (c) and the Japan Sea (d). The horizontal and vertical axes denote simulated andand JMA observed
values,respectively.
The
colorsindicate
the observed
and/orsimulated depths.
JMA observed
values,respectively. The
colorsindicate
the observed
and/or simulated
depths.

.
Figure 13. Scatter .plots showing Ωarg values for the subtropical (a) and subarctic (b) regions, theKuroshio-extension (c) and the Japan Sea (d). The horizontal and vertical axes denote the simulated
Figure 13.plots
Scatter plots showing
arg values for the subtropical (a) and subarctic (b) regions, theFigure 13. ScatterΩΩthethesubtropical(a)and subarctic
arg values
and JMA observedshowing
values, respectively.
The colorsfor
indicate
observed and/or
simulated
depths. As (b) regions,
Kuroshio-extension (c) and the Japan Sea (d). The horizontal and vertical axes denote the simulated
theobserved
values
for
Ω
arg
were
not
providedin(d).
the downloaded
JMA data, weandcalculated
the values
the Kuroshio-extension
(c)
and
the
Japan
Sea
The
horizontal
verticalaxes denote theand JMA observed values, respectively. The colors indicate the observed and/or simulated depths. As
for Ωarg from the observed JMA values for temperature, salinity, DIC and ALK.
simulated andobserved
respectively.
colors JMA
indicatetheobserved
and/or simulated
theJMAobserved
values forvalues,Ωarg were notprovided in theThe
downloaded
data, wecalculated
the values
forobservedΩarg from theobserved
JMA
for temperature,salinity,
ALK.
depths. As thevaluesforΩ values
were
not provided
in DIC
theanddownloaded
JMA data, we calculated

5. Discussionarg

the values5. forΩarg from the observed JMA values for temperature, salinity, DIC and ALK.
Discussion
5.1. Sensitivity Experiment for pH and Ωarg Based on Constant Temperature, Salinity, DIC and ALK

To assess Experiment
which variables
importantin Constantdetermining
the oceanSalinity,
acidification
indices
5.1. Sensitivity
for pHwereand Ω
arg Based on
Temperature,
DIC and
ALK pH andΩarg, we performed a series of sensitivity experiments assuming a constant set of conditions for
To assess which variables were important in determining the ocean acidification indices pH andtemperature (= 21 C), salinity (= 34.44), DIC (= 1995.2 μmol/kg) and ALK (= 2260.9 μmol/kg). In theΩarg, we performed a series of sensitivity experiments assuming a constant set of conditions forsensitivity experiments, pH and Ωarg were calculated with the assumed fixed variables oftemperature (= 21 C), salinity (= 34.44), DIC (= 1995.2 μmol/kg) and ALK (= 2260.9 μmol/kg). In thetemperature, salinity, DIC and ALK, individually, in place of these original target variables of thesensitivity experiments, pH and Ωarg were calculated with the assumed fixed variables ofmodel outputs. The other variables, except the assumed constant variables, were the same as thoseSustainability 2019, 11, 2677

18 of 28

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 pHand Ωarg , we performed a series of sensitivity experiments assuming a constant set of conditions fortemperature (= 21 ◦ C), salinity (= 34.44), DIC (= 1995.2 µmol/kg) and ALK (= 2260.9 µmol/kg). In thesensitivity 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 determinedby 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 thedifferences enabled the individual effects of each target variable to be assessed. The temperature and
DIC had most influence on the pHinsitu 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 oftemperature and salinity, while detectable, were small compared with the effects of DIC and ALK(Figure 15e–h).
5.2. Effect of Seasonal Processes on pHinsitu
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 tothe dissociation equilibrium: H2 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: H2 O←H+ +OH− (H2 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 thepH 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 associatedwith the CO2 equilibrium in sea water, defined by: CO2 +H2 O↔H+ +HCO3 − . 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: CO2 +H2 O→H+ +HCO3 − (CO2 +H2 O←H+ +HCO3 − ). The decreasein 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: CO2 +H2 O←H+ +HCO3 − .
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, andthe 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.

Sustainability 2019, 11, 2677

19 of 28

From the sensitivity analyses, we inferred that overall seasonal variations in the pHinsitu(Figure 9a–d and Figure 10a) were largely governed by the balance between temperature and DIC, andto a lesser extent by ALK. In the southern region, the relationship with temperature was relatively
important, and thus the pHinsitu decreased from winter to summer. In contrast, the pHinsitu increasedslightly from winter to summer in the northern region, where the relationship with DIC and ALK
was slightly stronger. The pH25 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 insummer, as shown in Figure 8b,f, and the variations in the pHinsitu may actually be mainly governedby temperature
and DIC.
Sustainability 2019, 11, 2677
21 of 30

Figure 14. Surface distributions of the pH differences (ΔpH) in January, April, July, and October from

Figure 14. Surface distributions of the pH differences (∆pH) in January, April, July, and October fromsensitivity experiments using fixed values for temperature (a−d), salinity (e−h), DIC (i−l), and ALK
sensitivity
using
valuesfor temperature(a–d),
salinity(e–h),
DICthan(i–l),(m–p).experiments
When calculating
the fixed
pH in thesensitivity
experiments, thevaluesof variables
other
theand ALK(m–p).fixedWhencalculatingin thesensitivity
experiments,
the valuesof variables
other thanvariables
were the thesamepHas thosein themodel outputs
shown in Figure
4e–h, Figure
6e–h, Figure
the fixed
those inthesubtracting
model outputs
shown
infromFigures
4e–h, 6e–h, 7e–h
7e–hvariables
and Figurewere
8d–g.theΔpHsameis the asremainderafter
the model
outputs
the sensitivity
experimentfrombasic model
outputs
for each month.
and 8d–g.
∆pH isthetheremainder
after
subtracting
the model outputs from the sensitivity experimentfrom the basic model outputs for each month.
5.3. Seasonal Processes Affecting Ωarg Variations

5.3. Seasonal
Affectingarg Variations
TheProcesses
sensitivityof the Ωtargetvariables to Ωarg are shown in Figure 15. The temperaturedependence of Ωarg can be explained by the characteristics of the aragonite solubility product, karg (arg

The sensitivityof the target variables to Ωarg are shown in Figure 15. The temperature dependence= ([Ca2+][CO32−])/karg). The value of karg decreases as the water temperature increases [34]. Latitudinalof Ωargcanbeexplained
by of
thetemperature
characteristics
aragonitekarg (Ωarg =differences in the sensitivity
resultedofin thehigherand lowersolubility
Ωarg values product,in the southern2− ])/k ). The value of k([Ca2+ ][COdecreases
as
the
water
temperature
increases[34].
Latitudinalarg
arg respectively. Seasonal warming/cooling affected the seasonaland northern
regions (Figure 15a–d),
3
variationin Ωsensitivity
arg throughout
entire region.resulted in higher and lower Ωarg values in the southerndifferencesin theofthe
temperature
The reaction
the dissolution
marine carbonate[3,6], CO
2+H2O+CO3−↔2HCOaffected3−, was related
and northern
regionsfor(Figure15a–d),ofrespectively.
Seasonal
warming/coolingthe seasonal
to latitudinal differences in the sensitivity of DIC (Figure 15i–l). As the DIC concentration was
variation in Ωarg throughout the entire region.
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 theentire 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

Sustainability 2019, 11, 2677

20 of 28

The reaction for the dissolution of marine carbonate [3,6], CO2 +H2 O+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 a decrease in
ALK in the northern region from winter to summer (Figure 8a–c), the Ωarg values tended to decrease
Sustainability 2019, 11, 2677
22 of 30
(Figure 15m–o). The sensitivity of salinity to Ωarg (Figure 15e–h) was low, and similar to that of the
ALK in the northern region from winter to summer (Figure 8a–c), the Ωarg values tended
pHinsitua decrease(Figure in
14e–h).
todecrease(Figure15m–o).
The sensitivityof salinity
to Ωarg(Figurewas low,and similar
From the sensitivity
analysis,
we inferredthat DIC
andALK15e–h)
had most
impacton thetoseasonalthat of the pHinsitu (Figure 14e–h).
variability in Ωarg across the entire region throughout the year, with less of a temperature effect.
From the sensitivity analysis, we inferred that DIC and ALK had most impact on the seasonal
However,the summer increase in Ω region(Figure10c) was mainly the result of the decrease in DIC invariability in Ωarg across the entirearg
throughout the year, with less of a temperature effect.
summer(Figurethe
9i–l
and Figure
becauseprobably
overestimateddecrease
However,summerincrease10c),in Ωarg(Figure the10c)model
was mainly
the resultof the decreasea in
DIC in in ALK
in summerandthe decreasein ALK
the increase
temperaturein summer
togethersummer(Figure9i–l and Figure10c),andbecausethe modelinprobably
overestimateda decreasein ALKcaused ain summer
decreasein Ωarg .and the decrease in ALK and the increase in temperature in summer together caused adecrease in Ωarg.

15. Surface
distributionsof differencesininthe
thearagonite
aragonite saturation
saturation (ΔΩarg) in January, April,
FigureFigure
15. Surface
distributionsof differences(∆Ωarg ) in January, April, July,
July, and October from sensitivity experiments using fixed values for temperature (a–d), salinity (e–and October from sensitivity experiments using fixed values for temperature (a–d), salinity (e–h), DICh), DIC (i–l), and ALK (m–p). When calculating the Ωarg in the sensitivity experiments, the values of(i–l), andALK (m–p). When calculating the Ωargsamein thesensitivity experiments, the values of variables
variables other than the fixed variables were theas those in the model outputs shown in Figure
other than
fixedvariables
outputs
shown in Figures
4e–h, 6e–h,
4e–h,theFigure6e–h,Figure were
7e–h theand same
Figureas8.thoseΔΩarg inis thethe model
remainder
after subtracting
the model
7e–h and8. ∆Ω
is theremainder
after subtracting
the
model
outputsfromthe sensitivity experiment
outputsfromsensitivity
experimentfrom the basic
model
outputs
for eachmonth.
arg thefrom the basic model outputs for each month.

6. Conclusions
We developed a new biogeochemical model (JCOPE_EC) for carbon processes, based on thephysical 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Sustainability 2019, 11, 2677

21 of 28

6. Conclusions
We developed a new biogeochemical model (JCOPE_EC) for carbon processes, based on thephysical 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 seasonalvariability in these biogeochemical variables was similar to the observed (climatological) variability,
although the modelled seasonal variability in chlorophyll-a and ALK deviated somewhat fromthe 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 CaCO3 and also ALK reduction (Figures 5g and8f). 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 pHinsitu and Ωarg values in some parts in the Okhotsk Sea and the Japan Sea fromwinter to spring were lower than expected (Figure 9a–d,i–l). This is probably a consequence of theuncertainty 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 pHinsitu were largely governed by the correspondingbalances 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 biogeochemicalvariables at the basin scale. However, the JCOPE_EC was designed to use operational physical datafrom 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 spring2015 in Figures 16 and 17 show that the JCOPE_EC was able to represent the distribution ofbiogeochemical 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 todownstream regions. Prominent areas of enhanced chlorophyll-a concentrations were evident, andtheir 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 associatedpH 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 performa nowcast/forecast experiment targeting ocean acidification indices (pHinsitu 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 somebiogeochemical variables, which highlights that the present model is diagnostic. To investigate themechanisms 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 andthe 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.

Sustainability 2019, 11, 2677
Sustainability
11, 2677
Sustainability
2019,2019,11, 2677

24 of 30
22 of 28
24 of 30

Figure 16. Reproduced snapshots of chlorophyll-a concentrations with current velocity (a) and DIN
Figure16. Reproduced
snapshots
concentrationswithcurrent
velocity(a)DINand DIN
Figure16. Reproduced
snapshotsofofchlorophyll-achlorophyll-a concentrationswithcurrent
velocity(a) and(b) on 1 April.
1 April.
(b) on(b)1 on
April.

Figure 17. Surface distributions of DIC (μmol/kg) (a), ALK (μmol/kg) (b), pHinsitu (c) and aragonite
Figure17. Surfaceof DIC(μmol/kg)(µmol/kg) (a),(µmol/kg)(b), pHand(c)
aragonite
Figure17.
Surface
distributions(a),ALK
ALK(μmol/kg)(b),insitu
pH(c)insituand aragonite
saturation(Ωarg;; d)ondistributions
1 April. of DIC
saturation Ωarg; ; (d) on 1 April.
saturation (Ωarg;; d) on 1 April.

Author Contributions: Conceptualization, T.T. and Y.M.; Methodology, Y.M. and M.I; Validation, M.I.; Formal
Author Contributions: Conceptualization, T.T. and Y.M.; Methodology, Y.M. and M.I; Validation, M.I.; Formal

analysis,
M.I.;M.I.;
Writing—Original
Draft
Preparation,
M.I.;
Writing—Review
andM.I;
Editing,
M.I. and
AuthorContributions:Conceptualization,
T.T.
andM.I.;
Y.M.;
Methodology,
andValidation,
M.I.;Y.M.;
Formalanalysis,
Writing—Original
Draft
Preparation,
Writing—Review
and Y.M.
Editing,
M.I.
and
Y.M.;
Visualization,
Visualization,
M.I.; Supervision,
Y.M.
andPreparation,
X.G.; Project administration,
T.T.; Funding
Acquisition,
T.T.
M.I.;
Supervision,
Y.M. and X.G.;
Project
administration,
T.T.; Funding
Acquisition,
T.T.andand.
Y.M.
analysis,
M.I.;
Writing—Original
Draft
M.I.;
Writing—Review
Editing,
M.I. and.
andY.M.
Y.M.;
Funding:This
research was
fundedby theSasakawa
Peace Foundationof the Ocean
Policy Research
Visualization,
M.I.;research
Supervision,
Y.M.
and
X.G.;
Project
administration,
T.T.;
Funding
T.T.Instituteand. Y.M.
Funding:This
was fundedby theSasakawa
Peace
Foundationof theOcean Acquisition,
PolicyResearch
Institute (OPRI-SPF).

(OPRI-SPF).
Funding:
This research was funded by the Sasakawa Peace Foundation of the Ocean Policy Research Institute(OPRI-SPF).
Acknowledgments: We acknowledge the consistent support from the Sasakawa Peace Foundation of the Ocean
Policy Research Institute (OPRI-SPF). This study is part of the Japan Coastal Ocean Predictability Experiment

Acknowledgments: We acknowledge the consistent support from the Sasakawa Peace Foundation of the Ocean(JCOPE) promoted by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC). NCEP/NCAR
Policy Research Institute (OPRI-SPF). This study is part of the Japan Coastal Ocean Predictability Experiment
reanalysis data were downloaded from the websites of National Centers for Atmospheric Research:(JCOPE) promoted by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC). NCEP/NCAR
https://climatedataguide.ucar.edu/climate-data/ncep-ncar-r1-overview. World Ocean Atlas 2001 (WOA01) andreanalysis
data were
downloaded
fromwebsitesCenters for AtmosphericResearch:2013 (WOA13)
data were
downloaded
fromthewebsitesof theofUSNational
National Oceanographic
Data Center (NODC):
https://climatedataguide.ucar.edu/climate-data/ncep-ncar-r1-overview.
World Ocean Atlas 2001 (WOA01) andhttps://www.nodc.noaa.gov/OC5/WOA01/pr_woa01.html
andhttps://www.nodc.noa.gov/OC5/woa13/,

Sustainability 2019, 11, 2677

23 of 28

Acknowledgments: We acknowledge the consistent support from the Sasakawa Peace Foundation of the
Ocean Policy Research Institute (OPRI-SPF). This study is part of the Japan Coastal Ocean Predictability
Experiment (JCOPE) promoted by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC).
NCEP/NCAR reanalysis data were downloaded from the websites of National Centers for Atmospheric Research:https://climatedataguide.ucar.edu/climate-data/ncep-ncar-r1-overview. World Ocean Atlas 2001 (WOA01) and2013 (WOA13) data were downloaded from websites of the US National Oceanographic Data Center (NODC): https://www.nodc.noaa.gov/OC5/WOA01/pr_woa01.html and https://www.nodc.noa.gov/OC5/woa13/, respectively.
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 forcarbon 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 datafor 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 cdom and Iopt 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 preliminarymodel at which chlorophyll-a reached a maximum in the subtropical region (10◦ –20◦ N) were muchshallower (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 cdom and Iopt that varied by latitude (Equations (39) and (40)) (Figure 3a,b).
Using this approach, the model reproduced the vertical distribution and maximum depthof 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 dependencesof cdom and Iopt were nonlinear, we manually adjusted the parameter values after comparing themodeled and observational results.
Previous NPZD models of the marginal seas of Japan have used various constant Iopt values.
For example, Kawamiya et al. [19–21] used an Iopt value of 48.8 W/m2 for station Papa (50.1◦ N,
144.9◦ W) [19,20] and the North Pacific [21], Yamanaka et al. [48] used an Iopt value of 104.7 W/m2
for the North Western Pacific, Sasai et al. [25] used an Iopt value of 100 W/m2 for stations S1 and K2,
Kishi et al. [22] used an Iopt value of 104.7 W/m2 for the subarctic region of the North Pacific, Guo andYanagi et al. [27] used an Iopt value of 150.0 W/m2 for the East China Sea, and Onitsuka and Yanagi [26]used an Iopt value of 70.0 W/m2 .
Edwards et al. [61] characterized growth-irradiance relationships from information about
308 growth-irradiance experiments performed on 119 species of marine phytoplankton sampledfrom all major functional groups, and reported that the Iopt values did not vary among taxonomic
groups. However, all traits varied considerably within most groups, and optic-ocean isolates tended
to have lower Iopt 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 Iopt .
Various studies have used constant values ranging from 0.035 to 0.04 m−1 for cdom in the North
Pacific [19–22,25,48]. These cdom values represent light dissipation coefficients, which are linked tothe distribution of chromophoric dissolved organic matter. Nelson and Sigel [63] showed that theabsorption 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 thesubarctic region. This indicates that the distribution of chromophoric dissolved organic matter makemore invisible the subsurface layer in the subarctic region than in the subtropical region. The latitudinaldifference for cdom used in the present study represents a change from 0.015 m−1 in the subtropicalSustainability 2019, 11, 2677

24 of 28

region to 0.045 m−1 in the subarctic region (Figure 3a), which is consistent with observed distributions
26 of 30of chromophoric dissolved organic matter [63].

Sustainability 2019, 11, 2677

3) along longitude
3 ) along
FigureA1. Vertical
Verticaldistribution
distributionof chlorophyll-a
165°E
from JMA
Figure A1.
of chlorophyll-a(mg/m(mg/mlongitude 165◦ E from
JMA observations
observations
in
2015
(a)
and
model
outputs(b−e).
The
figureshows
thesensitivity
experimentresults
in 2015 (a) and model outputs (b–e). The figure shows the sensitivity experiment results with (b) constant
2
−1
2
2
−1
2
with
(b)70constantIopt c(=dom70(=W/m
domconstant(= 0.04 m
Ioptlatitudinal(= 70 W/mdifferences) and latitudinal
Iopt (=W/m ) and0.04 )mand); c(c)Iopt );(=(c)
70constantW/m ) andfor cdom ;differences
for cdifferences
dom; (d) latitudinaldifferences
Iopt(=andconstant
cdom(e)(= latitudinal0.04 m−1); and(e) latitudinal(d) latitudinal
for Iopt andconstantfor
cdom
0.04
m−1 ); anddifferencesfor Iopt
differences
forresults
Iopt andofcdom
resultsof theas 4e–h).
shown in Figure 4e–h).
and cdom (thethe(themodel
outputs
asmodelshownoutputs
in Figure

References
1.

Baur, J.F.; Cai, W.J.; Raymond, P.A.; Bianchi, T.S.; Hopkinson, C.H.; Pierre, A.; Regnier, P. The changing
carbon cycle of the coastal ocean. Nature 2013, 504, 61–70, doi:10.1038/nature12857.

Sustainability 2019, 11, 2677

25 of 28

References
1.
2.

3.

4.

5.

6.

7.
8.

9.

10.

11.

12.
13.
14.
15.

16.
17.
18.

19.

Baur, J.F.; Cai, W.J.; Raymond, P.A.; Bianchi, T.S.; Hopkinson, C.H.; Pierre, A.; Regnier, P. The changing
carbon cycle of the coastal ocean. Nature 2013, 504, 61–70. [CrossRef]
Bednarsek, N.; Feely, A.R.; Reum, J.C.P.; Peterson, B.; Menkel, J.; Alin, R.S.; Hales, B. Limacina helicina shelldisoolution as an indicator of declining habitat suitability owing to ocean acidification in the California
Current Ecosystem. Proc. Biol. Sci. 2014, 281, 20140123. [CrossRef] [PubMed]
Feely, A.R.; Sabine, L.C.; Barne, H.R.; Millero, J.F.; Dickson, G.A.; Wanninkhof, R.; Murata, A.; Miller, A.L.;
Greeley, D. Decadal changes in the aragonite and calcite saturation state of the Pacific Ocean. Glob. Biogeochem.
Cycles 2012, 26, GB3001. [CrossRef]
Stocker, T.F.; Qin, D.; Plattner, G.-K.; Tignor, M.M.B.; Allen, S.K.; Boschung, J.; Nauels, A.; Zia, Y.B.;
Midgley, P.M. (Eds.) Intergovernmental Panel on Climate Change (IPCC), Climate Change 2013: The PhysicalScience Basis. In Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on
Climate Change; Cambridge University Press: Cambridge, UK; New York, NY, USA, 2013; pp. 1–1535.
Kroeker, J.K.; Kordas, L.R.; Crim, R.; Hendrik, E.I.; Ramajo, L.; Singh, S.G.; Duarte, M.C.; Gattuso, P.J. Impactsof ocean acidification on marine organisms: Qualitifying sensitivities and interaction with warming. Chang.
Biol. 2013, 19, 1884–1896. [CrossRef] [PubMed]
Feely, A.R.; Alin, R.S.; Carter, B.; Bednarsek, N.; Hales, B.; Chan, F.; Hill, M.T.; Gaylord, B.; Sanford, E.;
Byrne, H.R.; et al. Chemical and biological impacts of ocean acidification along the west coast of North
America. Estuar. Coast. Shelf Sci. 2016, 183, 260–270. [CrossRef]
Wittmann, C.A.; Portner, H. Sensitivities of extant animal taxa to ocean acidification. Nat. Clim. Chang. 2013,
3, 995–1001. [CrossRef]
IGBP; IOC; SCOR. Ocean Acidification Summary for Policymakers. In Proceedings of the Third symposiumon the Ocean in a High-CO2 World International Geoshpere-Bioshere Programme, Stockholm, Sweden,
14 November 2013.
Bender, M.; Doney, S.; Feely, R.A.; Fung, I.; Gruber, N.; Harrison, D.E.; Keeling, R.; Moore, J.K.; Sarmiento, J.;
Sarachik, E.; et al. A Large–Scale CO2 Observation Plan: In Situ Oceans and Atmosphere (LSCOP); NOAA OAR
Special Report; National Oceanic and Atmospheric Administration: Silver Spring, MD, USA, 2002; pp. 1–201.
Johnson, S.K.; Plant, N.J.; Coletti, J.L.; Jannasch, W.W.; Sakamoto, M.C.; Rise, C.S.; Swift, D.D.; Williams, L.N.;
Boss, E.; Haentjens, N.; et al. Biogeochemical sensor performance in the SOCCOM profiling float array.
J. Geophys. Res. Ocean 2017, 122, 6416–6436. [CrossRef]
Doney, S.C.; Lima, I.; Feely, R.A.; Glover, D.M.; Lindsay, K.; Mahowald, N.; Moore, J.K.; Wanninkhof, R.
Mechanisms governing interannual variability in upper-ocean inorganic carbon system and air-sea CO2
fluxes: Physical climate and atmosphere dust. Deep-Sea Res. II 2009, 56, 640–655. [CrossRef]
Metz, N.; Tilbrook, B.; Bakker, D.; Quere, L.; Doney, C.; Feely, R.; Hood, M.; Dargaville, R. Global changes inocean carbon: Variability and vulnerability. EOS Trans. 2007, 88, 287.
Le Quere, C.; Orr, J.C.; Monfray, P.; Aumont, O. Interannual variability of the oceanic sink of CO2 from 1979
through 1997. Glob. Biogeochem. Cycles 2000, 14, 1247–1265. [CrossRef]
McKinley, G.A.; Follows, M.J.; Marshall, J. Mechanisms of air–Sea CO2 flux variability in the equatorial
Pacific and the North Atlantic. Glob. Biogeochem. Cycles 2004, 18, GB2011. [CrossRef]
McKinley, G.A.; Takahashi, T.; Butenhuis, E.; Chai, F.; Christian, J.R.; Doney, S.C.; Jiang, M.S.; Lindsay, K.;
Moore, J.K.; LeQuere, C.; et al. North Pacific carbon cycle response to climate variability on seasonal todecadal timescales. J. Geophys. Res. Oceans 2006, 111, C07S06. [CrossRef]
Obata, A.; Kitamura, Y. Interannual variability of the sea-air exchange of CO2 from 1961 to 1998 simulatedwith a global ocean circulation-biogeochemistry model. J. Geophys. Res. 2003, 108, 3337. [CrossRef]
Wetzel, P.; Winguth, A.; Maier-Reimer, E. Sea–to–air CO2 flux from 1948 to 2003: A model study. Glob.
Biogeochem. Cycles 2005, 19. [CrossRef]
Schmittner, A.; Oschiles, A.; Matthews, H.D.; Galbraith, E.D. Future changes in climate, ocean circulation,
ecosystems, and biogeochemical cycling simulated for a business-as-usual CO2 emission scenario until year
4000 AD. Glob. Biogeochem. Cycles 2008, 22, GB1013. [CrossRef]
Kawamiya, M.; Kishi, M.; Yamanaka, Y.; Suginohara, N. An ecological-physical coupled model applied tostation Papa. J. Oceanogr. 1995, 51, 635–664. [CrossRef]

Sustainability 2019, 11, 2677

20.
21.

22.

23.

24.
25.

26.
27.
28.
29.
30.

31.

32.
33.
34.

35.

36.

37.
38.

39.

26 of 28

Kawamiya, M.; Kishi, M.J.; Yamanaka, Y.; Suginohara, N. Obtaining reasonable results in different oceanic
regimes with the same ecological-physical coupled model. J. Oceanogr. 1997, 53, 397–402.
Kawamiya, M.; Kishi, M.; Suginohara, N. An ecosystem model for the North Pacific embedded in a generalcirculation model Part I: Model description and characteristics of spatial distributions of biological variables.
J. Mar. Syst. 2000, 25, 129–157. [CrossRef]
Kishi, M.J.; Kashiwai, M.; Ware, D.M.; Megrey, B.A.; Eslinger, D.L.; Werner, F.E.; Noguchi-Aita, M.;
Azumaya, T.; Fujii, M.; Hashimot, S.; et al. NEMURO—A lower trophic level model for the North Pacificmarine ecosystem. Ecol. Model. 2007, 202, 12–25. [CrossRef]
Fujii, M.; Yamanak, Y.; Nojiri, Y.; Kishi, J.M.; Chai, F. Comparison of seasonal characteristics in biogeochemistry
among the subarctic North Pacific stations described with a NEMURO-based marine ecosystem model. Ecol.
Model. 2007, 202, 52–67. [CrossRef]
Yoshie, N.; Yamanaka, Y.; Rose, K.A.; Eslinger, D.L.; Ware, D.M.; Kishi, M.J. Parameter sensitivity study of
the NEMURO lower trophic level marine ecosystem model. Ecol. Model. 2007, 202, 26–37. [CrossRef]
Sasai, Y.; Yoshikawa, C.; Smith, L.S.; Hashioka, T.; Matsumoto, K.; Wakita, M.; Sasaoka, K.; Honda, C.M.
Coupled 1-D physical-biological model study of phytoplankton production at two contrasting time-series
stations in the western North Pacific. J. Oceanogr. 2016, 72, 509–526. [CrossRef]
Onitsuka, G.; Yanagi, T. Difference in ecosystem dynamics between the northern and southern parts of theJapan Sea: Analyses with two ecosystem models. J. Oceanogr. 2005, 61, 415–433. [CrossRef]
Guo, X.; Yangai, T. The role of the Taiwan Strait in an ecological model in the east china sea. Acta Oceanogr.
Taiwanica 1998, 37, 139–162.
Zhao, L.; Gao, X. Influence of cross-shelf water transport on nutrients and phytoplankton in the East China
Sea: A model study. Ocean Sci. 2011, 7, 27–43. [CrossRef]
Luo, X.; Wei, H.; Liu, Z.; Zhao, L. Seasonal variability of air-sea CO2 fluxes in the Yellow and East China Seas:
A case study of continental shelf sea carbon cycle model. Cont. Shelf Res. 2015, 107, 69–78. [CrossRef]
Bopp, L.; Resplandy, L.; Doney, S.C.; Dunner, J.P.; Gehlen, M.; Halloran, P.; Heinze, C.; Ilyinas, T.; Seferian, R.;
Tjiputra, J.; et al. Multiple stressors of ocean ecosystems in the 21th century: Projections with CMIP5 models.
Biogeosciences 2013, 10, 6225–6245. [CrossRef]
Moss, R.H.; Nakicenovic, N.; O’Neill, B.C. Towards New Scenarios for Analysis of Emissions, Climate
Change, Impacts, and Response Strategies. In Proceedings of the Intergovernmental Panel on Climate
Change, Geneva, Switzerland, 2 September 2008; p. 132.
Taylor, E.K.; Stouffer, J.R.; Meehl, A.G. A Summary of the CMIP5 Experiment Design. Available online:
https://pcmdi.llnl.gov/mips/cmip5/experiment_design.html (accessed on 1 February 2019).
Taylor, E.K.; Stouffer, J.; Meehl, A.G. An overview of CMIP5 and the experiment design. Bull. Am. Meteorol.
Soc. 2012, 93, 485–498. [CrossRef]
Watanabe, S.; Hajime, T.; Sudo, K.; Nagashima, T.; Takemura, T.; Okajima, H.; Nozawa, T.; Kawase, H.;
Abe, M.; Yokohata, T.; et al. MIROC-ESM 2010: Model description and basic results of CMIP5-20c3m
experiments. Geosci. Model. Dev. 2011, 4, 845–872. [CrossRef]
Miyazawa, Y.; Zhang, R.; Guo, X.; Tamura, H.; Ambe, D.; Lee, J.S.; Yoshinari, H.; Setou, T. Water mass
variability in the western north Pacific detected in a 15-year eddy resolving ocean reanalysis. J. Oceanogr.
2009, 65, 737–756. [CrossRef]
Miyazawa, Y.; Yamashita, N.; Taniyasu, S.; Yamazaki, E.; Guo, X.; Varlamov, M.S.; Miyama, T. Ocean
dispersion simulation of perfluoroalkyl substances in the Western North Pacific associated with the Great
East Japan Earthquake of 2011. J. Oceanogr. 2014, 70, 535–547. [CrossRef]
Ishizu, M.; Miyazawa, Y.; Tsunoda, T.; Guo, X. Development of a marine carbon model coupled with anoperational ocean model product for ocean acidification studies. Kaiyo Mon. 2018, 50, 217–225.
Miyazawa, Y.; Varlamov, M.S.; Miyama, T.; Guo, X.; Hihara, T.; Kiyomatsu, K.; Kachi, M.; Kurihara, Y.;
Murakami, H. Assimilation of high-resolution sea surface temperature data into an operational
nowcast/forecast system around Japan using a multi-scale three-dimensional variational scheme. Ocean Dyn.
2017, 67, 713–728. [CrossRef]
Heuven, V.S.; Pierrot, D.; Rae, J.W.B.; Lewis, E.; Wallance, D.W.R. Matlab program developed for CO2 system
calculation, ORNL/CDIAC-105b. In Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory;
U.S. Department of Energy: Oak Ridge, TN, USA, 2011. [CrossRef]

Sustainability 2019, 11, 2677

40.

41.
42.
43.

44.
45.
46.

47.
48.

49.
50.
51.
52.
53.
54.
55.

56.

57.

58.

59.

60.

27 of 28

Lewis, E.; Wallace, D.W.R. Program Developed for CO2 System Calculations. ORNL/CDIAC-105. In Carbon
Dioxide Information Analysis Center, Oak Ridge National Laboratory; U.S. Department of Energy: Oak Ridge, TN,
USA, 1998.
Orr, J.C.; Najjar, R.; Sabine, C.L.; Joos, F. Abiodic-HOWTO; Internal OCMIP Report; LSCE/CEA Saclay:
Gif-sur-Yvette, France, 1999; 25p.
Ewen, T.L.; Weaver, A.J.; Eby, M. Sensitivity of the inorganic ocean carbon cycle to future climate warming in
the UVic coupled model. Atmos. Ocean 2004, 42, 23–42. [CrossRef]
Smith, L.S.; Yamanaka, Y.; Pahlow, M.; Oschlies, A. Optimal uptake kinetics: Physiological acclimation
explains the pattern of nitrate uptake by phytoplankton in the ocean. Mar. Ecol. Prog. Ser. 2009, 384, 1–12.
[CrossRef]
Ishizu, M.; Richards, J.K. Relationship between oxygen, nitrate and phosphate in the world ocean based on
potential temperature. J. Geophys. Res. 2013, 118, 3586–3594. [CrossRef]
Redfield, A.C.; Ketchum, B.H.; Richards, A.F. The influence of organisusm on the Composition of Sea–Water.
In The Sea; Hill, M.N., Ed.; Interscience: New York, NY, USA, 1963; Volume 2, pp. 26–27.
Takatani, Y.; Enyo, K.; Iida, Y.; Kojima, A.; Sasano, D.; Kosugi, N.; Midorikawa, T.; Suzuki, T.; Ishii, M.
Relationships between total alkalinity in surface water and seas surface dynamic height in the Pacific Ocean.
J. Geophys. Res. Oceans 2014, 119, 2806–2814. [CrossRef]
Kantha, L.H. A general ecosystem model for applications to primary productivity and carbon cycle studies
in the global oceans. Ocean Model. 2004, 6, 285–334. [CrossRef]
Yamanaka, Y.; Yoshie, N.; Fujii, M.; Aita, N.M.; Kishi, M. An ecosystem model coupled withnitrogen-silicon-carbon cycles applied to Station A7 in the Northwestern Pacific. J. Oceanogr. 2004,
60, 227–241. [CrossRef]
Mussi, A. The solubility of calcite and argonate in seawater at various salinities, temperature and one
atmosphere total pressure. Am. J. Sci. 1983, 283, 780–799.
Zeebe, R.E.; Wolf-Gladrow, D. CO2 in Seawater: Equilibrium, Kinetics, Isotopes, 1st ed.; Elsevier: New York, NY,
USA, 2001; Volume 65, p. 360.
Mackenzie, F.T.; Lerman, A.; Anderson, A.J. Past and present of sediment and carbon biogeochemical cycling
models. Biogeosciences 2004, 1, 11–32. [CrossRef]
Brocker, W.S.; Peng, T.H. The role of CaCO3 compensation in the glacial to interglacial atmospheric CO2change. Glob. Biogeochem. Cycles 1987, 1, 5–29.
Wanninkhof, R. Relationship between wind speed and gas exchange over the ocean. J. Geophys. Res. Oceans
1992, 97, 7373–7382. [CrossRef]
Weiss, R. The Provisions of Social Relationships. In Doing unto Others; Rubin, Z., Ed.; Prentice Hall:Englewood Cliffs, UK, 1974; pp. 17–26.
Yasunaka, S.; Nojiri, Y.; Nakaoka, S.; Ono, T.; Whitney, F.; Telszewski, M. Mapping of sea surface nutrients in
the North Pacific: Basinwide distribution and seasonal to interannual variability. J. Geophys. Res. Oceans
2014, 119, 7756–7771. [CrossRef]
Goyet, C.; Healy, R.; Ryan, J. Global distribution of total inorganic carbon and total alkalinity below the
deepest winter mixed layer depths, ORNL/CDIAC-127, NDP-076. In Carbon Dioxide Information Analysis
Center, Oak Ridge National Laboratory; U.S. Department of Energy: Oak Ridge, TN, USA, 2000; pp. 1–40.
[CrossRef]
Key, M.R.; Kozyer, A.; Sabine, L.C.; Lee, K.; Wanninkhof, R.; Bullister, L.J.; Feely, A.R.; Millero, J.F.; Mordy, C.;
Peng, T.H. A global ocean carbon climatology: Results from global data analysis project (GLODAP). Glob.
Biogeochem. Cycles 2004, 18, GB4031. [CrossRef]
Yasunaka, S.; Nojiri, Y.; Nakaoka, S.; Ono, T.; Mukai, H.; Usui, N. Monthly maps of sea surface dissolved
inorganic carbon in the North Pacific: Basin-wide distribution and seasonal variation. J. Geophys. Res. Oceans
2013, 118, 3843–3850. [CrossRef]
Takahashi, T.; Sutherland, C.S.; Chipman, W.D.; Goddard, G.J.; Cheng, H.; Newberger, T.; Sweeney, C.;
Munro, R.D. Climatological distributions of pH, pCO2 , total CO2 , alkalinity, and CaCO3 saturation in the
global surface ocean, and temporal changes at selected locations. Mar. Chem. 2014, 164, 95–125. [CrossRef]
Yoshie, N.; Guo, X.; Fujii, N.; Komorita, T. Ecosystem and nutrient dynamics in the Seto Inland Sea, Japan,
Interdisciplinary Studies on environmental Chemistry—Marine Environmental modeling and analysis.
Model. Anal. Mar. Environ. Probl. 2011, 5, 39–49.

Sustainability 2019, 11, 2677

61.
62.
63.
64.
65.

66.
67.

28 of 28

Edward, F.K.; Thomas, K.M.; Klausmeier, A.C.; Litchman, E. Light and growth in marine phytoplankton:
Allometric, taxonomic, and environmental variation. Limnol. Oceanogr. 2015, 60, 540–552. [CrossRef]
Lalli, M.C.; Parsons, R.T. Biological Oceanography: An Introduction, 2nd ed.; Butterworth-Heinemann: Oxford,
UK, 1997; 320p.
Nelson, B.N.; Sigel, A.D. The global distribution and dynamics of chromophoric dissolved organic matter.
Annu. Rev. Mar. Sci. 2013, 5, 447–476. [CrossRef] [PubMed]
Jiang, L.Q.; Feely, A.R.; Carter, R.B.; Greeley, J.D.; Gledhill, K.D.; Arzayus, M.K. Climatological distributionof aragonite saturation state in the global oceans. Biogeochem. Cycles 2015, 29, 1656–1673. [CrossRef]
Yara, Y.; Vogt, M.; Fujii, M.; Yamano, H.; Hauri, C.; Steinacher, M.; Gruber, N.; Yamanaka, Y. Ocean
acidification limits temperature-induced poleward expansion of coral habitats around Japan. Biogeosciences
2012, 9, 4955–4968. [CrossRef]
Campbell, W.J. The lognormal distribution as a model for bio-optical variability in the sea. J. Geophys. Res.
1995, 100, 13237–13254. [CrossRef]
Tsuda, A.; Takeda, S.; Saito, H.; Nishioka, J.; Nojiri, Y.; Kudo, I.; Kiyosawa, H.; Imai, I.; Ono, T.; Shimamoto, A.;et al. A mesoscale iron enrichment in the western subarctic Pacific induces large centric diatom bloom.
Science 2003, 300, 958–961. [CrossRef] [PubMed]© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open accessarticle distributed under the terms and conditions of the Creative Commons Attribution(CC BY) license (http://creativecommons.org/licenses/by/4.0/).