Marine biogeochemical cycling and oceanic CO2 uptake simulated by the NUIST Earth System Model version 3 (NESM v3)

Abstract. In this study, we evaluate the performance of the Nanjing
University of Information Science and Technology (NUIST) Earth System Model
version 3 (hereafter NESM v3) in simulating the marine biogeochemical cycle
and carbon dioxide (CO2) uptake. Compared with observations, the NESM v3
reproduces the large-scale patterns of biogeochemical fields reasonably well
in the upper ocean, including nutrients, alkalinity, dissolved inorganic,
chlorophyll, and net primary production. Some discrepancies between model
simulations and observations are identified and the possible causes are
investigated. In the upper ocean, the simulated biases in biogeochemical
fields are mainly associated with shortcomings in the simulated ocean
circulation. Weak upwelling in the Indian Ocean suppresses the nutrient
entrainment to the upper ocean, thus reducing biological activities and
resulting in an underestimation of net primary production and the chlorophyll
concentration. In the Pacific and the Southern Ocean, nutrients are
overestimated as a result of strong iron limitation and excessive
vertical mixing. Alkalinity is also overestimated in high-latitude oceans
due to excessive convective mixing. The major discrepancy in biogeochemical
fields is that the model overestimates nutrients, alkalinity, and dissolved
inorganic carbon in the deep North Pacific, which is caused by the excessive
deep ocean remineralization. The model reasonably reproduces present-day
oceanic CO2 uptake. Model-simulated cumulative oceanic CO2 uptake
is 149 PgC between 1850 and 2016, which compares well with data-based
estimates of 150±20 PgC. In the 1 % yr−1 CO2 increase
(1ptCO2) experiment, the diagnosed carbon-climate
(γ=-7.9 PgC K−1) and carbon-concentration sensitivity parameters (β=0.88 PgC ppm−1) in the NESM v3 are comparable with those in Coupled Model
Intercomparison Project phase 5 (CMIP5) models (β: 0.69 to
0.91 PgC ppm−1; γ: −2.4 to −12.1 PgC K−1). The nonlinear interaction between
carbon-concentration and carbon-climate sensitivity in the NESM v3 accounts
for 10.3 % of the total carbon uptake, which is within the range of CMIP5
model results (3.6 %–10.6 %). Overall, the NESM v3 can be
employed as a useful modeling tool to investigate large-scale interactions
between the ocean carbon cycle and climate change.



