Impact of Increased Horizontal Resolution of an Ocean Model on Carbon Circulation in the North Pacific Ocean

The impact of resolving western boundary currents and mesoscale eddies on a carbon circulation simulation for the North Pacific Ocean is investigated to evaluate the merits of using high‐resolution ocean biogeochemical models for climate projections. Simulations by a 100‐km resolution global ocean biogeochemical model with and without embedding a 10‐km resolution model in the North Pacific Ocean are compared. The major improvement in the high‐resolution simulation is the representation of the Kuroshio, its extension current, and the recirculation gyres formed to its south and north, resulting in a proper representation of the North Pacific subtropical mode water (STMW) and an increase in storage of the anthropogenic CO2 (Canth) in STMW by about two‐thirds. The larger storage rate in STMW is caused by supply of a larger amount of warm surface water containing rich Canth to the formation region by the intensified Kuroshio. A huge buoyancy loss from this warm water results in the increased formation of STMW that occupies a vast area in the western subtropical gyre. The surface uptake of Canth in the formation region of STMW is slightly increased but is largely comparable to that of the low‐resolution model. Moreover, there is no structural difference in Canth uptake in other parts of the subtropical region. Thus, the improvement of Canth distribution can be understood as a redistribution of water mass in the subtropical gyre by the improved circulation. The present assessment motivates the use of a high horizontal resolution ocean model in next‐generation projection experiments with carbon cycles.

This study investigates the impact of using a high-resolution ocean model on C anth circulation and storage in the North Pacific Ocean with a view of adopting high-resolution ocean models for the next generation long-term projections of climate and ocean change and for the current activities to assess global carbon budgets (e.g., DeVries et al., 2023;Friedlingstein et al., 2022).Since the use of globally high-resolution ocean biogeochemical models for long-term projections would be unfeasible in the near future, a regional high-resolution model is used for the North Pacific Ocean and is embedded within a global low-resolution ocean model.Considering a possible use of such a global ocean model embedded with a high-resolution regional ocean model as an oceanic component of an Earth system model in the future projection experiments, a two-way conservative exchange of materials between these models is implemented because global conservation is an important element in longterm projections.It should be noted that successful application of such a system to a long-term simulation of ocean biogeochemical cycles is very rare.

10.1029/2023MS003720
3 of 38 As a first step, a long-term ocean biogeochemical simulation is conducted in an ocean-only mode, that is, the model is forced at the surface by a set of surface atmospheric variables based on observations and reanalyzes.By making use of the framework for the Ocean Model Intercomparison Project (OMIP; e.g., Tsujino et al., 2020), a long-term simulation for ocean biogeochemical cycles is first conducted using a CMIP6-class 100-km resolution ocean model.A 10-km resolution model is embedded in the North Pacific Ocean for the last 80 years and run in parallel with the simulation without embedding.The results of the last few decades are compared after the spin-up of physical and biogeochemical fields in the high-resolution model.Following a brief evaluation of longterm mean physical and biogeochemical fields for both simulations, a water mass transformation (WMT) framework (Groeskamp, Griffies, et al., 2019) is used to assess the column inventory, storage rate in layers separated by densities, and the regional budgets for C anth to clarify the processes that bring about the difference between the two simulations.It should be noted that there are still not many studies applying the WMT framework for carbon cycles (e.g., Aldama-Campino et al., 2020;Groeskamp et al., 2016;Iudicone et al., 2011;Zhai et al., 2017).The present study would be one of the first applications of this method to a model simulation of carbon cycles in the subtropical North Pacific Ocean, where much C anth has been accumulated, with eddying horizontal resolution.
The rest of the paper is organized as follows.Section 2 explains the model, experiment, and analysis method, with the formulation for the WMT framework given in Appendix A. Section 3 evaluates the general features of the simulated fields.Section 4 investigates changes in C anth circulation and storage in the North Pacific Ocean due to increased horizontal resolution of the model.Section 5 discusses the redistribution nature of the change in C anth distribution in the high-resolution model.Finally, Section 6 summarizes this study.The feasibility of the embedded system of global and regional ocean models for use in long-term simulations of ocean biogeochemical fields is assessed throughout this paper.

Model
We use the ocean component of the Meteorological Research Institute Earth System Model version 2 (MRI-ESM2; Yukimoto et al., 2019) as the base global ocean model.The source code is taken from the Meteorological Research Institute Community Ocean Model version 4 (MRI.COMv4; Tsujino et al., 2017), which is an ocean general circulation model formulated on a general orthogonal curvilinear coordinate (Murray, 1996) in the horizontal direction and on a z* coordinate (Adcroft & Campin, 2004) in the vertical direction and is discretized on the Arakawa B-grid frame.The global model uses a tripolar grid system to avoid singularity at the North Pole.The horizontal resolution is 1.0° in the zonal direction and 0.5° with an enhancement to 0.3° toward the Equator in the meridional direction.There are 60 vertical levels with an enhancement in the upper layer and an additional bottom boundary layer (Nakano & Suginohara, 2002) at the sea floor in the northern North Atlantic and in the Southern Ocean around Antarctica, where the deep and bottom waters are formed.Urakawa et al. (2020) described the configuration and assessed the performance of this global ocean model in terms of physical fields.The performance of ocean biogeochemical cycles of MRI-ESM2 in CMIP6 experiments can be found, for example, in Kwiatkowski et al. (2020), Planchat et al. (2023), Séférian et al. (2020), andVaittinada Ayar et al. (2022).
We embedded a high-resolution regional model for the North Pacific Ocean into the base global model to test the impact of increasing horizontal resolution of an ocean model on the biogeochemical cycles in the North Pacific Ocean.The base system for the present experiment has been developed for a coastal ocean monitoring and prediction system around Japan (Sakamoto et al., 2019).The domain for the regional model is 99°E-85°W zonally and 15°S-63°N meridionally, which is depicted as the rectangle in Figure 1d.The horizontal resolution is 1/11° and 1/10° for zonal and meridional directions, respectively.The same 60 vertical levels as the global model are used.At the boundary between the models, all prognostic variables are exchanged between the models to fill two grid points outward from the boundary of both models.Each model is then advanced in time with a common time interval.In this procedure, the horizontal fluxes of tracers and water volume at the side boundary of the high-resolution model are adjusted to match that of the low-resolution model, conserving the total tracer mass and water volume of the system.By imposing such a conservative exchange, the embedded model can be suitable for ocean-climate studies requiring long-term simulations.
The same set of schemes is used for the dynamical frame of both high-and low-resolution models.The advection scheme for momentum realizes diagonally upward and downward momentum advection at the bottom relief and 10.1029/2023MS003720 4 of 38 the conservation of pseudo-enstrophy, which is particularly suited for high-resolution simulations as demonstrated by Ishizaki and Motoi (1999).A second order moment conservation scheme by Prather (1986) is used for temperature and salinity advection.Regarding the subgrid scale mixing, the generic length scheme by Umlauf and Burchard (2003) is used to express the local turbulent vertical mixing.The vertical mixing coefficient for tracers is enhanced to 1.0 m 2 s −1 for gravitationally unstable vertical stratification.The background vertical mixing coefficient for tracers, which accounts for mixing due to breaking of internal tides, is given following Decloedt and Luther (2010).Horizontal viscosity is Laplacian and biharmonic type for low-and high-resolution models, respectively.The isopycnal diffusion (Redi, 1982) and the Gent and McWilliams (1990) parameterization for eddy-induced transport expressed as skew flux in the diffusion tensor as introduced by Griffies (1998) are used for the whole domain in the low-resolution model and only around the side boundaries in the high-resolution model.
The biogeochemical processes consist of a carbon cycle model with the carbonate chemistry and the surface gas exchange parameterization following the OMIP-BGC protocols (Orr et al., 2017) and a simple nutrient, phytoplankton, zooplankton, detritus-type ocean ecosystem model as used by Nakano et al. (2011 and2015).It should be noted that this study returned to a slower detritus sinking velocity (2 m day −1 ) originally used by Nakano et al. (2011) and for the CMIP5 experiments from a faster one (7 m day −1 ) used for CMIP6 by following Schmittner et al. (2009), to realize more surface CO 2 outgassing in the eastern equatorial Pacific Ocean as observed.A relatively low-cost scheme (piecewise parabolic method; Colella & Woodward, 1984) is used to advect biogeochemical tracers.