3
Given the importance of carbon cycle feedback in current and future global climate change, it is necessary to include the representation of the global carbon cycle in climate system models (Denman et al., 2007).
The first and second generations of the NUIST climate system model have good skill in simulating internal climate modes, global monsoon, and projecting future climate change (Li et al., 2018;Yang and Wang, 2018a;Yang et al., 2018b). However, the previous generations of NESM do not include an active 5 ocean biogeochemical model and cannot be used to study the ocean carbon cycle (Cao et al., 2015).
Recently, the third version of the NUIST earth system model (NESM v3) was developed as a newly registered model to CMIP phase 6 (CMIP6, Cao et al., 2018), which couples the Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES v2) to represent the ocean biogeochemical processes (Aumont et al., 2015). 10 The objective of this manuscript is to evaluate the performance of NESM v3 in simulating marine carbonrelated biogeochemical tracers and oceanic anthropogenic CO2 uptake. As a newly developed earth system model and a new member of the CMIP community, it is essential to evaluate the model's ability against observational data. First, we analyze the present-day mean state of macro-nutrients, chlorophyll, net primary production (NPP), sea-air CO2 flux, dissolved inorganic carbon (DIC), and alkalinity against 15 available observations and observational-based estimates. Then, we evaluate the amount, spatial distribution, and sensitivity parameters of the anthropogenic carbon uptake. The amount and spatial distribution of oceanic anthropogenic carbon uptake during the historical period and future scenarios are compared with observations and CMIP5 model results. The sensitivity parameters and nonlinearity of oceanic CO2 uptake diagnosed from the NESM v3 are compared with those from CMIP5 models. We also 20 provide a supplementary material that includes comparisons of biogeochemical fields in the NESM v3 and IPSL-CM5A-LR (hereafter IPSL) that shares the same marine biogeochemical component (PISCES) as used in NESM v3.
The outline of the paper is the following. In Section 2, we describe the NESM v3 with a focus on the ocean carbon cycle component, as well as the setup of model simulations. The results of model 25 simulations are analyzed in Section 3. Conclusions and discussions are presented in Section 4. Detailed descriptions of the physical components, major improvements, and model tuning procedures of the NESM v3 are documented in Cao et al. (2018). Here we give a brief introduction.
In this study, we use the low-resolution version of the NESM v3. The atmospheric resolution is T31L31 which has a horizontal resolution of ~ 3.75° latitude by 3.75° longitude and 31 layers. The atmospheric 10 model and land surface model are originally adopted from ECHAM v6.3. The detailed information is shown in Stevens et al. (2012) and Giorgetta et al. (2013). The sea-ice model includes four ice layers and one snow layer with a multi-layer thermodynamic scheme (Hunke et al., 2010;Cao et al., 2018). Ocean model has the ORCA2 global ocean configuration, which is a type of tripolar grid. It is based on a 2 degree Mercator mesh and has 31 layers with the thickness of the ocean layer increasing from 10 m in the 15 upper ocean to 500 m at 5000 m depth. A local transformation is applied in the tropics to refine the resolution to up to 0.5 degree at the equator. In the ocean model, the incoming solar radiation can penetrate to the upper ocean layers as deep as 391 m, and a bio-model penetration parameterization scheme is used to calculate the distribution of solar radiation (Lengaigne et al., 2009). The ocean background vertical diffusivity is modified in the NESM v3, where the constant value is replaced by latitude-dependent values 20 (Jochum et al., 2009). The parameterization scheme of the vertical diffusivity is detailed in the supplementary, and the global distribution of vertical diffusivity is also shown (Fig. S1). Compared to the original vertical diffusivity coefficient constant of 0.12 cm 2 /s, the coefficients of the tropical ocean are reduced and that of the middle and high latitude oceans are increased, especially in the middle latitude oceans (24°N~33°N and 24°S~33°S). Also, the NESM v3 incorporates the parameterization of brine 25 rejection in the ocean model and the reference sea-ice salinity is set as 4 PSU as suggested by Smith et al (2010).  (Aumont et al., 2002). It can be used for both regional and global simulations of lower trophic levels of the marine ecosystem and ocean carbon cycle (Bopp et al., 2005;Resplandy et al., 2012;Séférian et al., 5 2013). In the current version, there are 24 prognostic tracers in total, including dissolved inorganic and organic carbon, alkalinity, chlorophyll, and nutrients. We use the same biogeochemical parameters as that used in Aumont et al. (2015). The only exception is the advection scheme for passive tracers. Here we use the Total Variance Dissipation (TVD) formulation instead of the Monotone Upstream Scheme for Conservative Laws (MUSCL) formulation to keep the advection scheme to be consistent with the one 10 used in the physical ocean model. Both TVD and MUSCL schemes have good performance in biogeochemical modeling. The MUSCL scheme has better performance in resolving the small scales, while TVD scheme minimizes systematic error through numerical diffusion, and is a better option for coarse-resolution models (Lévy et al., 2001a).
Two different types of phytoplankton: nanophytoplankton and diatoms, and two size classes of 15 zooplankton: mesozooplankton and microzooplankton are presented in the model. The life cycle of phytoplankton is regulated by processes of growth, mortality, aggregation, and grazing by zooplankton (Aumont et al., 2015). The growth rate of phytoplankton is determined by temperature, photosynthetically active radiation (PAR), and availability of nutrients, including phosphate, nitrate, silicate, iron, and ammonium. The mortality rate of phytoplankton is set as a constant and is identical for nanophytoplankton 20 and diatoms. The aggregations of nanophytoplankton, which transform the dissolved organic carbon (DOC) to the particular organic matter (POM), only depend on the shear rate, as the main driver of aggregation is the local turbulence. In the NESM v3, this shear rate is set to 1 s -1 in the mixed layer and 0.01 s -1 below. The same is assumed for diatoms, while the aggregations of diatoms are further enhanced by the nutrients co-limitation. For all species, the phosphate, nitrate, and carbon are linked by a constant 25 Redfield ratio. In NESM v3, the Redfield ratio of C: N: P is set to be 122:16:1 (Takahashi et al., 1985) and the -O/C ratio is set to 1.34 (Kortzinger et al., 2001). In contrast, the Fe / C, chlorophyll / C, and silicon / C ratio are prognostically simulated by the model based on the external concentrations of the 6 limiting nutrients as in the quota-approach (McCarthy, 1980;Droop, 1983;Aumont et al., 2015).
The remineralization of semi-labile DOC can occur in either oxic water or anoxic water depending on the local oxygen concentration, and their degradation rates are specified and identical for oxic respiration and denitrification. Detritus is represented by different types, including POM, calcite, iron particles, and biogenic silica. The POM is partitioned into two size classes: a smaller class (POC: 1-100μm) and a larger 5 class (GOC: 100-500μm). The sinking speed of GOC (50-200 m d -1 ) increases with depth and is much faster than that of POC (3 m d -1 ). A fraction of phytoplankton would be turned to the POM through the processes of mortality and aggregation. The fate of mortality and aggregation of nanophytoplankton depends on the proportion of the calcifying organisms. For nanophytoplankton, it is assumed that half of the calcifying organisms are associated with the calcifying organisms. Because the density of the calcite 10 is larger than that of organic matter, 50% of the dying calcifiers are routed to the fast-sinking particles.
The same is assumed for the mortality of diatoms, and 50% of the dying diatoms are turned to the POM due to the larger density of biogenic silica compared to that of organic matter. The degradation rate of the POM depends on the local temperature with a Q10 factor (temperature dependence ratio) of 1.9.
The geochemical boundary condition accounts for the external nutrient supply from five different sources, 15 including atmospheric dust deposition of iron and silicon, river recharge of nutrients, dissolved carbon, and alkalinity, atmospheric deposition of nitrogen, and sediment mobilization of sedimentary iron. In the NESM v3, atmospheric deposition and river recharge are prescribed and sediment mobilization is parameterized. At the bottom of the ocean, different sediment parameterization schemes are applied to biogenic silica, POM, and particulate iron. The amount of permanently buried biogenic silica is assumed 20 to balance the external source, and the burial efficiency of POM is determined by the organic carbon sinking rate at the bottom follows the algorithm proposed by Dunne et al. (2007). All the particulate iron would be buried into the sediment once they reach the ocean bottom. The amount of the unburied calcite and biogenic silica would dissolve back into the ocean water instantaneously. In this study, the initial conditions of the biogeochemical model have been adopted from World Ocean Atlas 2009 (WOA09, 25 Garcia, et al., 2010) and Global Ocean Data Analysis Project (GLODAP, Key et al., 2004;Sabine et al., 2004) dataset.
Carbonate chemistry is formulated based on the Ocean Carbon-Cycle Model Intercomparison Project 7 (OCMIP-2) protocol (more information can be accessed at http://ocmip5.ipsl.jussieu.fr/OCMIP/). The quadratic wind-speed formulation proposed by Wanninkhof (1992) is used to compute the air-sea exchange of carbon and oxygen.

Simulations
In total, there are eight different simulations in this study, including one fully coupled spin-up simulation 5 for 2000 years, one PI-control run (CTRL) for 251 years, three transient runs driven by forcing conditions from historical observational data and shared socio-economic pathway scenarios (SSP5-8.5) from 1850 to 2100, i.e., fully-coupled (FC), biogeochemically-coupled (BC), and radiatively-coupled (RC), and three idealized 1ptCO2 runs for 140 years, i.e., FC-1%, BC-1%, and RC-1%. The following is the detailed experimental design. Following the protocol of CMIP6 historical and SSP5-8.5 experiment design (Eyring et al., 2016;Jones et al., 2016), the model is further integrated with changing conditions, including GHGs, ozone, aerosol, Also, we performed a similar set of simulations to have a direct comparison with CMIP5 results, i.e., an idealized 1ptCO2 run, in which atmospheric CO2 concentration increases at a rate of 1% per year starting from the end state of the pre-industrial simulation with other GHGs concentration remaining at the preindustrial level. The simulation lasted for 140 years until atmospheric CO2 concentration has quadrupled. 5 Following Friedlingstein et al., 2006 andArora et al., 2013, we performed three types of experiments (biogeochemically coupled, radiatively coupled, and fully coupled) for each scenario to separate the effect of atmospheric CO2 and global warming on the ocean carbon cycle: 1) BC simulations in which the code of the ocean carbon cycle sees changing atmospheric CO2, but the code of atmospheric radiation sees a constant pre-industrial concentration of CO2. In this way, the ocean 10 carbon cycle is only affected by changing atmospheric CO2, but no direct effect of CO2-induced warming; 2) RC simulations in which the code of the ocean carbon cycle sees pre-industrial atmospheric CO2, but the code of atmospheric radiation sees the changing concentrations of atmospheric CO2. In this way, the ocean carbon cycle is only affected by CO2-induced warming, but no direct effect of changing atmospheric CO2. 3) FC simulations in which both the codes of the ocean carbon cycle and atmospheric radiation see the changing concentrations of atmospheric CO2. In this way, the ocean carbon cycle is affected by changes in both atmospheric CO2 and CO2-induced warming.