Experiment
Basically, the OMIP-BGC protocols for ocean biogeochemical simulation (Orr et al., 2017) were followed with some minor modifications.Taking 366 years as a unit cycle, the ocean biogeochemical fields were integrated for four cycles (1,464 years) as a spin-up.The surface atmospheric state was taken from a repeat year from 1 May 1990 to 30 April 1991 (Stewart et al., 2020), of the JRA55-do surface data set (Tsujino et al., 2018).The atmospheric CO 2 concentration was kept at a preindustrial (year 1850) value of 284.32 ppm.The physical state of the ocean was started from a state of rest with temperature and salinity from WOA13v2 (Locarnini et al., 2013;Zweng et al., 2013) and reset to this initial state at the beginning of each cycle (every 366 years).This treatment was performed to prevent the physical state from being drifted farther from the observation-based initial state during the long-term spin-up.On the other hand, the biogeochemical fields are continuous during the entire integration, with the initial state taken from GLODAPv2 (Lauvset et al., 2016) for preindustrial DIC and alkalinity, and from WOA13v2 for oxygen (Garcia et al., 2014a) and nutrients (Garcia et al., 2014b).The last 366 years followed the OMIP-2 protocol (Tsujino et al., 2020).The model was forced by repeating the interannually varying atmospheric state from the 61-year JRA55-do data set  for six times (366 years; nominally starting from year 1653 and ending in 2018).Therefore, the total integration years for the biogeochemical state is 1,830 years.Figure 1a depicts the globally integrated annual mean surface CO 2 uptake during the entire integration.The drift keeps decreasing, but a small CO 2 loss with its absolute value <0.2 PgC year −1 remains (i.e., the total natural carbon content in the model that sees a preindustrial atmospheric CO 2 concentration is decreasing).In the present simulation, the formation of North Atlantic Deep Water (∼14 Sverdrups [Sv]; 1 Sv = 10 6 m 3 s −1 ; Figure S1 in Supporting Information S1) was smaller than observation (∼18 Sv; Smeed et al., 2019) and that the southward transport of CO 2 would be smaller than estimated in the Atlantic Ocean, which could be one of the reasons for the global carbon content to slowly decrease.Considering that the main purpose of this study is to investigate the impacts of increasing horizontal resolution of an ocean model in the North Pacific Ocean, the issue of small drift is not serious.This issue is left for future studies.with xCO 2 that starts to increase from year 1850 for the global model without embedding and (b) after the embedding of the high-resolution North Pacific Ocean model from year 1938.Details about the experiment are described in Section 2.2.In (a), surface uptake during the spin-up and the last 366 years  for the natural dissolved inorganic carbon (DIC) tracer that sees the constant xCO 2 of 284.32 ppm is depicted with the blue line.The one for the contemporary DIC that sees an increasing xCO 2 is depicted with the red line (also in (b)).In (b), surface uptake for the embedded run is depicted with the darkblue line.Means of observationbased products (green) and ocean biogeochemical models (yellow) provided by Global Carbon Project (2022) are also depicted with uncertainty bounds represented by ±1 standard deviation in shades.(c-f) Long-term mean , annual mean surface uptake of CO 2 (mol m −2 year −1 ) for the global model (c) without and (d) with embedding, (e) difference between with and without embedding, and (f) an observation-based estimate (JMA-MLR) by Iida et al. (2021).In (d) and (e), the embedded region is represented by the purple rectangle.For (c) and (d), root-mean-square error and mean bias north of 12°N relative to JMA-MLR are put on top of panels.The dashed rectangle in (c)-(f) is used for a quantitative assessment of simulations (see text in Section 3).

Methods for Analysis
We will specifically investigate the change in C anth uptake and storage due to increasing horizontal resolution of an ocean model in this study.For this purpose, a natural DIC tracer that sees a constant preindustrial reference atmospheric mole fraction value of CO 2 (xCO 2 = 284.32ppm) is included in addition to the contemporary DIC tracer that sees an increasing atmospheric CO 2 concentration.Here C anth is evaluated as the difference between the contemporary and natural DIC tracers.It behaves as a passive tracer (Sarmiento et al., 1992).The analysis of Nakano et al. (2015) was followed, and a budget analysis of C anth was conducted for the subdivided regions of the North Pacific Ocean for vertical layers separated by potential density (σ 0 ) at 22.0, 23.0, 24.0, 24.5, 25.0, 25.5, 26.0, 26.5, 27.0, and 27.5.The formulation for the present budget analysis is based on de Szoeke and Bennett (1993) and Groeskamp, Griffies, et al. (2019) and given in Appendix A. The budgets for the rate of C anth storage for a density class bounded by ρ u (z) < ρ < ρ l (z) (with u and l representing upper and lower, respectively) over a fixed region are considered based on the following formula for a conservative tracer ϕ as used by Nakano et al. (2015): (1) (Units: m 3 • s −1 • (units of tracer ϕ)), where tildes denote that the calculation is performed on density coordinates.In Equation 1,  F represents thickness (z ρ )-weighted tracer flux on density coordinates (Units: where (u, v) where ω tot and ω gm are diapycnal fluxes (Units: kg m −3 • s −1 ) due to total velocity field (resolved plus parameterized) and Gent-McWilliams scheme, respectively, and  F vdiff is flux due to vertical diffusion.As shown in Equation 3, the conversion term involves many terms and sources of errors.Therefore, its value can be alternatively estimated as a residual using other three terms, namely, storage rate, transport convergence, and surface flux, as given by the left-hand side, the first term on the right-hand side, and the last term on the righthand side of Equation 1, respectively, as in Nakano et al. (2015).It should be noted that the residual also contains errors originating from those three terms.As an extension to the analysis of Nakano et al. (2015), the components of the conversion term in Equation 3 are evaluated explicitly to understand the reason for the differences caused by increasing the horizontal resolution of the model.
Diabatic processes, such as surface buoyancy flux, vertical diffusion for density, and lateral mixing (cabbeling, thermobaricity) cause the part of the conversion term due to the convergence of diapycnal water mass transport represented by the first term on the right-hand side of Equation 3. The diapycnal water mass transport is split into its components as which can be quantified separately using sampling of explicit mixing terms and surface fluxes for temperature and salinity.Those components can be used to split the first term on the right-hand side of Equation 3 into its components as (5) The computation of the terms on the right-hand side of Equation 3 should be tested by comparing their sum with the top-down quantity calculated as residual by using Equation 1.The difference, or errors, can be caused by various numerical treatments during the processing and by the numerical diffusion inherent in tracer advection schemes used by the model.The overall computation errors will be evaluated as the difference between the conversion terms estimated as the residual using Equation 1 and as the sum of its components evaluated explicitly using Equations 3 and 5.
Finally, to consider the budgets of the column inventory, Equation 1 is integrated vertically over the whole water column.The conversion term becomes zero upon vertical integration since it involves only internal exchange processes, therefore, the following formula is obtained for time evolution of the column inventory: The analysis is conducted using 5-day mean diagnostics sampled during runtime.The potential density (σ 0 ) is calculated first, and the depths of the isopycnal surfaces are calculated by interpolation to density levels finer than those used for considering budgets as listed above.Thickness weighted averaging is then used to conservatively partition the diagnostic quantities in the vertical column on the z* coordinate into density layers.The diagnostics range from the basic prognostic variables to their tendency and transport terms, reducing the need for offline diagnostic calculations using the model outputs and the model source code.
Simulation results are evaluated by comparing with observation-based estimates first, then the simulations without and with embedding, differing in the horizontal resolution in the North Pacific Ocean, are compared.We mainly assess whether important features in the western North Pacific are qualitatively improved by shifting to the high-resolution.Though it may not seem quantitative, this approach would be justified because the main purpose of this study is to investigate the impacts of increasing horizontal resolution of an ocean model.
To quantitatively grasp model biases, area mean root-mean-square errors and biases in simulations relative to observation-based estimates for the North Pacific Basin north of 12°N are shown on top of panels.In addition, focused comparisons are made for surface CO 2 flux and C anth inventory.The issues of reducing overall model biases are left for future studies.
A 31-year period from 1987 to 2017 is adopted to assess long-term mean fields and budgets of inorganic carbon variables unless otherwise noted since the observation-based DIC and C anth distribution from GLODAPv2 used for evaluation is normalized to the year 2002.The 31-year period spans 15 years equally before and after the year 2002 and the 31-year averaging is intended to suppress the effects of decadal variability on the assessment of the mean fields as much as possible.
For comparison, model diagnostics from both the low-and high-resolution models are regridded on a standard 1° latitude and 1° longitude grid.Most reference data sets used in this study are provided on the standard 1° × 1° grid and are readily used for model evaluation, except for the sea surface height (SSH) provided by Copernicus Marine Environment Monitoring Service which is daily and on 0.25° latitude × 0.25° longitude grid.The SSH data is averaged for each month and then regridded on the standard 1° × 1° grid using a Gaussian filter with a half width of 1.5° as in Tsujino et al. (2020).

Difference of Simulated Fields Due To Increased Horizontal Resolution of the Model
The globally integrated surface CO 2 fluxes show no notable difference by embedding the high-resolution model in the North Pacific Ocean (Figure 1b).Our simulations are left behind the recent rapid increase of the surface uptake seen in the observation-based products.However, the absence of the rapid increase is also found for the majority of ocean biogeochemical models that contributed to GCB 2022 (Friedlingstein et al., 2022).Overall, our model can be regarded to show a standard behavior compared to the majority of models in terms of the globally integrated surface CO 2 flux.
Surface CO 2 fluxes of the high-and low-resolution models differ locally in the embedded region (Figure 1e).Somewhat by chance, these differences seem to compensate one another upon integration over the embedded region.Notable differences are found in the subpolar (negative, outgassing), the western to central midlatitude (positive, uptake), and the eastern tropical (positive) regions.Some of these differences tend to reduce the biases of the low-resolution model, that is, the unrealistic uptake in the subpolar region and the smaller uptake in the western to central midlatitude region (Figures 1c-1f).In the region of the most rigorous CO 2 uptake in the midlatitude North Pacific in the observation-based estimate (dashed rectangle in Figure 1f), the mean uptake is 3.11 mol m −2 year −1 , while it is 2.71 and 3.27 mol m −2 year −1 in the low-and high-resolution model, respectively.Though the simulated biases relative to the observation-based estimate are within the uncertainty range of the observation-based estimate (∼25%; Iida et al., 2021), the bias is reduced from −13% to +5% by shifting to the high-resolution.These changes are explained by those of basic fields as described in the following.
Physical fields are notably different between the high-and low-resolution models in the midlatitude western North Pacific (Figures 2c,2g,2k,and 2o).Both the low Sea surface temperature (SST) bias in the western to central midlatitude region and the high SST bias to the east of Japan in the low-resolution model almost disappear in the high-resolution model (Figures 2a-2c) due to the improvement of the Kuroshio and its extension as represented by the SSH fields (Figures 2e-2g).Specifically, a proper representation of the twin relative recirculation gyres in the Kuroshio Extension region, as visualized by the pair of zonally elongated negative and positive anomalies along about 35°N in the difference of SSH (Figure 2g), is a prerequisite for proper representation of physical fields in this region (Nakano et al., 2008).It has a major impact on the representation of STMW (Hanawa & Talley, 2001) as reported by Nishikawa et al. (2010).The broad region of deep winter mixed layer depth (MLD) in the low-resolution model to the east of Japan is separated into north and south in the high-resolution model (Figures 2i-2k The partial pressure of CO 2 (pCO 2 ) of the surface water (Figure 3) and the DIC distribution along 164.5°E (Figure 4) explain the changes of surface CO 2 fluxes by shifting to the high-resolution shown in Figure 1.The higher DIC concentration at the surface (Figure 4c) explains the higher pCO 2 in the subpolar region in the high-resolution model (Figure 3c).The lower pCO 2 to the north of the Equator in the eastern tropical Pacific is due to the spread of the low pCO 2 water from the north by the tropical instability waves as indicated also by the higher SST in the high-resolution model (Figure 2c).The lower pCO 2 in the Kuroshio Extension region is due to the lower DIC, which is caused by the increased transport of low DIC and high alkalinity (Alk), that is, high buffer capacity subtropical water by the Kuroshio in the high-resolution model (Figures 4c and 4g).
General features of the observed distribution of vertically integrated net primary production (intPP), such as high productivity around the equator and in mid-latitudes, are reproduced in both simulations (Figures 3e-3h).Differences between the high-and low-resolution models of pCO 2 and intPP show similar features.Since a reduced primary production would generally increase surface DIC and therefore pCO 2 , the difference fields imply that difference of intPP does not significantly impact that of pCO 2 in our simulations.
After the long-term spin-up for more than 1,500 years, DIC and Alk appear to have large biases (Figure 4).The negative DIC and Alk biases in the thermocline depths north of 10°N are particularly notable.However, it should be noted that even state-of-the-art Earth system models show biases in DIC and Alk that locally reach up to 0.1 mol m −3 (e.g., Yool et al., 2020).The bias in our simulations is considered to be due to the weak upwelling of the mid-depth to deep water in the North Pacific Ocean.The simulated upwelling from the mid-depth to the surface is indeed generally weak in the North Pacific Ocean in the meridional overturning stream function (10°-60°N and 1,000-3,000 m) (Figure 5a).Upwelling is intensified above 1,000 m of the subpolar region (40°-55°N) in the high-resolution model (Figure 5b), and the negative DIC and Alk biases are reduced in the thermocline depths (Figures 4c and 4g).Although this may simply imply that the choice of parameter for the Gent-McWilliams parameterization is not appropriate for the low-resolution model as mentioned earlier, the positive aspect would be that the process relevant to the improvement (i.e., upwelling) is identified using the high-resolution model that can express the relevant processes explicitly (i.e., without parameterization).The present simulations tend to have higher DIC and Alk in the deep to bottom layer, where observation implies inflow of low DIC and Alk water from the south.Again, the bias of the model is presumably due to the weak deep meridional overturning circulations of the Pacific Ocean of about 5 Sv compared to an observed value of 8 Sv at 10°N (Wijffels et al., 1996).The deep circulation is slightly stronger by about 1 Sv in the high-resolution model, which is presumably due to the improved representation of the bottom topography, but this does not seem to have significant impact on distribution of biogeochemical variables except for Alk.On the other hand, the improvement of vertical distribution of these biogeochemical tracers can be also achieved by adjusting biological processes.This is another important aspect of ocean biogeochemical modeling, but consideration about correcting ocean biochemical fields by improving biological models is beyond the scope of this study.
Both the negative and positive biases of nitrate distribution in the subpolar and subtropical upper oceans, respectively, in the low-resolution model is improved in the high-resolution model (Figures 4i-4k).Since the availability of simulated nitrate largely restricts the primary production in our biogeochemical model, this change in In summary, many improvements, at least qualitative, may be expected in simulation of both physical and biogeochemical fields by taking a high horizontal resolution for the North Pacific Ocean.Also, no notable discontinuity is found around the boundary between the embedded high-resolution regional model and the outer low-resolution global model (e.g., Figures 1, 4, and 5, with the boundary depicted as the purple rectangle in Figure 1d), showing that embedding can be a useful tool.Furthermore, managing an 80-year integration without notable problems for the set of models with a two-way exchange demonstrates the ability of such a model system for use in long-term projections of climate and ocean change.
4. Uptake, Distribution, and Storage of Anthropogenic CO 2

General Features
The column inventory of C anth is large in the northwestern subtropical region in GLODAPv2 (Figure 6d), but negative bias is found there in the simulations (Figures 6a and 6b).This bias is reduced in the recirculation gyre of the Kuroshio and its extension in the high-resolution model (Figures 6b and 6c).The column inventory in the high-resolution model is smaller than the low-resolution model on both sides of the Equator.This is due to In the first and second column, bias relative to the observation-based estimates is depicted for the low-and high-resolution model, respectively.On top of panels, root-mean-square error and mean bias north of 12°N relative to the observation-based estimates are put.In the third column, difference between the high-and low-resolution models is depicted.In the fourth column, the reference data based on observations is depicted.DIC and Alk from GLODAPv2 and NO 3 from WOA13v2.
the difference in the representation of the subtropical ridge: the old water with high DIC concentration is lifted close to the sea surface in the high-resolution model, disturbing C anth uptake and storage in the surface layer in this latitude band (Figure 4c).The storage rate (Figures 6e-6h) basically shows the same feature as the total column inventory except for the patchy features caused by taking the difference between instantaneous fields (first of January between 2018 and 1987), implying that the total column inventory generally reflects the recent trend of the C anth storage.As shown in Figure 6l, the latitudinal peaks in the zonally integrated storage and surface uptake rates of C anth in the Pacific Basin do not coincide, implying the importance of redistribution by the ocean circulation.Compared with the storage rate, which shows the shift toward the higher latitude in the subtropics (15°-35°N) in the high resolution model, the latitudinal distribution of surface uptake rate does not change remarkably by shifting to the high-resolution except for the smaller uptake on both sides of the Equator and the larger uptake in northern high latitudes by the high-resolution model.
In the latitude-depth section of C anth along 164.5°E (Figure 7), the increase of C anth concentration for the latitude band around 30°N corresponds to that of the column inventory in the high-resolution model.Peaks of C anth concentration increase are found for the density of STMW (σ 0 = 25.0-25.5)and the denser water (σ 0 = 26.0-27.0).The density of STMW in the high-resolution model is lower because it is formed in the lower latitudes than the low-resolution model owing to the separation of the Kuroshio from the eastern coast of Japan and the formation of the intense Kuroshio Extension at the correct latitude.The Kuroshio Extension also clearly separates the The difference in the regionally integrated column inventory between the high-and low-resolution models keeps increasing in the northwestern subtropical region (Figures 8a and 8b).Thus, the difference between the averages for 1987-2017 shown in Figure 6c is not affected by variability on decadal to multi-decadal timescales and can be regarded robust.The column inventory in the northwestern subtropical region is increased by about 10% by shifting to the high-resolution.When partitioned into density classes, the contribution to the total difference is largest from the STMW density class (σ 0 = 25.0-25.5),and the storage in this density class is increased by about 66%.In the southwestern subtropical region, there is almost no difference in the column inventory, but the storage in the STMW density class is increased by 40%.Instead, the storage is reduced in the upper (σ 0 = 23-24) and the lower (σ 0 = 26-27) layers.The lower layer corresponds to the CMW density range.
The time series of C anth storage in the STMW density class of the low-and high-resolution model clearly separates from each other in both northern and southern subtropical regions (Figures 8f and 8g).The volume (Figures 8h  and 8i) of this density range increases by shifting to the high-resolution in both regions.The initial adjustment is largely completed in the first 20 years after the start of the high resolution model.After the adjustment, the TSUJINO ET AL. 10.1029/2023MS003720 14 of 38 volume varies with decadal timescales but no notable drift is found, implying that the differences of C anth storage between the high-and low-resolution models are largely caused by the differences of water volume.The decadal variabilities in the low-and high-resolution models are in phase, with the high-resolution model varying with a larger amplitude.These features of volume variability are reflected in the variability of C anth storage.The differences in the C anth storage between the low-and high-resolution models start to increase from around 1980 and keep increasing to the present.This is caused by the increase of C anth concentration (Figures 8j and 8k).
In the horizontal distribution within density layers (Figures 9 and 10), C anth is accumulated mainly in mode waters over the subtropical region, which are characterized by large thicknesses, in both models, that is, STMW in the western subtropical gyre in σ 0 = 25.0-25.5 and CMW in the northern subtropical gyre in σ 0 = 26.0-26.5.The increase in C anth storage in the STMW layer by shifting to the high-resolution is notable (Figure 9c).It occupies the western subtropical gyre, resulting in the smaller accumulation in the lower layer of the high-resolution model: C anth in CMW is displaced to the outskirt of the subtropical gyre (Figures 10a-10c).These changes can be also largely explained by those of layer thicknesses.Note also that more C anth spreads southeastward from the eastern coast of Japan to the northwestern part of the subtropical gyre in σ 0 = 26.0-26.5 in the high-resolution model (Figure 10c).Since the difference of thickness between the high-and low-resolution models are small (Figure 10f), this feature is understood to be due to the higher C anth concentration in the high-resolution model.Spreading of the recently ventilated water from the northern side of the Kuroshio is suggested, which is consistent with Nakano et al. (2021).This is also seen in the latitude-depth section, contributing to the increase of column inventory in the high-resolution model around 30°N (Figure 7c).
Overall, in the region of the main C anth storage in the midlatitude North Pacific in GLODAPv2 (dashed rectangle in Figure 6d), the mean storage is 45.9 mol m −2 , while the simulated bias is −18.5 and −13.8 mol m −2 in the low-and high-resolution model, respectively.Though the simulated biases relative to GLODAPv2 are outside the  uncertainty range of GLODAPv2 (0.7-1.3 mol m −2 for 1,000 m water column; Lauvset et al., 2016), the bias is reduced from −40% to −30% by shifting to the high-resolution.

Budget Analysis
Results of a WMT anaylsis (Groeskamp, Griffies, et al., 2019) for the separated density classes in the northwestern and southwestern subtropical regions are summarized in Figures 11 and 12, respectively.The column integrated values relevant to time evolution of the column inventory (Equation 6) are put on top of panels and their geographical view is shown in Figure 13.The rate of increase of the column inventory in the northwestern subtropical region in the high-resolution model (50.6 TgC year −1 ) is higher by about 10% than the low-resolution model (45.9 TgC year −1 ).This relative fractional increase due to shifting to the high-resolution is larger than any other regions (Figure 13c).Among the possible sources of this larger rate, the inflow from the southwestern subtropical region through the section along 28°N is the most notable.Although the outflows through other boundaries are also larger in the high-resolution model, the difference of the vertically integrated transport convergence (Figure 11c) largely explains that of the storage rate for this region.However, considering the large difference in the surface heat flux (Figure 2o), it would be an unexpected result that the surface uptake (Figure 11b) does not explain the difference.The surface uptake is rather small in the high-resolution model.
Notable changes occur for the budget of the C anth storage rate (Equation 1) in the density class of STMW (σ 0 = 25.0-25.5) in both the northwestern and southwestern subtropical regions (see Figures 14c and 14d for specific values).The rate of increase of C anth storage in the STMW layer due to conversion in the northwestern subtropical region is about three times larger in the high-resolution model than the low-resolution model.Most of this excessive formation enters the southwestern subtropical region.The storage rate in the northwestern subtropical region in the high-resolution model is higher by about 30% than the low-resolution model.The loss due to conversion in the southwestern subtropical region is about six times larger in the high-resolution model.The storage rate is higher by about 80%.The general increase of transports through the boundaries between the separated regions in the subtropics is due to the intensified Kuroshio recirculation gyre in the high-resolution model (e.g., Figures 2g, 9a, and 9b).
The increase of the surface uptake due to increasing the horizontal resolution alone cannot explain that of the storage rate for the STMW density class (σ 0 = 25.0-25.5) in the northwestern subtropical region (Figures 11b,  14c, and 14d).It could be possible to understand that the difference in the surface uptake (1.4 TgC year −1 ) explains about 70% of the increase in the storage rate (2.0 TgC year −1 ); however, for this explanation to be correct, all surface uptake needs to be stored in this region.This explanation is unlikely, given that the conversion rate is three times larger (Figure 14d).Indeed, the convergence of diapycnal transport and the surface uptake take place in the same region, that is, south of the Kuroshio Extension (Figures 15e-15h).Thus, it would be more appropriate to understand the result as follows.In the high-resolution model, the conversion and surface uptake contributed to the C anth contents in the newly formed water (99.2TgC year −1 = conversion + surface uptake) by 79% and 21%, respectively.About 8% (8.3 TgC year −1 ) of this newly formed C anth is stored in this region and is used to increase the C anth concentration of this layer, assuming that the volume has almost no trend.In the same manner, in the low-resolution model, the conversion and the surface uptake contributed almost equally to the C anth contents in the newly formed water (42.2TgC year −1 ), and about 15% (6.3 TgC year −1 ) is stored within this region.The C anth concentration is similar for both models (Figure 8j) since the volume of this layer is larger in the   TSUJINO ET AL. 10.1029/2023MS003720 22 of 38 high-resolution model (Figure 8h).A consideration presented here is for the budget of the STMW density class in the entire northwestern subtropical region, which is not necessarily filled with STMW.The next section provides some detailed discussion focusing on C anth budget for STMW recirculating in the western subtropical region.
In the transport convergence for the northwestern subtropical region (Figure 11c), the excess outflow of C anth in the density range σ 0 = 25.0-26.5, which includes STMW (σ 0 = 25.0-25.5), is compensated by the excess inflow in the lighter density range, specifically in σ 0 = 23.0-23.5.Inflow at 28°N (Figure 11e) is indeed large for this density range, implying that this is caused by the increase of transport of the Kuroshio in the East China Sea by increasing the horizontal resolution.However, no notable difference is found in the surface uptake between the low-and high-resolution models in the density range lighter than the STMW density class in the northwestern subtropical region (i.e., north of 28°N; Figure 11b).Thus, the source for the change of the storage rate in the STMW layer by increasing the horizontal resolution is supplied by the enhanced transport of C anth already acquired by the low density upper layer water in low latitudes (Figure S3 in Supporting Information S1).
The conversion rate (Figure 11d) is split into its components as formulated by Equation 3and the conversion rate due to convergence of diapycnal transport caused by explicit mixing is further split into its components as given by Equation 5, and they are shown in Figures 11i-11m and 11n-11p, respectively.The conversion rate in the density range of STMW is dominated by the convergence of diapycnal transport caused by the surface buoyancy flux (Figures 11i and 11o).The diapycnal transport due to the surface buoyancy flux corresponds to the formation of STMW.As noted earlier, the increased transport of warm water by the Kuroshio and the larger heat loss to the atmosphere in the high-resolution model (Figure 2o) explain the large buoyancy loss and thus the diapycnal volume transport from the low-to high-density layer.The diapycnal transports due to vertical diffusion in the high-resolution model (Figure 11n) and the Gent-McWilliams scheme in the low-resolution model (Figure 11j) result in negative conversion rates, partly offsetting the formation due to convergence of diapycnal transport caused by surface buoyancy flux in both models.The contribution from the Gent-McWilliams skew flux term plus the isopycnal diffusion term for C anth to the STMW density range (Figure 11l) is relatively minor in the low-resolution model.The overall effect of vertical diffusion term for C anth evaluated by taking the difference between vertical diffusion terms for contemporary DIC and natural DIC is to diffuse high concentration in the lighter water close to the surface to the denser layer away from the surface (Figure 11k), with no notable difference between the low-and high-resolution models.As expected, the contribution of the isopycnal diffusion term in the low-resolution model to the convergence of diapycnal transport is relatively minor (Figure 11p).
In the southwestern subtropical region, more conversion from the STMW density class (σ 0 = 25.0-25.5)to the lighter classes is found in the high-resolution model (Figures 12d,15a,and 15b).This is caused by the diapycnal transport due to vertical diffusion (Figure 12m).The converted water is transported southward and eventually contributes to the C anth convergence into the density class of σ 0 = 23.0-24.0,which is comparable to the convergence of diapycnal transport caused by the surface buoyancy flux in the high-resolution model (Figures 12m  and 12n).The upward flux of C anth is seen at the latitude band of 20°-25°N at σ 0 = 24.0 in the high-resolution model and this is more than that of the low-resolution model (Figures 16c and 16d).Eddy kinetic energy is high in the latitude band of around 20°N where the eastward flowing, baroclinically unstable, subtropical counter current (Qiu, 1999) exists in the high-resolution model (Figure S2 in Supporting Information S1).Combined with the downward flux from the upper surface σ 0 = 23.0 due to the surface buoyancy loss (Figure 16b), this results in the more convergence for σ 0 = 23.0-24.0(Figures 16e and 16f).Part of the converted water in σ 0 = 23.0-24.0flows out to the tropics through 12°N while the majority flows into the northwestern subtropical region through 28°N (Figures 12e and 12f, respectively).Although the formation, or conversion, rate of C anth in σ 0 = 23.0-24.0 is larger in the high-resolution model (∼100 TgC year −1 ) than the low-resolution model (∼60 TgC year −1 ), the storage rate in this density class is slightly higher in the low-resolution model (10.5 TgC year −1 and 9.5 TgC year −1 for the low-and high-resolution models, respectively).In the high-resolution model, much of the newly formed C anth -rich water is transported to the north and used as the source water for STMW rather than being stored locally.
It should be noted that the part of the conversion rate that is not explained by the explicit diffusion terms (Figures 11m and 12l), which is evaluated as residuals, is much smaller than the dominant terms such as convergence of diapycnal transport (Figures 11i and 12h).The numerical errors are generally larger in the high-resolution model and the distribution of signs across density classes is similar to that of the convergence of diapycnal transport due to vertical diffusion (Figures 11n and 12m).This implies that the numerical errors are presumably due to the numerical diffusion inherent in the tracer advection algorithm caused by large velocities in the high-resolution model, which is the common problem in depth-coordinate high-resolution models (e.g., Ilicak et al., 2012).

Discussion
The present simulation and analysis have shown that a larger amount of C anth is stored in the northwestern subtropical North Pacific in the high-resolution model, which is closer to observation-based estimates.Moreover, this change in storage does not accompany notable structural changes in the distribution and magnitude of surface uptake of C anth (Figure 6k), although a large difference is found for the surface heat flux (Figure 2o).This difference between the changes of surface heat and C anth fluxes is basically due to the difference of typical timescales for the adjustment of the surface state, which are estimated about 1 month and 6 months for heat and CO 2 , respectively, for a surface mixed layer of 40-m depth under the surface wind with a speed of 7.5 m s −1 (Sarmiento & Gruber, 2006).An insufficient local supply of C anth at the sea surface to the enhanced WMT due to the enhanced buoyancy loss is compensated by the enhanced transport of C anth -rich water from low latitudes by the Kuroshio.
In other words, the present improvement in C anth distribution should be understood as horizontal and vertical redistribution within the basin.This section discusses how this change of C anth is achieved and maintained.

C anth Budget for the STMW Recirculation
In the western North Pacific, C anth is mainly accumulated in STMW.It is well known (e.g., Suga & Hanawa, 1995) and shown in the present simulation that STMW is accommodated in the recirculation gyre that is formed south of the Kuroshio and circulates on almost closed stream lines (Figures 9a and 9b).Here we consider C anth budget for the STMW in the recirculation gyre (the STMW recirculation) in the high-resolution model using Figure 14d.The eastern (170°W) and southern (12°N) boundaries of the southwestern subtropical region are open.The contour lines of the Montgomery potential (Figure 9b) imply that all of the southward transport at the southern boundary (39.0 TgC year −1 ) has come through the eastern boundary.The recirculation component at the eastern boundary is estimated as the residual (11.9 TgC year −1 ).The same amount is considered to flow eastward across the eastern boundary of the northwestern subtropical region in the recirculation, which leaves 14.3 TgC year −1 as the component flowing outside the recirculation gyre.Using this value, the sum of transports exiting the northwestern subtropical region without being imported to the STMW recirculation is estimated to be 19.2TgC year −1 by considering the exiting transports through the Tsushima Strait (9.7 TgC year −1 ), the Tsugaru Strait (−3.9 TgC year −1 ), and the western part of the 42°N section (−0.7 TgC year −1 ).Furthermore, the fraction (γ) of the total C anth input that is imported to the STMW recirculation is defined as follows: where with the right-hand side defined by the third and fourth terms on the right-hand side of Equation 1.The subscript NW denotes the northwestern subtropical region.For the high-resolution model, that is, about 80.6% (80.0 TgC year −1 ) of the total C anth input (99.2 TgC year −1 ) in the northwestern region is imported to the STMW recirculation.The total C anth export from the STMW recirculation is estimated as the sum of the export due to diapycnal transport which is dominated by the upward transport caused by vertical diffusion and the surface uptake in the southwestern subtropical region: TSUJINO ET AL.Thus, the storage rate for the STMW recirculation is estimated as follows: (Canth storage rate) STMW recirc.= (total Canth input) − (total Canth export) = 80.0 − 57.6 = 22.4 TgC year −1 . (13) This is quite consistent with the simulated storage rate (22.4(=8.3+ 14.1) TgC year −1 ) in the northwestern and southwestern subtropical regions.The same consideration can be applied to the low-resolution model (Figure 14c), which gives γ ∼ 0.562, that is, 56.2% (23.7 TgC year −1 ) of the total input (42.2(=22.5 + 19.7) Tg C year −1 ) is imported to the STMW recirculation.With the surface uptake of 0.6 TgC year −1 in the southwestern subtropical region, the total input of 24.3 TgC year −1 into the northwestern and southwestern subtropical region explains 15.1(=6.3+ 8.8) TgC year −1 and 9.2 TgC year −1 of storage rate and export, respectively, very well.
A simplified C anth budget model for the STMW recirculation in the western subtropical region is formulated by extracting the essence of the above quantitative estimation deduced from the WMT analysis.By expressing the import (formation) and export (dissipation) diapycnal volume transport for the STMW recirculation using a common transport value (I stmw ) to keep its volume, C anth concentration of the upper layer water imported to the STMW layer due to density conversion as   import anth , and that of the STMW layer water exported to the upper layer due to density conversion as   export anth , the import and export C anth transports are expressed as follows: and The surface uptake of C anth by the STMW recirculation is the fraction γ as defined by Equation 7 of the surface flux through the outcrop in the northwestern subtropical region that is calculated as the fourth term on the right hand side of Equation 1, which is symbolically represented as   stmw anth , The fractions (γ) are 0.806 and 0.562 for the high-and low-resolution models, respectively.Using these sources and sinks, the C anth storage rate with mean concentration   mean anth in the STMW recirculation, whose volume is V, is roughly expressed as follows: This simplified budget model is applied to the high-resolution model and compared with the budget for the STMW recirculation deduced from the WMT analysis in the previous paragraph using Table 1.The first term on the right-hand-side, that is, a net input to the STMW recirculation from the overturning circulation, is estimated to be 4.8 TgC year −1 , by taking I stmw ∼ 2.8 Sv from the amount of transport from the STMW layer to the surface layer in the southwestern subtropical region (Figure S5 in Supporting Information S1) and by estimating the value in the parenthesis to be 54.7 − 50.2 = 4.5 mmol m −3 with   import anth taken as C anth for σ 0 = 23.0-24.0 in the northwestern subtropical region and with   export anth taken as the mean C anth for σ 0 = 25.0-25.5 in the northwestern and southwestern subtropical region (Table 1).The net input of 4.8 TgC year −1 from the simplified model is comparable to the value (62.9 − 58.1 = 4.8 TgC year −1 ) deduced from the WMT analysis (Table 1), in which the input (62.9) is estimated as 80.6% of the total production due to conversion (78.1 TgC year −1 ).
Similarly for the low-resolution model, the upward volume transport is estimated to be less than the high-resolution model by 2.0 Sv based on Figure 5d, which gives I stmw ∼ 0.8 Sv.With   import anth −  export anth ∼ 5.4 mmol m −3 (=54.0 − 48.6) estimated in the same way as the high-resolution model (Table 1), a net input of 1.7 TgC year −1 is obtained for the first term on the right-hand-side of Equation 17.This value is smaller than that deduced from the WMT analysis (3.4 TgC year −1 ; Table 1); however, it should be noted that the upward export transport (9.2 TgC year −1 ; Figure 14d) from the STMW recirculation in the WMT analysis has been possibly underestimated because this has been offset by the downward diapycnal transport from above which is not necessarily within the STMW recirculation (e.g., the downward flux in the East China Sea to the west of the Kuroshio seen in Figure 15a).A larger export would result in a smaller storage rate for the STMW recirculation in the above quantitative estimation deduced from the WMT analysis.However, the C anth storage rate of the STMW recirculation in the low-resolution model could be actually smaller than the sum of those in the southwestern and northwestern subtropical regions since STMW occupies relatively small area in the σ 0 = 25.0-25.5 layer (Figures 9a and 9d).Thus, the simple budget model may be regarded to represent the STMW recirculation of the low-resolution model well, too.
The implication from the simplified model Equation 17 would be as follows.In addition to the incorporation of the larger fraction (γ) of surface uptake into the STMW recirculation in the high-resolution model due to the structural change in the STMW circulation, the presence of a stronger overturning circulation (I stmw ) would increase the storage rate of STMW if the change of ) due to increased horizontal resolution is small.Therefore, it would be necessary to consider next why the tendency of   import anth is kept almost comparable to that of the low-resolution model despite the change of circulation in the high-resolution model.However, before considering the surface layer budget, the importance of the structural change in the STMW formation and circulation is emphasized to keep the level of   export anth in the high-resolution model, where the STMW volume is also increased.In the high-resolution model (Figures 14b and 14d), 255.1 TgC year −1 of C anth is transported downward across the surface of σ 0 = 25.0, and 62.9 TgC year −1 is incorporated in the STMW recirculation, which is about 25%.In the low-resolution model (Figures 14a and 14c), 171.2 TgC year −1 of C anth is transported downward, and 12.6 TgC year −1 is incorporated in the STMW recirculation, which is about 7%.The separation of the Kuroshio off the eastern coast of Japan and the proper representation of the Kuroshio Extension and the intense recirculation gyre are critically important for the large amount of C anth to be isolated and used to increase   export anth .

C anth Budget for the Surface Layer
In the surface layer from the sea surface to σ 0 = 25.0 (Figures 14a and 14b), the C anth convergence due to diapycnal transport from below in the southwestern subtropical region, the northward transport at 28°N, and the divergence due to downward diapycnal transport in the northwestern subtropical region are all enhanced by at least  This might be an underestimation due to the presence of the gain due to the surface buoyancy loss from above (cf.i).g Surface uptake in the SW subtropical region (0.5) is subtracted from the sum of storage rates in the SW and NW subtropical regions.h Surface uptake in the SW subtropical region (0.6) is subtracted from the sum of storage rates in the SW and NW subtropical regions.i This might be an overestimation because the STMW recirculation occupies a small area in the NW and SW subtropical regions.(cf.f).

Table 1
Comparison Between a Simple Tracer Budget Model (Equation 17) and a Quantitative Estimation Deduced From the Water Mass Transformation (WMT) Analysis (Figures 14c and 14d The C anth budget for the surface layer of the southwestern subtropical region implies that the difference in the northward transport of C anth (269.4 − 190.0 = 79.4TgC year −1 ) at 28°N cannot be entirely explained by that of the transport from below (50.5 − (−7.1) = 57.6TgC year −1 ).Since almost no difference is seen in the surface uptake, the excessive supply from the eastern boundary seems to be necessary to maintain C anth in the southwestern subtropical region at a high level.The excessive westward transport from the eastern to the southwestern subtropical regions (182.5 − 156.3 = 26.2TgC year −1 ) at 170°W is achieved at the expense of the smaller downward transport to the higher density classes σ 0 > 25.0 such as the eastern STMW (Hautala & Roemmich, 1998) in the high-resolution model, which can be confirmed using Figures 15a and 15b.

How Should the Column Inventory Be Understood?
The difference of the column inventory between the low-and high-resolution models shown in Figure 6c would be explained as follows.In the high-resolution model, more C anth is accumulated and isolated in the STMW recirculation in the northwestern subtropical region, increasing the column inventory.The increased C anth is supplied from the surface water in the southern part of the subtropical gyre, at the expense of a less amount of formation of the eastern STMW in the eastern subtropical region.The increased volume of STMW in the high-resolution model occupies a vast area of the western subtropical gyre, pushing aside the upper and lower layer waters.The inventory and storage rate of C anth in the lower-and higher-density ranges indeed become smaller in the southwestern subtropical region (Figures 8e and 12a).The Kuroshio displaces the reduced part of the upper layer water northward to increase STMW through the formation in the northwestern subtropical region.In the lower layer, CMW is displaced toward the outskirt of the subtropical gyre, reducing the inventory below STMW.As a result, the difference in the column inventory is nearly zero in the southwestern subtropical region (Figures 6c and 8c).
The deficit of C anth in the eastern STMW due to the smaller downward transport of C anth from the surface layer in the high-resolution model is compensated mainly by the shift of the newly formed CMW to the outskirt of the subtropical gyre and partly by the southward Ekman transport of the excess C anth uptake in the subpolar region in the high-resolution model (Figures 6k and 13c).As a result, almost no difference is seen in the column inventory in the northern part of the eastern subtropical region (Figure 6c).The net increase of the southward transport through the section of 42°N-Soya Strait (4.9 TgC year −1 ) is almost equivalent to the increase of the column inventory in the northwestern subtropical region (4.7 TgC year −1 ; Figure 13), which fills the loss of C anth column inventory due to the surface transport from the eastern subtropical region to the southwestern subtropical region.
Figure 17 summarizes the change of the column inventory of C anth by increasing the horizontal resolution.Previous studies (e.g., Iudicone et al., 2016;Nakano et al., 2015) have proposed the concept that the surface layer part of the shallow MOC in the tropical and subtropical ocean is preconditioning the high C anth inventory in mode waters by transporting surface C anth -rich water from the low latitudes to the northern part of the subtropical region.The present result is consistent with this concept and adds a new finding that the high-resolution model realizes a distinct formation for STMW and CMW in the North Pacific Ocean, with an enhanced formation of STMW accompanied by an enhanced MOC.This increases C anth inventory in STMW, making the overall distribution of C anth inventory closer to observation, specifically in the northwestern subtropical region.The change in the distribution is explained as the redistribution of surface and subsurface C anth within the subtropical gyre and not as the change in the amount of surface C anth uptake.

Summary and Conclusion
In this study, simulations by a global ocean biogeochemical model with and without embedding a high-resolution model in the North Pacific Ocean are compared to investigate the impact of resolving boundary currents, oceanic fronts, and mesoscale eddies on the carbon circulation in the North Pacific Ocean.Considering practical applications to climate projections and possibly to ongoing assessments of carbon budgets (e.g., DeVries et al., 2023;Friedlingstein et al., 2022;Hauck et al., 2020), the embedded model is run for the recent 80 years after branching off from a low-resolution global model that has been spun-up without embedding.The physical and biogeochemical fields for the high-resolution model adjust in the first 20 years, while the last 40 years are compared with the low-resolution counterpart.The results suggest that the basic physical and biogeochemical fields may be improved in high horizontal resolution models, especially for the region where boundary currents, fronts, and eddies are involved, as expected (e.g., Hewitt et al., 2020;Kiss et al., 2020).An improved surface CO 2 flux is also found in the subpolar region of the North Pacific due to the higher surface DIC concentration realized by upwelling of the high DIC water from the subsurface.Although ecosystem processes would be also important for the balance of surface DIC in high-latitudes (e.g., Takahashi et al., 2002), the higher surface DIC caused by the improved physical processes would be a robust step toward the overall improvement.The present result is consistent with the concept that C anth -rich water from low latitudes is preconditioning the high C anth inventory in STMW (e.g., Iudicone et al., 2016;Nakano et al., 2015) and adds a new finding that this holds for the enhanced component of the STMW formation in the high-resolution model that accompanies the eddy-induced, enhanced MOC.It is emphasized that the WMT framework (Groeskamp, Griffies, et al., 2019) is very useful in drawing these conclusions based on quantitative analysis.Specifically, an estimate on the C anth budget for the STMW recirculation was established, which is one of the major regions of C anth storage in the world ocean (e.g., Carter et al., 2019;Gruber et al., 2019a), using a simulation result of a high-resolution model where the Kuroshio, its extension, recirculation, and STMW are represented realistically, for the first time in this research field.The improvements due to taking a high horizontal resolution for an ocean biogeochemical model presented here would motivate climate modelers to use high-resolution ocean biogeochemical models in the next-generation climate projections.The projection of acidification rate of subsurface waters in the western subtropical region as assessed by Resplandy et al. (2013) and Watanabe and Kawamiya (2017), would be reinforced by such a model.The change in the horizontal resolution may impact projections of carbon cycles when the difference between the distributions of the C anth storage in low-and high-resolution models is more pronounced in the future.This study has also demonstrated the feasibility of embedding regional ocean biogeochemical models into a global ocean model for use in climate projection experiments.As long as the use of globally 10-km class ocean biogeochemical models is not feasible, an embedding method could be another possible option for modelers who are rather concerned about the regional change of ocean biogeochemistry, especially where boundary currents and mesoscale phenomena are involved.
As a final note, a global ocean carbon cycle simulation requires a long-term spin-up before the ocean carbon cycle reaches a quasi-steady state.Some deviations from the observed distributions, at least locally, for some biogeochemical properties are unavoidable during the spin-up even for the state-of-the-art models (e.g., Stock et al., 2020;Yool et al., 2020).Both the physical and ecosystem processes can be improved to reduce these deviations.Considering the possible merits of high horizontal resolution as demonstrated in this study, the improvement of ecosystem processes should be preferably sought within a simple model framework to limit the overall computational cost (e.g., Dunne et al., 2020aDunne et al., , 2020b)).This should be a subject of future development efforts.
Let us consider to transform the above set of equations in depth coordinates (x, y, z, t) to density coordinates , where tildes denote that the calculation of partial derivatives is performed along sloping isopycnals.By using the following transformation rule, ) we obtain, Note that in the above equations, and ) can be considered represent thickness of density layer.However, since we take z as positive upward and z ρ is generally negative, we should regard −z ρ as thickness.The velocity in the advection term in Equation A13 can be considered as transport because it is multiplied by thickness z ρ .The property ω defined in Equation A14 (Units: kg m −3 • s −1 ) is related to the diapycnal transport across density surfaces e ≡ z ρ ω (Units: m s −1 ), which is important in considering MOC and WMT.From Equations A8 and A14, we can obtain diapycnal flux ω as We regard Equation A17 as implying that dissipation of tracer in an isopycnal layer bounded by ρ u (z) < ρ < ρ l (z) (the left-hand-side) can be evaluated using the depth integral of tracer dissipation term based on depth coordinate (the right-hand-side), This indicates that the sampling of tracer dissipation terms at run time can be used to evaluate dissipation in a density layer.The evaluation can be done without leak or source.Transformation of fluxes in depth coordinates to those in density coordinates as in Equation A16 is not necessary.