Evaluation data
Global ocean distribution data of nutrients concentrations, including nitrate, phosphate, and silicate, are 20 from the World Ocean Atlas 2018 (WOA18, Garcia, et al., 2018). In this study, we assume the density of the seawater in the model and observations are the same, and then the unit of observations is converted from umol/kg to mmol/m3 by multiplying the modeled density (in the unit of kg/m3). Geographic distributions of DIC, alkalinity, and anthropogenic carbon is taken from the GLODAP v2 (Key et al., 2015;Lauvset et al., 2016). Both WOA18 and GLODAP v2 data have a horizontal resolution of 1°×1° 25 with 33 levels and represent the climatology in recent decades. We compared modeled chlorophyll in recent decades with the SeaWiFS dataset (NASA Goddard Space Flight Center, 2014), GlobColour 9 merged data (Maritorena, et al., 2010), and Ocean Colour Climate Change Initiative (OCCCI) merged data (http://www.oceancolour.org/). OCCCI and GlobColour incorporate the same datasets, while their uncertainty information and algorithms are not the same.
Mode-simulated NPP is compared with moderate Resolution Imaging Spectroradiometer (MODIS) estimated marine NPP based on three different algorithms, including the Standard Vertically Generalized 5 Production Model (VGPM), Eppley-VGPM, and the carbon-based Production Model (CbPM). The datasets can be accessed at http://www.science.oregonstate.edu/ocean.productivity/index.php. In the VGPM and Eppley-VGPM, NPP is described as the product of chlorophyll and photosynthetic efficiencies Falkowski, 1997a, 1997b), while the Eppley-VGPM emphasizes the photoacclimation effect at high SSTs, i.e., the growth rate is higher at high-temperature regions (Eppley, 10 1972;Morel, 1991). In the CbPM, NPP is described as the product of carbon biomass and growth rate (Behrenfeld et al. 2005;Westberry et al. 2008). All three datasets have a horizontal resolution of 1/12°×1/12° from 2003 to 2014. The distribution of observed surface ocean sea-air CO2 flux for the reference year of 2000 is taken from Takahashi et al (2009) and has a spatial resolution of 4° latitude by 5° longitude. 15 To have a direct comparison between the NESM v3 output and observations, we interpolated all modeled results and observations to a 1°×1° grid using the distance-weighted average remapping method, except for the sea-air CO2 flux. Due to the low resolution of observational sea-air flux, we interpolated the modeled sea-air CO2 flux to a 4° × 5° grid.

Nutrient decomposition
To figure out whether the ocean circulation or biogeochemical processes causes the bias of nutrients in the model simulation, we decomposing phosphate to its preformed and regenerated components (Weiss et al., 1970;Duteil et al., 2012). The regenerated phosphate is released through the remineralization processes of organic matter, and the preformed phosphate is the remaining biotically unutilized surface 25 phosphate, which is transported into the ocean interior by ocean circulation. The regenerated and preformed phosphate are computed as: Where AOU is the apparent oxygen utilization, which represents the biological consumption of oxygen.
It is computed as the difference between oxygen saturation and simulated oxygen concentration. :− 2 represents the oxidation ratio of phosphate and oxygen, which is set to 1/163 in the NESM v3. P represents 5 the simulated phosphate concentration.

Nutrient limitation
In the NESM v3, the nutrients limitation coefficient (0~1) is computed from the Michaelis-Menten equation as follow: Where MM is the Michaelis-Menten coefficient, N is the nutrient concentration, and K is the halfsaturation coefficient, which is parameterized based on half-saturation constant and concentrations of nutrients, phytoplankton, and diatoms (Aumont et al., 2015).
We calculated the annual mean nutrient limitation coefficient of each nutrient (phosphate, nitrate, silicate, and iron) and then considered the nutrient with the lowest limitation coefficient as the most limiting factor. 15 Temperature and light are assumed to be the most limiting factor when the annual mean nutrient coefficients are greater than 0.9 for all nutrients (Moore et al., 2013).

Carbon-concentration and carbon-climate sensitivity parameters
Following Arora et al. (2013), we diagnose the carbon-climate and carbon-concentration sensitivity parameters from two types of experiments performed by a subset of CMIP5 models, i.e., BC simulations 20 and RC simulations.
In the biogeochemically-coupled simulations where the ocean carbon uptake is only affected by changing atmospheric CO2, the relationship between atmospheric CO2 concentration and sea-air CO2 flux can be simplified as:

25
Where ′ represents oceanic carbon uptake change in the biogeochemically coupled simulation. In the 11 radiatively-coupled simulations where the oceanic carbon uptake is only affected by temperature change, the relationship between temperature and sea-air CO2 flux can be simplified as: Where ′ represents oceanic carbon uptake change in the radiatively coupled simulation.

Nutrients
Nutrients play vital roles in the ocean biogeochemical cycle. A lack of nutrients would limit the growth of phytoplankton. Figure  water to the surface (Whitney, 2011). Also, the strong iron limitation (Fig. 4) that reduces the biological uptake of the macro-nutrients is one of the main causes of the high nutrients level in the Southern Ocean.
The relatively high concentration of nutrients is simulated in the subarctic Pacific Ocean, and the mideastern Pacific Ocean, which is about 50% of the concentration in the Southern Ocean. Relatively low 25 concentration of nutrients, less than 20% of the concentration in the Southern Ocean, is simulated in 12 subtropical regions. Some noticeable discrepancies between model simulations and observations are also found ( Fig. a3, b3, c3). Phosphate and nitrate are overestimated in the Southern Ocean and the Pacific Ocean but are underestimated in the Indian Ocean, Subarctic Pacific, and Middle-low latitude Atlantic.
Except for the Indian Ocean and Subarctic Pacific, silicate is overestimated over all of the global ocean, resulting in the high global mean concentration (50% higher than the WOA18 observations). To further analyze the possible reasons for discrepancies in nutrients distributions, we decompose phosphate to its preformed and regenerated components (Weiss et al., 1970;Duteil et al., 2012) and 15 compare the results with the WOA18 observations (Fig. 3). For the global, Atlantic, and Pacific Ocean, the preformed phosphate diagnosed from the model accounts for 51%, 47%, and 57% of the total phosphate inventory, and the result diagnosed from the WOA18 is 57%, 55%, and 64%, respectively. A The NESM v3 simulates a dipole pattern of preformed phosphate bias above 1500 m depth, which is probably associated with the overestimated vertical mixing that brings too much nutrient-rich water from 25 the intermediate ocean to the surface (Fig. a2, a4, and a6). The noticeable positive biases of regenerated phosphate are found in the deep Northern Pacific. In the model simulation, the high regenerated phosphate water in the North Pacific is too deep, and the biases resemble the difference found in latitudinal-depth 13 distributions of nutrients. In the deep ocean, preformed phosphate is only affected by ocean circulation, while regenerated phosphate is affected by both circulation and remineralization. In the deep ocean, the NESM v3 simulates the preformed phosphate well but overestimates the regenerated phosphate, suggesting that the overestimated nutrients in the North Pacific deep ocean are mainly caused by biological processes. 5 Next, we present the NESM v3-simulated patterns of nutrient limitation. As shown in figure 4, the limiting patterns of nanophytoplankton and diatoms are similar in the mid-low latitude oceans. Iron is the most limiting nutrient for both nanophytoplankton and diatoms in the equatorial Pacific Ocean and the Southern Ocean. Nitrate is the most limiting factor in the subtropical Pacific Ocean, and phosphate is the most limiting factor in the Indian Ocean and middle-low latitude of the Atlantic Ocean. At high latitude oceans, 10 nanophytoplankton is mostly limited by the available light and temperature, while diatoms are mostly limited by silicate. The NESM v3 simulated limiting pattern is generally consistent with the results diagnosed from IPSL-CM4A-LOOP (Schneider et al., 2008), except that the iron limitation diagnosed from the NESM v3 is stronger in the Pacific and the Southern Ocean.

Biological Production
15 Figure 5 shows the modeled spatial distribution of annual mean surface chlorophyll concentration from 1998 to 2014 compared with SeaWiFS observational data, GlobColour merged data, and OCCCI merged data.
In the NESM v3, chlorophyll in both nanophytoplankton and diatoms are parameterized based on the photo-adaptive model (Geider et al., 1997) in which chlorophyll is regulated by the chlorophyll-to-carbon Ocean, maritime continent, and the tropical Atlantic Ocean. This underestimation is partly associated with the deficiencies in modeled coastal dynamics, which is usually not represented well by the relatively coarse global ocean models (Aumont et al., 2015). It is reported that the observed chlorophyll distribution is better reproduced when PISCES is coupled to a higher resolution ocean circulation model (Lee et al., 5 2000;Hood et al., 2003;Kone et al., 2009). Also, we inspect an underestimation of chlorophyll over the Northern Indian Ocean. This is associated with the underestimation of nutrients over the Indian Ocean ( Fig. 1) that increases nutrients limitation and inhibits growth.
In the Southern Ocean where the seawater is typically characterized by high nutrients and low chlorophyll (Lin et al., 2016), noticeable discrepancies are seen among different observational datasets that are 10 associated with different algorithms used for these products. For example, in the Southern Ocean, chlorophyll derived from reflectance by standard algorithms tend to be underestimated by a factor of about 2 to 2.5 (Kahru and Mitchell, 2010). In the Southern Ocean, the NESM v3 overestimates the chlorophyll concentration to the east of 150°E and underestimates it near the International Date Line (180°E). In the Atlantic part of the Southern Ocean, the modeled chlorophyll concentration is within the 15 range of observational estimates, higher than the SeaWiFS but lower than the GlobColour and OCCCI. on the MODIS observations. Similar to the CbPM, the NESM v3 also simulates the NPP as the product of phytoplankton biomass and growth rate. However, the calculation of growth rate in the NESM v3 is 20 more complex than that in CbPM, which involves chlorophyll, nutrients availability, temperature, respiration, and PAR.
There are significant differences between the NESM v3 simulated NPP and VGPM because the formulation of growth rate in the NESM v3 follows Eppley (1972), i.e., the growth rate is higher at hightemperature regions. Therefore, the NESM v3 estimates more NPP in low latitude oceans and less at high 25 latitude oceans than that in the VGPM (Fig. 6e). The climatology of the NESM v3 simulated vertically integrated NPP resembles Eppley-VGPM and CbPM estimates (Fig. 6a, c, d). The main observed spatial pattern of high concentrations of NPP in the eastern equatorial Pacific and middle-latitude oceans and low 15 concentration of NPP in the middle-low latitude oceans, the Southern Ocean, and high latitude oceans are reproduced by the NESM v3. Also, the NESM v3 reproduced the high concentrations of NPP in lowlatitude coastal regions to some extent. Comparison between three observational data-based estimates and the NESM v3 also suggests that temperature-dependence is important to produce the meridional distribution pattern of marine NPP with the high level of NPP over the high-latitude oceans and low level  (Longhurst et al., 1995;Antoine et al., 1996;Behrenfeld and Falkowski, 1997b;Behrenfeld et 15 al., 2005). Global NPP simulated by CMIP5 models also shows a wide range of values from 30.9 to 78.7 PgC yr -1 . NESM v3 simulated global NPP is within the range of data-based estimates and current CMIP5 model estimates. Of the NESM v3 simulated global ocean NPP, 20% is contributed by diatoms, and 80% is contributed by nanophytoplankton. For comparison, from the data-based estimate, 7% to 32% of the total NPP is associated with diatoms (Uitz et al.,2010;Hirata et al., 2011), while ocean 20 biogeochemical models estimate that 15% to 30% global NPP is from diatoms (Aumont et al., 2003;Dutkiewicz et al., 2005;Yool et al., 2011).  The model's skill in simulating alkalinity is moderate (PCC = 0.56). The observed global spatial pattern of alkalinity is generally simulated by the NESM v3. The model reproduced the observed high alkalinity in the subtropical surface oceans and low alkalinity near the maritime continent, and the modeled global upper ocean mean alkalinity only has a minor negative bias of 0.45%. The major discrepancies are seen in the Southern Ocean and the subarctic Pacific with a positive bias of more than 80 mmol/m 3 . In highlatitude oceans, convective mixing of alkalinity-rich deep water is an important factor of changing upper ocean alkalinity, and SST can be used as a proxy of the convective mixing change (Lee et al., 2006). An 5 underestimation of SST of 1 o C is simulated at high latitude oceans (figures not shown), indicating a stronger convective mixing, which may explain the overestimated alkalinity at high latitude oceans. The alkalinity has a negative bias of more than 60 mmol/m 3 near the maritime continent, where the alkalinity concentration is usually related to salinity (Lee et al., 2006). Cao et al. (2018) found that the NESM v3 simulates excessive precipitation over the maritime continent, which results in an underestimation of the  3.4 Assessment of biogeochemical fields by Taylor diagram Figure 9 respectively compares the spatial patterns of the NESM v3 and IPSL simulated biogeochemistryrelated fields with corresponding observations using a Taylor diagram (Taylor, 2001). In summary, modelsimulated statistical patterns of the upper ocean nutrients compare well with observations, while the simulated spatial patterns of chlorophyll, primary production, and alkalinity show larger discrepancies 5 from observations. It is noted that chlorophyll and NPP are not directly observed but diagnosed from the observation-based data, and thus their estimations are subject to considerable uncertainties. Compared with biogeochemical fields in IPSL, the NESM v3 has comparable skill in reproducing spatial distributions of nutrients and chlorophyll, but less skill in reproducing DIC and alkalinity with relatively larger SDs. 10 We also examine other CMIP5 models results that documented in previous studies (Moore et al., 2013;Anav et al., 2013;Tjiputra et al., 2013;Séférian et al., 2013). The skill of the models are different according to the biogeochemical fields examined. For example, PCCs between nutrients and observations in Community Earth System Model (CESM) are about 0.8, which is lower than that in the NESM v3, but CESM has a better representation of chlorophyll distribution (Moore et al., 2013). A further study is 15 needed to rank those models, which is beyond the scope of this study. Nevertheless, the NESM v3 shows comparable skill in simulating upper ocean biogeochemical fields in the present-day with other CMIP5 models.

Oceanic anthropogenic CO2 uptake during the historical period
In this section, we compare the NESM v3 simulated anthropogenic carbon uptake during the historical 20 period (FC) against available observational data-based estimates.
First, we compare the NESM v3 simulated sea-air CO2 flux against available observations for the reference year of 2000 (Takahashi et al., 2009). As shown in Fig. 10, the NESM v3 realistically reproduces the large-scale pattern of observed sea-air CO2 flux with CO2 outgassing in the equatorial oceans and uptake in the mid-to-high latitude oceans (PCC=0.71 and SDs =1.04). For both observation and model 25 results, strong CO2 uptake is found in the North Atlantic where sea surface temperature is low and the formation of deep water is active. Compared to the data-based estimates, modeled sea-air CO2 flux is 18 overestimated in the tropical Pacific, the Southern ocean, and the Northern Pacific (near 30°N), while the strongest underestimates of modeled sea-air CO2 flux are seen in the high-latitude oceans ( Fig. 9c and   9d). The globally integrated ocean uptake flux from observation is 1.6 ± 0.9 PgC per year in the year 2000 (Takahashi et al., 2009), while the value is 2.8 PgC per year from the model simulation. The difference between model and observation is mainly originated from positive bias (~1 PgC yr -1 ) in the pre-industrial 5 steady-state oceanic CO2 uptake due to the 3-dimensional correction of nutrient and alkalinity in the PISCES model (Séférian et al., 2013). The 3-dimensional correction refers to that the global inventory of nutrient and alkalinity are restored toward the observations every January 1 (Aumont et al., 2015). In the NESM v3, the pre-industrial steady-state of total oceanic CO2 uptake is 1.  . 11a) and from 1992 to 1996 (Fig. 11c) are compared with GLODAP v2 (Fig. 11b) and GLODAP v1 (Fig. 11d), respectively. NESM v3 reasonably captures carbon, while it is 17.6% in the observation ( Fig. 11a and 11b). In other oceans, the large inventory is mainly found in the middle-latitude areas near 30°N and 30°S. In the Southern Hemisphere Oceans, 58.9% of the global oceanic anthropogenic DIC inventory is simulated, compared to the value of 62.6% in the observation.
The most noticeable discrepancy between the GLODAP v2 and model simulation around 2002 is found 5 in the south of 50°S. Only 8.3% of the global oceanic anthropogenic DIC inventory is simulated, while the value is 15.5% in the observation. However, we noticed that the vertically integrated anthropogenic DIC concentration is also low in the southern ocean south of 50°S in the GLODAP v1 and only 9.9% of the global inventory is stored in this region. It is noted that anthropogenic DIC in the GLODAP is diagnosed by a crude application of the transit time distribution method, and thus the results are subject 10 to considerable uncertainties (Lauvset et al., 2016). and high latitude oceans where vertical mixing is strong . Similar to the vertically integrated inventory of anthropogenic DIC (Fig. 11), the major discrepancy of anthropogenic DIC in the latitudinal-depth distribution is also found in the Southern Atlantic south of 50°S.
3.6 Sensitivity of the oceanic CO2 uptake to increasing atmospheric CO2 and global warming The ocean carbon cycle is regulated by changes in atmospheric CO2 and physical climate (Doney et al., 25 2004). Increasing atmospheric CO2 affects oceanic CO2 uptake directly. Meanwhile, global warming also affects the ocean carbon cycle via changes in climatic fields (Gregory et al., 2005;Pierce et al., 2012). In 20 this section, we first presented the NESM v3 simulated physical climate change and oceanic CO2 uptake under the historical and SSP5-8.5 scenario. Then, we present NESM v3 simulated carbon cycle sensitivity parameters and their nonlinearity in the SSP5-8.5 and 1ptCO2 runs.
3.6.1 NESM v3 simulated physical climate change under historical and SSP5-8.5 scenario Figure 13 shows the NESM v3 simulated changes (relative to control simulation) in global annual mean  (Collins and Knutti, 2013;Knutti and Sedláček, 2013). It is noted that the CMIP6 input forcing is used in this study and the atmospheric CO2 concentration at the end of the 21 st century in SSP5-8.5 is about 10% higher than the concentration in the RCP 8.5 scenario.
In the year 2100, SAT change is 6.3 (6.7 and 0.8) K in the FC (RC and BC) simulation. With the increasing atmospheric temperature, the global ocean also becomes warmer in FC and RC simulations, 15 reducing CO2 solubility and oceanic CO2 uptake.
Modeled MLD decreases since the 1980s. In the year 2100, the MLD change is -8.5 (-8 and -1) meter in the FC (RC and BC) simulation. The reduction of mixed layer depth indicates a more stratified upper ocean. A substantial weakening of AMOC intensity in the FC simulation is projected for the 21 st century, which is associated with ocean surface warming and increased freshwater input into the North Atlantic 20 (Gregory et al., 2005). In the pre-industrial period, the model-simulated AMOC index at 30°N is 17.5 Sv (1Sv =10 6 m 3 s -1 ), within the range of 14 to 31 Sv from CMIP5 models (Weaver et al., 2012 Ocean Circulation and Heatflux Array) is 17.5 ± 3.8 Sv (Rhein et al., 2013). By the year 2100, the 25 simulated intensity of AMOC declines to 8.0 Sv. In the FC simulation, the simulated 54% weakening of AMOC by the end of this century is at the higher end of what is simulated by CMIP5 models that range 21 from 15% to 60% under the RCP 8.5 scenario (Cheng et al., 2013). The higher atmospheric CO2 concentration at the end of 2100 in the SSP5-8.5 may partly explain the larger AMOC change in this study.
Also, Cao et al. (2018) pointed out that the equilibrium climate sensitivity to CO2 forcing in the NESM v3 is about 10% higher than the CMIP5 ensemble.
3.6.2 NESM v3 simulated oceanic CO2 uptakes under historical and SSP5-8.5 scenario 5 Figure 14 shows the time evolution of the oceanic CO2 uptake from the BC, RC, FC, and the linear sum of BC and RC. In the BC simulation, the global ocean absorbed a total amount of 662 PgC of anthropogenic CO2 from the atmosphere by the year 2100. In the RC simulation, the increased sea surface temperature, enhanced ocean stratification, and the weakened AMOC all act to decrease CO2 uptake (Cox et al., 2000;Zickfeld et al., 2008;Roy et al., 2011;Goris et al. 2015). As a result, global warming alone 10 causes the ocean to release CO2 into the atmosphere by reducing CO2 solubility, enhancing the oceanic pCO2 and altering biological rates (Steinacher et al., 2010;Olonscheck et al., 2013;Lewandowska et al., 2014;Cao et al. 2017). By the year 2100, the modeled cumulative CO2 uptake is -35.9 PgC. In the FC simulation, oceanic CO2 uptake is affected by both the increase in atmospheric CO2 and global warming.
By the end of the 21st century, simulated cumulative oceanic CO2 uptake since the pre-industrial era is 15 567 PgC, which is within the ranges from 420 PgC to 600 Pg C from CMIP5 models results under the RCP 8.5 scenario . The sum of the simulated oceanic CO2 uptake from the BC and RC simulations (626 PgC) is larger than that from the FC run (567 PgC), indicating that the effect of increasing atmospheric CO2 (carbon-concentration sensitivity) and the effect of global warming (carbonclimate sensitivity) on the oceanic CO2 uptake are not exactly additive. This nonlinearity was also found 20 in previous studies (Boer and Arora, 2009;Gregory et al., 2009;Schwinger et al., 2014). The NESM v3 simulated nonlinearity (i.e., BC+RC-FC) is 59 PgC by the end of the 21st century, which is larger than the absolute value of the radiative effect on oceanic carbon uptake (-35.9 PgC).
To better understand oceanic CO2 uptake in response to changing atmospheric CO2 and global warming in model simulation, Figure 15 shows the spatial distribution of anthropogenic sea-air CO2 flux at the end 25 of the 21 st century (averaged over the year 2091 to 2100) under the SSP5-8.5 scenario from FC, RC, and BC simulations, and the difference between FC simulation and the sum of RC and BC simulations.