A2. Specific Treatments for the Present Analysis
Next, we consider to apply the above formulation to the tracer budget analysis for our model results, taking into account of the specific settings and schemes used for the model calculation.Readers are referred to Groeskamp, Griffies, et al. (2019) for general discussions about the tracer budget analysis on density coordinates.

A2.1. Equation for Density Evolution and Diapycnal Transport
When Gent and McWilliams (1990) parameterization for eddy induced transport is expressed as skew flux in the diffusion tensor as introduced by Griffies (1998), the evolution of density on depth coordinates is expressed as follows: where the first term on the right hand side is due to convergence of Gent-McWilliams skew flux, the second term is isopycnal diffusion, the third term is vertical diffusion, and the fourth term is convergence of buoyancy flux due to external forcing.Note that the advection term on the left hand side is due to the resolved velocity field, which does not contain parameterization for eddy-induced transport.By transforming Equation A20 to density coordinates, we obtain diapycnal transport: where κ gm is coefficient for this parameterization (Units: m 2 s −1 ).
Comparing Equations A20 and A22, we have, Thus, we can evaluate all relevant terms of diapycnal transport: z ρ ω from Equation A21, z ρ ω tot from Equation A32, and z ρ ω gm from Equation A33.

A2.2. Equation for Tracers
The evolution of tracer on depth coordinates is expressed using skew flux as follows: (A39) When the skew flux is used to express the Gent and McWliiams parameterization for eddy-induced velocity, the component of skew flux is not the advection flux due to eddy induced transport.This means that it is not necessarily appropriate to use Equation A39 for the purpose of tracer budget analysis based on the data sampled during runtime.In other words, it is very difficult to estimate advection flux due to eddy induced velocity.Also, it is very difficult to accurately compute ω from u and v while it is relatively easy to compute ω tot from sampled mixing terms of temperature and salinity.Based on this consideration, and using ( gm ( gm from Equations A35 and A38, we write Equation A39 as We can compute tracer budgets based on runtime sampling of each term in Equation A41.

A2.3. Tracer Budget
For considering tracer budget for a density range over a fixed region based on Equation A41, we should extract horizontal fluxes from the flux divergence terms to evaluate fluxes through lateral boundaries.Relevant terms are the first and second advection terms and fifth (skew flux) and sixth (isopycnal mixing flux) terms on the right hand side.For the fifth and sixth terms, we should evaluate the sum of the interior divergences and the boundary fluxes separately.Integrating Equation A41 for a density range and over a fixed horizontal region, whose infinitesimal element is   Ã (Units: m 2 ), we obtain budgets of the storage rate for the tracer (Units: (units of tracer) • m 3 • s −1 ), where and n is a unit normal vector that orients outward from the region boundary and ds is a line element that orients along the boundary (Units: m).The conversion term represents ) . (A45) Figure 1.

Figure 1 .
Figure1.Time series of globally integrated, annual mean surface uptake of CO 2 (PgC year −1 ) (a) during the long-term spin-up (189-1,652) and the last 366 years (1653-2018) with xCO 2 that starts to increase from year 1850 for the global model without embedding and (b) after the embedding of the high-resolution North Pacific Ocean model from year 1938.Details about the experiment are described in Section 2.2.In (a), surface uptake during the spin-up and the last 366 years for the natural dissolved inorganic carbon (DIC) tracer that sees the constant xCO 2 of 284.32 ppm is depicted with the blue line.The one for the contemporary DIC that sees an increasing xCO 2 is depicted with the red line (also in (b)).In (b), surface uptake for the embedded run is depicted with the darkblue line.Means of observationbased products (green) and ocean biogeochemical models (yellow) provided by Global Carbon Project (2022) are also depicted with uncertainty bounds represented by ±1 standard deviation in shades.(c-f) Long-term mean, annual mean surface uptake of CO 2 (mol m −2 year −1 ) for the global model (c) without and (d) with embedding, (e) difference between with and without embedding, and (f) an observation-based estimate (JMA-MLR) byIida et al. (2021).In (d) and (e), the embedded region is represented by the purple rectangle.For (c) and (d), root-mean-square error and mean bias north of 12°N relative to JMA-MLR are put on top of panels.The dashed rectangle in (c)-(f) is used for a quantitative assessment of simulations (see text in Section 3).
The embedding of the high-resolution North Pacific model within the global model started at the beginning of the year 1937.The high-resolution model started from a state of rest with a horizontally uniform vertical distribution of temperature and salinity, which was calculated by horizontally averaging those of the global model, and horizontally interpolated fields from the global model for other tracers (biogeochemical and inert ones) at the starting time of the high-resolution model.In January 1937, tracer fields were nudged very strongly with a 0.2-day timescale toward those temporally interpolated from the time series of monthly means of the global model.By the beginning of February 1937, the tracer fields of the high-resolution model became almost those of the global model, and the flow fields were established by the geostrophic adjustment and the surface forcing.The nudging for tracer fields was stopped on 1 February 1937, and the set of the embedded high-resolution North Pacific and the outer global models was forced only at the surface and integrated until the end of year 2018.Exchange of variables at the boundary between the models was one-way (only from the global to the North Pacific model) for the first 2 months (January and February 1937) and two-way after March 1937.
).The southern MLD maximum formed to the south of the Kuroshio Extension in the high-resolution model is related to the formation of STMW.The differences of surface heat flux almost exactly reflect those of SST (Figures2m-2o).The larger heat release in the midlatitude western to central North Pacific in the high-resolution model implies enhanced formation of mode waters due to density change of surface water masses.Looking at other regions, the positive bias of SSH in the subpolar region of the low-resolution model is reduced in the high-resolution model, implying that the subpolar gyre deepens in the high-resolution model.Alternatively, this would imply that the flattening effect for slanting isopycnal surfaces by the Gent-McWilliams parameterization is excessive in this region of the low-resolution model.The low SST bias immediately north of 10.1029/2023MS003720 9 of 38 the Equator in the low-resolution model is reduced in the high-resolution model presumably due to the proper representation of tropical instability waves in the high-resolution model.