22
In the BC simulation, the oceanic anthropogenic CO2 uptake is 8.0 Pg C per year at the end of the 21 st century. The ocean absorbs atmospheric CO2 in most regions except for a few scattered grid points at the mid-latitudes with slight CO2 outgassing. The strongest CO2 uptake of about 150 g C m -2 yr -1 is found in the North Atlantic, subarctic Pacific, and the Southern Ocean. Results from the RC simulation show CO2 outgassing in large parts of the global ocean as a result of global warming. In the Arctic Ocean, warming 5 induces a net uptake of CO2 of 0.07 PgC yr -1 because the reduced sea-ice extent under global warming allows more open seawater to absorb atmospheric CO2. In the Northern Atlantic, the capacity of the ocean uptakes CO2 is significantly suppressed due to the reduced AMOC.
The FC simulation shows the combined effects of the increasing atmospheric CO2 and global warming (Fig. 14c). Oceanic CO2 uptake is simulated in most regions with the strongest uptake in the Southern 10 Ocean, indicating the dominant role of the increasing atmospheric CO2. CO2 outgassing is seen in the subtropical Pacific, indicating that the radiative effect dominates the response of oceanic CO2 uptake in this region. The nonlinearity of oceanic carbon uptake sensitivity during the 2090s is shown in Figure   15d. In the NESM v3, a relatively large nonlinearity is simulated in the Northern Atlantic north of 45°N (19.8% of the total nonlinearity) and the Southern Ocean south of 40°S (35.3% of the total nonlinearity), 15 which is consistent with the findings of previous studies (Zickfeld et al., 2011;Schwinger et al., 2014).
The background simulation effects can partly explain the nonlinearity. Compared with the RC simulation, more carbon is subject to the impact of climate change in the FC simulation. As a consequence, in the FC simulation, the increased temperature would have a larger effect on CO2 solubility and buffer factor (Yi et al., 2001). Also, reduced ocean circulation and increased ocean stratification would slow down the 20 transport of anthropogenic CO2 from the surface to the deep ocean. Thus, compared to the BC simulation, slowing ocean ventilation in the FC simulation would cause a larger reduction in oceanic CO2 uptake.
The oceanic carbon uptake in the FC simulation is lower than the sum of the BC and RC simulations, which is consistent with other CMIP5 models (Schwinger et al., 2014). The above results also indicate that oceanic CO2 uptake in high-latitude oceans is more sensitive to both the increasing atmospheric CO2 25 concentration and global warming than low-latitude oceans, as well as their nonlinear interactions.
3.6.3 carbon-concentration and carbon-climate sensitivity parameters diagnosed from the SSP5-8.5.