Figure 2 .
Figure 2. Evaluation of the physical properties in terms of long-term mean.Sea surface temperature (SST) for 1989-2018 (first row: a-d), sea surface height (SSH) for 1993-2018 (second row: e-h), mixed layer depth (MLD) in winter (mean of February and March) for 1989-2018 (third row: i-l), and surface heat flux (Q tot ) for 1988-2017 (bottom row: m-p).Note that the periods for averaging are adjusted according to the available period of the reference observation-based data.In the first and second column, bias relative to the observation-based estimates (SST, SSH, and Q tot ) or long-term mean (MLD) is depicted for the global model without and with embedding, that is, the low-and high-resolution model, respectively, for the North Pacific Ocean.On top of panels, root-mean-square error and mean bias north of 12°N relative to the observation-based estimates are put (Sea of Japan is excluded for MLD).For the third column, difference between the high-and low-resolution models is depicted.For the fourth column, the reference data based on observations is depicted.SST from Program for Climate Model Diagnosis and Intercomparison (PCMDI) (Durack & Taylor, 2019) following a procedure described by Hurrell et al. (2008), SSH from the EU Copernicus Marine Environment Monitoring Service (CMEMS), MLD from de Boyer Montégut et al. (2004), and Q tot from J-OFURO3(Tomita et al., 2019).