23
In this section, we investigated oceanic CO2 uptake under the framework of the carbon-concentration and carbon-climate sensitivity parameters. Figure 16 shows the change in ocean carbon storage against the change in the atmospheric CO2 concentration (Fig. 16a) and the global annual mean surface temperature (Fig. 16b), respectively. The derived evolution of the carbon-concentration sensitivity parameter as a function of atmospheric CO2 5 concentration and carbon-climate sensitivity parameter as a function of the change in temperature is shown in Fig. 16c and 16d, respectively.
As shown in Fig. 16, in the BC and RC simulations, modeled ocean storage of anthropogenic CO2 scales roughly linearly with atmospheric CO2 and changes in global mean surface temperature.
Increasing atmospheric CO2 alone increases oceanic CO2 uptake whereas increasing temperature alone 10 decreases CO2 uptake. In the year 2100, the carbon-climate parameter is negative of -5.4 Pg C/K while the carbon-concentration parameter is positive of 0.79 Pg C/ppm. From 1850 to 2100, the carbon-climate parameter decreases with the increasing temperature change, indicating that with enhanced warming, each degree increase of surface temperature would induce more CO2 outgassing from the ocean (Fig. 16d). The Carbon-concentration parameter initially increases with atmospheric 15 CO2 and then decreases (Fig. 16c). The decreasing trend of is consistent with the slowdown of the increasing trend of the oceanic CO2 uptake at the end of the 21 st century as a result of decreased oceanic buffer ability due to the increasing DIC concentration. Similar trends of carbon-climate and carbonconcentration sensitivity parameters are also found in CMIP5 models . The increased sensitivity of CO2 outgassing to temperature and the decreased sensitivity of CO2 uptake to atmospheric 20 CO2 concentration indicate that the ocean's ability to absorb atmospheric CO2 would be weakened with increasing atmospheric CO2 and global warming.
In this section, we compared the carbon sensitivity parameters diagnosed from the 1ptCO2 experiment between the NESM v3 and CMIP5 models. The total CO2 uptake during the 140 years in FC-1% is 636 25 Pg C, while the results from CMIP5 models range from 533 to 676 Pg C. The sum of the total CO2 uptake in the RC-1% and the BC-1% is 65.7 Pg C, which is larger than that in the FC-1%. The simulated 24 nonlinearity (i.e. BC-1% + RC-1% -FC-1%) is about 10.3% of the total CO2 uptake in the FC-1%, which is at the higher end of the nonlinearity estimated for CMIP5 models ranging from 3.6%-10.6% (Schwinger et al., 2014).
Then, we compare NESM v3 simulated β and γ parameters with those of CMIP5 results. Figure 17 shows the simulated β and γ parameters in the 1ptCO2 runs. At the end of 1ptCO2 runs, the diagnosed value of β 5 from CMIP5 models ranges from 0.69 to 0.91 PgC/ppm with a multi-model mean value of 0.80 PgC/ppm.
For comparison, the β diagnosed from the NESM v3 simulations is 0.88 PgC/ppm at the end of the simulation. The declining trend in β is found after ~550 ppm, later than in Hist+RCP8.5 experiments (~400 ppm), which is consistent with the results in CMIP5 models. Compared with β, the γ from CMIP5 models has a much larger range, and the value at the end of the simulations ranges from -2.4 to -12.1 10 PgC/K. The larger spread of γ is associated with the spread of the model-simulated climate change and the dependency of carbon cycle processes on climate change. For comparison, our simulated γ parameter is -7.9 PgC/K.

Discussion and conclusion
In this study, we evaluate the performance of the NESM v3 in simulating the present-day ocean future scenario, the NESM v3 simulated the total amount of oceanic CO2 (567 PgC) is at the higher end of CMIP5 results (420~600 PgC). The sensitivity of oceanic CO2 uptake strongly depends on scenarios.
In the 1ptCO2 experiment, the NESM v3 simulated sensitivities of oceanic CO2 uptake, i.e., carbon- In the upper ocean, slight overestimations of nutrients are found in the Pacific and the Southern Ocean 25 (Fig. 1). In our simulations, the regions of overestimated nutrients (Fig.1) in general corresponds to regions with strong iron limitation (Fig. 4). The strong iron limitation in these areas limits biological activities, therefore reducing the uptake of nutrients by phytoplankton. Also, the overestimated nutrients 26 are probably associated with strong vertical mixing in the Pacific and the Atlantic, which is indicated by the dipole pattern of preformed phosphate bias above 1500 m depth (Fig. 3). In the Indian Ocean, the underestimation of nutrients is associated with the weak upwelling (figures not shown) that suppresses the nutrient entrainment to surface water. The low-level nutrients in the Indian Ocean reduce the biological activities and then result in underestimation of NPP and chlorophyll. Also, in a relatively coarse resolution 5 model, the underestimation of NPP and chlorophyll in the Indian Ocean could be associated with the poor descriptions of mesoscale and submesoscale processes (McGillicuddy et al., 1998;Lévy et al., 2001b).
The underestimation of alkalinity is simulated near the maritime continent, where the model underestimates the surface salinity due to excessive precipitation. In the high-latitude ocean, the model underestimates SST of about 1°C, indicating stronger convective mixing, which leads to the 10 overestimation of alkalinity.
As for vertical profiles of the biogeochemical field, the latitudinal-depth distribution of nutrients broadly agrees with observations, but the simulated high-concentration centers in the Northern Pacific are too strong and too deep (Fig. 2). Excessive remineralization is the main cause of the overestimated nutrients in the Northern Deep Pacific according to the results of nutrients decomposition (Fig. 3). Similar to the  Overall, compared with both observations and CMIP5 models, the NESM v3 does a good job in simulating ocean biogeochemical fields and oceanic carbon uptake. Despite these model-observation discrepancies, it is expected that NESM v3 can be used as a useful modeling tool to study interactive feedbacks between the ocean carbon cycle and climate change and the underlying mechanisms.
Code and data availability. 25 The      56 Figure 17. Same as Figure 16, but for the 1ptCO2 runs.