Figure 3 .
Figure3.Surface partial pressure of CO 2 (pCO 2 , units in μatm) for 1990-2018 (first row: a-d) and vertically integrated net primary production (intPP, units in molC m −2 year −1 ) for 2003-2018 (second row: e-h).In the first and second column, bias relative to the observation-based estimate (pCO 2 ) or long-term mean (intPP) is depicted for the low-and high-resolution model, respectively.For pCO 2 , root-mean-square error and mean bias north of 12°N relative to the observation-based estimate are put on top of panels.For the third column, difference between the high-and low-resolution models is depicted.The dashed lines in (c) and (d) indicate the longitude of 164.5°E along which meridional sections are drawn in Figure4for biogeochemical variables.For the fourth column, the reference data based on observations is depicted.Surface pCO 2 from JMA-MLR(Iida et al., 2021) and intPP from a mean of three MODIS-based products, the original and an Epply version of the Vertically Generalized Production Model (VGPM)(Behrenfeld & Falkowski, 1997) and Carbon-based Production Model (CbPM)(Westberry et al., 2008) provided by Oregon State University.

Figure 4 .
Figure 4. Evaluation of the biogeochemical variables in the western North Pacific Ocean along 164.5°E.Long-term mean for 1987-2017 is evaluated.Dissolved inorganic carbon (DIC, mol m −3 ) in the first row (a-d), total alkalinity (Alk, mol m −3 ) in the second row (e-h), and nitrate (NO 3 , μmol L −1 ) in the third row (i-l).In the first and second column, bias relative to the observation-based estimates is depicted for the low-and high-resolution model, respectively.On top of panels, root-mean-square error and mean bias north of 12°N relative to the observation-based estimates are put.In the third column, difference between the high-and low-resolution models is depicted.In the fourth column, the reference data based on observations is depicted.DIC and Alk from GLODAPv2 and NO 3 from WOA13v2.

Figure 5 .
Figure5.Long-term (1987Long-term ( -2017) )  mean Meridional Overturning Circulation (MOC) for the Indo-Pacific Ocean.Units in Sverdrups (1 Sverdrup = 10 6 m 3 s −1 ).Note that the total velocity, that is, resolved plus parameterized eddy induced velocity (EquationA24), is integrated zonally over the basin.(a) MOC on depth-coordinate of the high-resolution model and (b) its difference from the low-resolution model.(c, d) Same as (a, b) but for MOC on density coordinate (σ 0 ) lighter than σ 0 = 27.5 that corresponds to the depth range shallower than approximately 1,000 m.See FigureS1in Supporting Information S1 for MOC on depth-coordinate in other basins.Note that regions with negative values imply anti-clockwise circulations and vice versa.

Figure 6 .
Figure6.The column inventory in mol m −2 (first row: a-d), storage rate in mol m −2 year −1 (second row: e-h), and surface uptake rate in mol m −2 year −1 (bottom row: i-k) of C anth averaged for 1987-2017.In the first and second column, long-term mean is depicted for the low-and high-resolution model, respectively.For the column inventory, bias relative to the observation-based estimate (GLODAPv2) is depicted with root-mean-square error and mean bias north of 12°N put on top of panels.In the third column, difference between the high-and low-resolution models is depicted.In the fourth column, the reference data based on observations is depicted for the column inventory in 2002 (GLODAPv2) and the storage rate as represented by the rate of change between 1994 and 2007(Gruber et al., 2019a).In the bottom row (l), zonally integrated surface uptake (light blue and pink for the low-and high-resolution models, respectively) and storage (blue and red for the low-and high-resolution models, respectively) rates, both in units of TgC year −1 , in the Pacific Ocean averaged for 1987-2017 are depicted.The dashed rectangle in (a)-(d) is used for a quantitative assessment of simulations (see text in Section 4.1).

Figure 7 .
Figure 7. C anth concentration bias relative to the observation-based estimate (GLODAPv2) (mmol m −3 ; shade) and density (σ 0 ; contour) in the upper western North Pacific Ocean along 164.5°E for (a) low-and (b) high-resolution model.On top of panels, root-mean-square error and mean bias north of 12°N relative to the observation-based estimate are put.(c) Difference of C anth concentration between the high and low-resolution models (shade).Blue lines are density distribution of the low-resolution model depicted in (a) and red lines are that of the high-resolution model depicted in (b).Long-term mean for 1987-2017 is shown.(d) The C anth concentration from GLODAPv2 and density from WOA13v2.The dashed rectangle in (a)-(c) denotes an approximate region of existence of the subtropical mode water (STMW) in the high-resolution model and observation-based data.

Figure 8 .
Figure 8. Mean and trend of properties related to C anth for the western subtropical region.The domain of the high-resolution North Pacific model is divided into regions depicted in (a).Monthly time series of the column inventory of C anth in PgC (b, c), the 31-year (1987-2017) mean storage of C anth in PgC partitioned in density classes (d, e), monthly time series of the C anth storage in PgC (f, g), volume in km 3 (h, i), C anth concentration in mmol m −3 (j, k) in the subtropical mode water layer (σ 0 = 25.0-25.5),for the northwestern and southwestern subtropical region in the left and right panels, respectively.Blue and orange are used for the low-and high-resolution models, respectively, and red for the difference.

Figure 9 .
Figure 9. Horizontal distribution of C anth storage in mol m −2 (upper panels) and thickness in meter (lower panels) for the density layer of subtropical mode water (σ 0 = 25.0-25.5).In the first and second column, the long-term (1987-2017) means in shade and the Montgomery potential (units: Pa) with contours (only for upper panels) are depicted for the low-and high-resolution model, respectively.Contour interval for the Montgomery Potential is 1,000 Pa.For the third column, difference between the high-and low-resolution models is depicted.

Figure 10 .
Figure 10.Same as Figure 9 but for the density layer below subtropical mode water (σ 0 = 26.0-26.5).Contour interval for the Montgomery Potential is 500 Pa.

Figure 11 .
Figure11.Budget for the C anth storage rate separated in density classes in the northwestern subtropical region for 1987-2017 (Units in TgC year −1 ).In the first row, (a) storage, (b) surface uptake, (c) transport convergence, and (d) conversion rates.In the second row, transport through (e) the western part of 28°N, (f) the western part of 42°N, (g) the northern part of 170°W, and (h) the sum of Tsushima and Tsugaru Strait.For the top two rows, vertically integrated values for both models are put on top of each panel except for (d).In the third row, the conversion rate shown in (d) is split into components as formulated by Equation 3. Convergence of diapycnal transport of C anth due to (i) explicit mixing (ω tot ) and (j) the Gent-McWilliams scheme (−ω gm ).Conversion due to (k) vertical diffusion and (l) sum of Gent-McWilliams skew flux and isopycnal diffusion terms for C anth .In the bottom panels, the part unexplained by the explicit mixing terms shown in the third row is depicted on the lefthand side (m).For the rest of the row, convergence of diapycnal transport of C anth due to explicit mixing shown in (i) is further split into components as formulated by Equation 5: (n) vertical diffusion, (o) surface buoyancy flux, and (p) isopycnal diffusion terms.Blue and orange are used for the low-and high-resolution models, respectively.Note the difference in the vertical scales used in panels.

Figure 12 .
Figure 12.Same as Figure 11 but for the southwestern subtropical region.For the second row, transports through (e) the western part of the 12°N, (f) the western part of the 28°N, (g) the southern part of 170°W are depicted.Alphabetical title for each panel is shifted forward by one in the lower two rows.

Figure 13 .
Figure 13.Budgets of the column integrated C anth inventory rate (TgC year −1 ) for the separated regions depicted in Figure 8a for the period 1987-2017 evaluated according to Equation 6.(a) Low-and (b) high-resolution model, and (c) their difference (high minus low).Inventory rate in pink, surface flux in light-blue, and horizontal transport in yellow.On the right side, budgets for tropical, subtropical, and subpolar regions separated zonally by 12°S-Indonesian Through Flow, 12°N, 42°N-Soya-Tartar Strait, and Bering Strait sections are depicted.Distributions in density classes in the northwestern and southwestern subtropical regions are shown in Figures 11 and 12, respectively.Those in the eastern subtropical region are shown in Figure S4 in Supporting Information S1.Volume transports through the boundaries of the separated regions are shown in Figures S6 and S7 in Supporting Information S1 for the low-and high-resolution models, respectively.In the same way, C anth transports are shown in Figures S8 and S9 in Supporting Information S1.

Figure 14 .
Figure 14.Budget of the C anth storage rate (TgC year −1 ) for (a, b) the surface layer (from sea surface to σ 0 = 25.0) and (c, d) subtropical mode water layer (σ 0 = 25.0-25.5)for the period 1987-2017 evaluated according to Equation 1. Storage rate in pink, surface flux in light-blue, horizontal transport in yellow, and convergence of diapycnal transports across isopycnal surfaces bounding the layer in dark-blue.Left panels (a, c) are for the low-resolution model and right panels (b, d) are for the high-resolution model.For the sections along 170°W in (c) and (d), the transport is separated into the recirculation component (dark-yellow) and others (light-green).The dark red arrow connects the recirculation transport from the northern to southern section of 170°W and the dark blue arrow connects the southwestward transport from the southern 170°W section to the western 12°N section.See the main text in Section 5 for details about drawing the arrows.

Figure 15 .
Figure 15.Horizontal distribution for some quantities relevant to C anth budget for the density layer of σ 0 = 25.0-25.5.Diapycnal transport of C anth (a, b) for the upper surface (σ 0 = 25.0) and (c, d) for the lower surface (σ 0 = 25.5),(e, f) convergence of the diapycnal transport of C anth for this layer, and (g, h) surface uptake of C anth by this layer.Units in mol m −2 year −1 .The left and right columns for the low-and high-resolution model, respectively.Contours are Montgomery potential with 1,000 Pa interval.

TSUJINO
) for C anth Budget of the Subtropical Mode Water (STMW) Recirculation in the Density Range σ 0 = 25.0-25.5 year −1 in the high-resolution model, consistent with the above evaluation about the STMW recirculation.

A
larger C anth amount is stored in the northwestern subtropical gyre in the high-resolution model, which compares better with observation-based estimates.The improved North Pacific STMW that occupies a vast region south of the Kuroshio and its extension in the high-resolution model (e.g.,Nishikawa et al., 2010) contributes to the increased storage of C anth .The larger storage rate in STMW of the high-resolution model is due to two reasons: The enhanced formation of new water mass along the Kuroshio and the larger fraction of the new water mass being incorporated into STMW in the Kuroshio recirculation and isolated from outside.Source of C anth for the increased part of the new water mass is the C anth -rich surface water transported from the southern part of the subtropical gyre by the intensified Kuroshio.It is transformed into STMW through large buoyancy loss.Somewhat contrary to the expectation, the surface uptake of C anth in the region of large buoyancy loss is slightly increased but is largely comparable to that of the low-resolution model.Regarding the column inventory in the eastern subtropical region, the loss of the surface C anth storage from the eastern to northwestern subtropical region is compensated by the supply of the recently ventilated denser waters in the lower layers such as CMW and partly by the surface southward transport of the enhanced C anth uptake from the subpolar gyre, keeping the column inventory almost unchanged there from the low-resolution model.Thus, the change in the distribution of C anth storage in the high-resolution model can be understood as a redistribution of newly gained C anth in the subtropical gyre by the improved circulation.

Figure 17 .
Figure 17.Schematic explaining redistribution of C anth within the subtropical region due to change in the horizontal resolution of the model.Background shading is the difference of the mean (1987-2017) column inventory of C anth between the high-and low-resolution models shown in Figure 6c.

=
skew + iso + vdiff + surf.(A21)Thiswould constitute diapycnal velocity in the continuity Equation A11.On the other hand, when the Gent-McWilliams parameterization for eddy induced transport is expressed as velocity and introduced in the advection term as part of total velocity, Equation A20 takes the form + tot + tot + tot, (A23) and the total velocity is defined astot =  + gm, tot =  + gm, tot =  + gm.(A24)Velocity components for the eddy induced transport scheme by the Gent-McWilliams parameterization are given as follows: represents the zonal and meridional component of velocity, respectively,  F skew is the Gent-McWilliams skew flux,  F iso is flux due to isopycnal diffusion, n is a unit normal vector that orients outward from the region boundary and ds is a line element that orients along the boundary.In the last term, Surface uptake (21.1) in the NW subtropical region multiplied by the fraction (0.806) into STMW recirculation.b Estimated as a simple mean of C anth for σ 0 = 25.0-25.5 in the NW and SW subtropical regions.c Surface uptake (19.7) in the NW subtropical region multiplied by the fraction (0.562) into STMW recirculation.d Total conversion (78.1) in the NW subtropical region multiplied by the fraction (0.806) into STMW recirculation.e Total conversion (22.4) in the NW subtropical region multiplied by the fraction (0.562) into STMW recirculation.f a =  T   skew +  T   iso +  T   vdiff +  T   surf .Since the advection due to eddy induced velocity is expressed in terms of skew flux,gm + gm + gm = − T     iso +  T   vdiff +  T   surf .