Contribution of Biological Effects to Carbonate-System Variations and the Air–Water CO 2 Flux in Urbanized Bays in Japan

We evaluated the contribution of net biological effects (photosynthesis + respiration + decomposition) to the carbonate parameters and air–water CO 2 fluxes in Tokyo Bay, Ise Bay, and Osaka Bay in Japan. The carbonate parameters (fugacity of CO 2 , total alkalinity, and dissolved inorganic carbon) were measured mainly by cargo ships traveling between Japan and other countries. We used the measurement data from three inner bays and surrounding outer bays in Japan along with reference data from previous studies for complementary analysis. We found that (a) the inner and outer bays in this study were strong annual atmospheric CO 2 sinks, (b) the annual biological effect on the air–water CO 2 fluxes was about 6%–27% of the measured CO 2 fluxes, and it affected the seasonal variation of the CO 2 flux, the biological effect was largest in in Ise and Osaka intensity of the biological effect corresponded with nutrient concentrations, which seemed to be controlled by the wastewater treatment in urbanized areas around the bays. Our results suggest that labile carbon/ nutrient ratio should be a major consideration for evaluating the biological effect on the carbon cycle in urbanized bays, which will continue to expand globally.

have reported that urbanized inner bays in Japan are annual atmospheric CO 2 sinks (Endo et al., 2018;Fujii et al., 2013;Kubo et al., 2017). Especially in Tokyo Bay, CO 2 undersaturation was indicated to be the result of wastewater treatment by continuous measurement of carbonate system (Kubo & Kanda, 2020;Kubo et al., 2017;Kuwae et al., 2016). They suggested that the treatment process removes labile (easily decomposed to inorganic carbon like CO 2 ) carbon, yielding water with relatively less carbon than nutrients, which promotes primary production in the inner bay water. In addition, the organic matter remaining in the treated water is refractory (hardly decomposed to inorganic carbon) (Kubo et al., 2015), and thus, further mineralization and increase in the CO 2 concertation in the inner bay water are suppressed. However, the measurement area in these previous studies was limited within the inner bays; thus, the effect of the biological processes such as photosynthesis, respiration, and decomposition could not be distinguished from the effects of non-biological processes due the mixing of water outside of the inner bay (riverine and oceanic water). Therefore, the biological effect on the carbonate parameters and the air-water CO 2 flux has not been quantified precisely.
Here, we evaluated the biological effects of inner and outer bay water on temporal and spatial changes in the carbonate parameters and the air-water CO 2 flux by quantifying the non-biological effects of outside water in three bays in Japan. We discuss the biological effects on the flux in terms of the regulating factors and the extendibility or our observations to other areas and for global estimation.

Study Site
This study took place in Japan in Tokyo Bay (are: 1,400 km 2 ), Ise Bay (1,700 km 2 ), and Osaka Bay (1,300 km 2 ), and surrounding areas ( Figure 1); in this study, the term "bay" refers to the inner and outer areas of each bay, which were analyzed together for expedience. These bays have similar topographic conditions such as a southward entrance to the bay open to the Pacific Ocean (Kuroshio area) and a surrounding, highly urbanized coastal area. The populations of the urban areas surrounding Tokyo, Ise, and Osaka Bays were 37, 9, and 19 million, respectively, in 2014 (Global Metro Monitor; https://www.brookings.edu/research/ global-metro-monitor/). The largest rivers are Edo River (discharge: ∼100 m 3 s −1 ), Kiso River (∼200 m 3 s −1 ), and Yodo River (∼200 m 3 s −1 ) in Tokyo Bay, Ise Bay, and Osaka Bay, respectively; the discharge of the rivers is larger than the sum of other three largest rivers in each bay (Japan River Association; http://www.japanriver.or.jp/river_law/map.htm). The wastewater from the urban area into the inner bays was treated mainly by activated sludge method for decreasing the organic pollutant. There are 10 or more sewage treatment plants around the inner bays and the location of the plant with the largest amount of treatment is show in Figure 1.

Data Processing
The data used for this study are measurements of water temperature, salinity, and fugacity of CO 2 (fCO 2 ) in the main urbanized inner bays and in the surrounding outer bays, as observed by the National Institute for Environmental Studies (NIES), Tokyo University of Marine Science and Technology (TUMSAT) and Osaka City University (OCU). The other carbonate parameters-total alkalinity (TA) and dissolved inorganic carbon (DIC)-were estimated using an empirical relationship between TA and salinity, and the equilibrium calculation. The biological effect on carbonate parameters and the air-water CO 2 flux was calculated from the difference between the above estimated DIC and the value from the conservative mixing line between the oceanic and riverine endmembers. We quantified the air-water CO 2 flux and the magnitude of the biological effect to evaluate the contribution of the biological effect in the bays to the exchange of atmospheric CO 2 .
The water temperature, salinity, and fCO 2 in water and air were obtained from the NIES database (http:// soop.jp). The data are also available from the Surface Ocean CO 2 Atlas (SOCAT; https://www.socat.info), which has been a public database since 2011 and represents an international collaboration among research institutes. The NIES observations implemented as the Voluntary Observing Ship (VOS) programs by cargo ships were originally for understanding the global carbon cycle, but their data also include inner and outer bay measurements from 2005 to 2016 in Tokyo Bay and Ise Bay and from 2011 to 2016 in Osaka Bay, where the cargo ships taking measurements have anchor stations. All parameters were taken by the underway measurement during the cargo shipping. The water temperature was measured by a sensor attached at the depth of 7 m below the water surface. The salinity and fCO 2 in water were measured from the water sample, which was taken at 7 m depth. The fCO 2 in air was measured by a sensor attached at 30 m height above the water surface. The details on the measurements are shown in Table 1 and reported in Nakaoka et al. (2013) and the Metadata of the cargo ships (M/S Pyxis, M/S New Century 2, M/S Trans Future 5: https://www.socat. info/index.php/data-access/). We first extracted the data for 30-40°N and 130-145°E observed by the cargo ships as the original coastal data. TOKORO ET AL.  Tokyo Bay, (c) Ise Bay, (d) Osaka Bay. The black filled circles indicate the locations of the National Institute for Environmental Studies (NIES) data used in this study. The blue lines and dots indicate additional data collected by Tokyo University of Marine Science and Technology (TUMSAT) in Tokyo Bay, which was the dataset by TR/V Seiyomaru in Kubo et al. (2017), and by Osaka City University (OCU) in Osaka Bay which was the same dataset used in Endo et al. (2018). The red stars mark the the largest river mouths located in the inner part of each bay and used as the zero point for the distance parameter dist (Tokyo Bay,35.68°N,139.94 E;Ise Bay,35.04°N,136.74°E;Osaka Bay,34.68°N,135.41°E). The green diamonds are the sewage treatment plant with the largest amount of treatment in each bay (Morigasaki; 14.9 m 3 s −1 in Tokyo Bay, Uchide; 2.4 m 3 s −1 in Ise Bay, Toba; 8.6 m 3 s −1 in Osaka Bay, Japan Sewage Works Association, 2010).
Although the NIES data have the information including the surrounding outer bay, the landward information is low-resolution due to the difficulty of near-shore measurement by cargo ships. For the detailed analysis of the riverine effect, we included some previous data as more landward information than the NIES data. For Tokyo Bay, we added the data collected by TUMSAT from 2007 to 2010 to the NIES data. The data are almost the same as that used in the previous study (Kubo et al., 2017), but a part of data (measured by TR/V Hirodori) was not included due to the different data format from other larger data (measured by TR/V Seiyo-maru). In addition, the pCO 2 data were converted to fCO 2 using the empirical relationship incorporating temperature (Kӧrtzinger, 1999). Likewise, we added more landward measurement data from Osaka Bay collected by OCU in spring and autumn of 2014 (the same dataset used for Osaka Bay in Endo et al., 2018). The additional raw data for Osaka Bay included the water temperature, salinity, and DIC. We therefore estimated fCO 2 using the equilibrium calculation (Zeebe & Wolf-Gladrow, 2001; we used their "recommended" coefficients for the calculation) and the TA from the empirical relationship with salinity (Taguchi et al., 2009). The TUMSAT and OCU data were measured at near water surface using the research vessels. The details are shown in Table 1 and the related studies (Endo et al., 2018;Kubo et al., 2017).
Because the NIES data for the bays were distributed uni-dimensionally along the course of the cargo ships, the spatial information for the data in this study was standardized by the distance parameter dist mentioned below. The distance parameter dist (km) was calculated as follows: where lat and lon are the latitude and longitude of the measurement point in radians, respectively. dist0_lat and dist0_lon are the latitude and longitude of the point representing the mouth of the river with the highest flow among the rivers in the inner part of each inner bay (Edo River in Tokyo Bay, Kiso River in Ise Bay, and Yodo River in Osaka Bay; Figure 1). The data for which dist > 100 km were excluded from analysis as being out-of-range of the inner and outer bays, as determined by changepoint analysis using the salinity and DIC (see Text S1). The change points were also defined as the boundary between inner and outer bay (71-77 km in Tokyo Bay, 62-64 km in Ise Bay, and 33-48 km in Osaka Bay, see Text S1). In addition, the data from Tokyo Bay with longitude >140°E were excluded because even though these data were within the range of dist ≤ 100 km, and they were from locations on the opposite side of a peninsula from the inner bay ( Figure 1). In total, we analyzed 17,977 data points from Tokyo Bay (  (2) where month and day are the month (1-12) and day (1-28 or 30 or 31) of the measurement, respectively. Additionally, given the temporal and spatial heterogeneity of the measurement data, we used natural neighbor interpolation (Sibson, 1981) to interpolate a grid with values at intervals of 1 km (dist) and 0.1 (monthlydata). In order for the interpolation to reflect the seasonal cycle from the minimum and maximum values for monthlydata (on 1 January and 31 December, respectively), the interpolation was performed with the data from the latter half of the year (July to December) added before the start of original data and that from the first half (January to June) appended to the end.
In order to minimize the effect of the trend in fCO 2 due to anthropogenic CO 2 input to the ocean and atmosphere, we corrected for the increase rate in each bay and in the additional data, except for the OCU data because those measurements were only taken for a single year. The increase rate for each bay and in the additional data was estimated from the linear regression using annual average of fCO 2 in water and air in each bay during the whole period. The increase in fCO 2 was corrected to that in 2010 if the increase was significant. The linear rates of fCO 2 increase in water were 5.83 and 3.24 μatm year −1 in Tokyo and Ise Bay, respectively, and were significant (p < 0.001), whereas those in Osaka Bay were not significant (p = 0.64). Also, the additional TUMSAT data for Tokyo Bay showed no significant increase (p = 0.85). Therefore, the rate of increase was corrected to the base year of 2010 only for NIES data for Tokyo Bay and Ise Bay. For fCO 2 in air, the linear increases were 2.1, 2.4 and 2.4 μatm year −1 in Tokyo, Ise and Osaka Bay respectively, and were significant (p < 0.001). As the same as fCO 2 in water, the trend of TUMSAT data was not significant (p = 0.22). Thus, the rate correction for fCO 2 in air was done for NIES data for Tokyo, Ise, and Osaka Bay.

Calculation of the Biological Effect
Because fCO 2 in water is affected by temperature and salinity, the biological effect cannot be quantified using only carbonate-system parameters. Although some previous studies normalized fCO 2 using an empirical relationship between temperature and fCO 2 (e.g., Takahashi et al., 2009), such a technique cannot be applied to coastal carbonate analysis because the salinity varies widely and affects TA and DIC non-linearly. Instead, we evaluated the biological effect using DIC, which is a conservative parameter, and whose variability is direct result of the biological effect from processes such as photosynthesis, respiration, and decomposition ( Figure S1). The DIC was calculated using the carbonate equilibrium calculation (Zeebe & Wolf-Gladrow, 2001) with fCO 2 and TA estimated from the empirical relationship with salinity in each bay from a previous study (Taguchi et al., 2009).
For quantifying the biological effect in bays, we defined DIC b as the increase or decrease in DIC through biological activities. The DIC b was calculated as the difference between DIC estimated from the equilibrium calculation and DIC ab , which was the interpolated value between the oceanic and riverine endmembers in the coordinate system of salinity and DIC (Tokoro et al., 2014). DIC ab represents the DIC concentration resulting from the mixing of the endmembers without any additional change in the inner bays such as through biological activity.
The endmembers and the biological effect were determined as follows. (a) the salinity and DIC of the oceanic endmember were estimated as the value outside of the respective outer bays. In this study, those values were estimated as the monthly average of the salinity and DIC between dist of 90 and 100 km. (b) We assumed that biological activity reached a minimum (DIC b ≈ 0) in the winter period (here, the 3 months with the lowest average water temperatures), according to the positive correlation between production by seagrass and algae and water temperature (Davison, 1991;Lee et al., 2007;Tait & Schiel, 2013). The assumption supported the measurement in Tokyo Bay in which the Chlorophyll a concentration (chl.a) tended to lower during winter period (Kubo et al., 2017). Then, the least-squares line of the salinity and DIC through the oceanic endmember average value during the above winter period can be used to approximate the DIC ab in the winter period and the riverine endmember (DICr, μmol kg −1 ). That is, the former parameter is the approximation using the regression line of DIC versus salinity, and the latter parameter is the regression line value when the salinity is zero. These parameters can be estimated as follows: Here, Sal w and DIC w are the salinity and DIC in each bay during the 3 months with the lowest average water temperature. Sal ow and DIC ow are the mean values of the salinity and DIC, respectively, of the oceanic endmember during these three months. Sal r is the salinity of the riverine endmember and assumed to be zero in this study.
(3) DIC ab in each month was calculated again as the linearly interpolated value between the endmembers. To evaluate the uncertainty of the estimated riverine DIC, we defined the range of DIC as ±200 μmol kg −1 (see Text S2) and calculated the precision of the range for the following procedures. (4) DIC b was calculated as the difference between the estimated DIC and DIC ab . (5) fCO 2 was calculated from the equilibrium calculation using DIC ab as the fCO 2 without any biological effects (fCO 2ab ) and fCO 2b (the difference between fCO 2 and fCO 2ab ) as the fCO 2 change due to the biological effects in each inner bay.
The air-water CO 2 flux (F, μmol m −2 s −1 or mol m −2 year −1 ) was calculated as follows: where k is the gas transfer velocity (m s −1 ) and calculated as described in the next paragraph. S is the solubility of CO 2 in water (mol m −3 atm −1 ) and estimated using an empirical equation using water temperature and salinity (Weiss, 1974). fCO 2water and fCO 2air are the fugacity of CO 2 in water and air (μatm), respectively. A positive value for the flux indicates a CO 2 efflux to the atmosphere, and vice versa. The air-water CO 2 flux due to abiotic factors (F ab ) and the contribution of biological effect to the flux (F b ) were also calculated using fCO 2water calculated from fCO 2ab and fCO 2b, instead of the fCO 2water , respectively. Because the temporal and spatial information for fCO 2air from TUMSAT did not completely correspond to that of fCO 2water measurements, the data were approximated by the data with the same dist-monthlydata grid information, which was calculated by natural neighbor interpolation using the original TUMSAT fCO 2air data. Because the OCU data did not include fCO 2air measurements, we used the average value of the NIES data for Osaka Bay from 2011 to 2016 (392.6 μatm) for the calculation.
The gas transfer velocity k (here, the units are cm hour −1 ) was calculated as follows (Wanninkhof, 2014): where U 10 (m s −1 ) is the wind speed at the height of 10 m from the water surface. Sc is the Schmidt number, defined as the ratio of the CO 2 molecular diffusion coefficient to the dynamic viscosity of seawater; Sc of seawater can be calculated from an empirical equation using the water temperature (Jӓhne et al., 1987). The wind data were taken from the database of the NEDO Offshore Wind System (NeoWins; http://app10. infoc.nedo.go.jp/Nedo_Webgis/top.html), which is the open database of the wind simulator provided by the New Energy and Industrial Technology Development Organization (NEDO) for evaluating the efficiency of offshore wind power generation in Japan. The wind data were supplied as the wind rose data, monthly averaged data, and averaged wind profile data, with 500-m resolution. We took the monthly averaged data at 10-km intervals beginning at dist = 0 along the water course (Table S1). Because the monthly averaged data were supplied as the data for 60-140 m above the water surface, while the wind profile was at 10-200-m height, we calculated the monthly averaged data at 10-m height according to the power-law of wind profile as follows: where U 60 (m s −1 ) is the monthly averaged data at 60-m height. U p10 and U p60 are the wind-profile data at 10-m and 60-m height, respectively (Pagon, 1935). The wind data were applied to the carbonate measurement data for each of the 12 months and 10 km of dist. Because the wind data was estimated from the monthly averaged data, the effect of the episodic high wind events like Typhoon was eliminated from the flux calculation.

Data Interpolation and Error Range
Because the measurement data were not homogeneous temporally or spatially, the averages and deviations described in this study were calculated from the interpolated data, which is a homogeneous grid of the distance (every 1 km dist) and date (every 0.1-monthlydata) data (n = 100 × 120 = 12,000). The error values for the carbonate parameters are shown in Table S2.
The procedure used for error estimation is described in the supporting information (Text S3). The error propagation was complicated because of the non-linearity of the carbonate equilibrium calculation, so the error was separated into the random error mainly due to the TA estimation and the riverine error due to uncertainty in estimating the riverine DIC.

Temperature, Salinity, and Carbonate Parameters
The basic measurement parameters (water temperature, salinity, and fCO 2 ) in each bay were shown in Figure 2 and Table 2. The estimated TA and DIC were 2,152 ± 116 and 1,862 ± 129 μmol kg −1 (Tokyo Bay), 2,166 ± 115 and 1,883 ± 120 μmol kg −1 (Ise Bay), and 2,179 ± 64 and 1,893 ± 85 μmol kg −1 (Osaka Bay) (Ave. ± 1SD) ( Figure S2). The variability of water temperature and salinity reflected the typical seasonal trend on the North Pacific coastal in Japan such as the wet summer and dry winter. The variability of fCO 2 also showed the seasonal trend such as the spring-summer decrease at landward area but the temporal and spatial distribution was different among the bays.
The salinity and DIC of the oceanic endmembers for the calculation of abiotic and biotic parameters were estimated for each month as the average values between 90 and 100 km from the river mouth reference point in each bay (Table 2). These values were higher during winter and lower during summer and were consistent with the results from an empirical equation derived for the Kuroshio stream area (Ishii et al., 2011) ( Figure S3). The DIC of the riverine endmember was also shown in Table 2.
The DIC ab and DIC b were 1,892 ± 83 and −30 ± 69 μmol kg −1 in Tokyo Bay, 1,896 ± 97 and −13 ± 28 μmol kg −1 in Ise Bay, and 1906 ± 55 and −13 ± 39 μmol kg −1 in Osaka Bay, respectively ( Figure 3). The fCO 2ab and fCO 2b were 333.43 ± 21.88 and −16.69 ± 107.66 μatm in Tokyo Bay, 314.39 ± 13.78 and −21.64 ± 31.08 μatm in Ise Bay, and 311.23 ± 15.50 and −12.54 ± 44.47 μatm in Ise Bay, respectively. The negative DIC b indicates that the DIC decrease by biological activities (e.g., photosynthesis) exceeded the DIC increase by other activities (e.g., respiration, degradation). The DIC ab decreased within 60-70 km from land in summer in all bays. The average and deviation of DIC ab in Tokyo and Ise Bay were almost the same while those in Osaka Bay showed a little high average and small deviation. While the estimated DIC b showed a significant decrease within 70 km from land in summer as the same as DIC ab . Meanwhile, there was an increase in Tokyo Bay within about 10 km of land from October to January. The seasonal trend of variability in DIC b was almost the same among the bays. The magnitude and spatial distribution of the decrease in DIC b was the greatest in Tokyo Bay, where the increase in DIC b were evident mainly in autumn to winter. In the other two bays, the magnitude and spatial distribution of DIC b were almost the same and smaller than in Tokyo Bay.

Air-Water CO 2 Fluxes
The air-water CO 2 fluxes in the three bays indicated that these areas were annual atmospheric CO 2 sinks (−2.86 ± 2.70, −3.13 ± 1.50 and −2.37 ± 1.70 mol m −2 year −1 in Tokyo, Ise and Osaka Bays, respectively) ( Figure 4). The temporal and spatial distribution of the air-water CO 2 flux in Ise and Osaka Bays were similar and showed a seasonal pattern whereas the distribution in Tokyo Bay seemed to be more consistent TOKORO ET AL.
10.1029/2020JC016974 with that of DIC b than a seasonal pattern (Figure 3). The peak of CO 2 absorption was during February to April, and CO 2 absorption was lowest in August and September in Ise and Osaka Bays, and in October and November in Tokyo Bay.
The fluxes without (F ab ) and from biological activity (F b ) were −2.09 ± 1.92 and −0.77 ± 3.12 mol m −2 year −1 (Tokyo Bay), −2.75 ± 1.70 and −0.38 ± 1.09 mol m −2 year −1 (Ise Bay), and −2.24 ± 1.49 and −0.13 ± 1.09 mol m −2 year −1 (Osaka Bay), respectively ( Figure S4). The standard deviations of the biotic fluxes (F b ) indicated that the biological effect in Ise and Osaka Bays was weaker than in Tokyo Bay. The temporal and spatial distribution of F ab showed a seasonal pattern whereas the distribution of F b seemed to be related to the distribution of DIC b .

Discussion
The validity of our estimates of biological effects on DIC (DIC b ) and CO 2 flux (F b ) depended on the precision of the riverine DIC endmember because the abiotic DIC (DIC ab ) and CO 2 flux (F ab ) were determined from the riverine endmember, along with the oceanic endmember whose error was relatively small. Although we could not find reasonable reference data for riverine DIC in the bays, the reported riverine TA data by Taguchi et al. (2009) which are 1,006, 518 and 759 μmol kg −1 in Tokyo, Ise, and Osaka Bays, respectively, support the validity of the DIC estimation because water with a higher TA can contain more carbonate and bicarbonate ions as DIC under the same fCO 2 conditions. The range of estimated riverine fCO 2 in each bay also supports the validity of our estimates of riverine DIC and its precision (±200 μmol kg −1 ) (see Text S2).
The negative average annual biotic DIC in all three of the bays indicates that the ecosystems were net autotrophic. The magnitude of the biotic DIC seemed to be consistent with the nutrient concentrations and chl.a in the inner bay which was reported in the reference material of the Ministry of the Environment, Japan (https://www.env.go.jp/council/09water/y0917-07/ref02.pdf). In this report, the highest concentrations of total nitrogen, phosphorous, and chl.a were in Tokyo Bay and were almost the same in Osaka and Ise Bays (average total nitrate, 0.8, 0.3, 0.4 mg L −1 ; average total phosphate, 0.06, 0.03, 0.04 mg L −1 ; average chl.a, 10-20 mg m −3 , 5-10 mg m −3 , 5-10 mg m −3 , of 10-50 sampling points in the inner bays from 1981 to 2012 in Tokyo, Ise, and Osaka Bays, respectively) ( Table 3).
The relationship between the autotrophic condition and the nutrient concentration is consistent with the hypothesized mechanisms related to wastewater treatment (Kubo & Kanda, 2020;Kubo et al., 2017;Kuwae et al., 2016). Typical wastewater treatment method removes carbon in the form of sludge and CO 2 gas more TOKORO ET AL.    efficiently than nutrients such as nitrogen and phosphorus (Sedlak, 1991). Through these effluents, the balance of primary production and respiration in the inner bay is offset toward an excess of primary production and the resultant suppression of fCO 2 increase. Second, the effluent contains relatively refractory carbon, because labile organic matter has already been removed by treatment (Kubo et al., 2015). Therefore, respiration and mineralization rates of effluent are low, and subsequent fCO 2 increases are suppressed. The reported ratios of C:N:P in the treated water from the largest treatment plant in each bay were about 6:21:1, 26:81:1, and 10:28:1 in Tokyo Bay, Ise Bay, and Osaka Bay, respectively (Japan Sewage Works Association, 2010) ( Table 3). Note that the value of C was estimated from the COD as the approximated concentration of labile organic matters. This indicates that the decomposable organic carbon in treated water was far less than the organized DIC, which was estimated from the nutrient concentration and the Redfield ratio. Although COD do not equal to labile organic matters concentration strictly, the ratio supports the decrease in DIC and fCO 2 in the bays due to the wastewater treatment.
The nutrient concentrations in the three bays were probably determined by the volume of treated wastewater discharged into the bays, which is related to the size of the population of the surrounding urbanized area (Table 3). For example, the inflow of the treated water is reported as the main component of the total freshwater inflow to Tokyo Bay and thus the effect of treated water should be noticeable on the nutrient concentration in the bays (Kubo et al., 2015). Although the net primary production in the bays depends on several parameters such as the seawater residence time and vertical stratification, its general magnitude would be similar to that of the biotic DIC because it is similarly influenced by hydrographic conditions.
It is possible that the net primary production in Osaka Bay was underestimated because Osaka Bay connects with the Seto Inland Sea at its landward end and this topography results in strong tidal currents in the bay (Odamaki, 2002). Stronger tides would enhance the water exchange between the bay and the surrounding areas and weaken the biological effects on DIC and the air-water CO 2 flux (Table 3). In addition, the area at dist of 50-100 km coincides with narrow straits (Kitan Strait and Kii channel) where several large rivers flow into the bay. Thus, the calculation of the endmember effect in Osaka Bay might be biased compared with those for the other two bays.
The CO 2 fluxes found in this study indicate that the inner bays and surrounding outer bays in Japan are one of the largest atmospheric CO 2 sinks among the global coastal areas reviewed in previous studies (Borges et  Main river discharge (m 3 s −1 ) b ∼100 (Edo River) ∼200 (Kiso River) ∼200 (

Table 3
The List of Information Related to the Nutrients and Biotic Parameters in Each Bay more atmospheric CO 2 absorption than the average in these previous studies for estuaries (7.74-10.26 mol m −2 year −1 ) and marginal seas (continental shelves) (−1.64 to −1.06 mol m −2 year −1 ). The absorption in the bays was mainly based on the abiotic flux (73%, 88%, and 94% of the net CO 2 flux (F) in Tokyo, Ise, and Osaka Bays, respectively). The oceanic endmember from the Kuroshio stream area is the most plausible explanation for the CO 2 absorption (overall average CO 2 absorption of −2.23 mol m −2 year −1 ). That area has been reported as the one of the largest CO 2 sinks in the world because of the cooling effect of cold Oyashio water on warm Kuroshio water (Takahashi et al., 2002(Takahashi et al., , 2009. We suggest that the CO 2 absorption in the bays was enhanced by additional cooling due to the terrestrial effect during winter (Figure 2).
The largest biotic CO 2 flux in Tokyo Bay was consistent with the largest biotic DIC and nutrient concentrations in the bay. The higher CO 2 absorption in Ise Bay compared to Osaka Bay despite almost identical biotic DIC can be explained by the lower salinity and TA in the near-shore area in Ise Bay (Figure 3 and S2). Because lower-TA water has less buffering effect on fCO 2 , the decrease in fCO 2 in Ise Bay was greater than that in Osaka Bay even when the decrease in DIC was the same. Although the annual average of the biological effect was limited (27%, 12%, and 6% of the net CO 2 flux in Tokyo, Ise, and Osaka Bays, respectively), it affected the temporal and spatial distribution of the air-water CO 2 flux. Temporally, the pattern of F b was the opposite of F ab , with an influx (F b ) or efflux (F ab ) in summer and vice versa in winter ( Figure 5). The biological effect was the strongest in Tokyo Bay and less notable in the other bays, as with the biotic DIC. As a result, the seasonal variation of the CO 2 flux in Tokyo Bay was different from that in the other two bays despite having almost the same variation in the abiotic CO 2 flux ( Figure 5). For example, the peak CO 2 influx in winter in the bays continued into summer in Tokyo Bay, but the winter influx in Tokyo Bay was about two-thirds that in the other bays.
The difference in the spatial distribution of F b among the bays was less noticeable than that in the temporal distribution because of the offset of the influx in summer and efflux in winter in Tokyo Bay ( Figure 5). The magnitude of the influx peak at dist of about 20 km was almost the same among the bays whereas the efflux at dist > 50 km was the largest in Osaka Bay. This might be caused by the decomposition of the organic matter produced by the photosynthesis at about 20 km dist.
Meanwhile, the effect of CO 2 release to atmosphere in near-shore area was shown within ≤5 km area (Figure 4) although fCO 2 of the riverine endmember was estimated at more than 2,000-5,000 μatm in the case of the intermediate riverine DIC. As for Tokyo Bay, there was a tendency toward an efflux in the near-shore area both in F ab and F b ( Figure 5 and S4), corresponding to the decrease in the DIC ab in summer and the increase in DIC b in winter (Figures 3 and 4). Because the distribution of the decrease in DIC ab was consistent with the salinity distribution (Figures 2 and 3), an increase in riverine flow from precipitation might cause a CO 2 release as an abiotic factor. On the other hand, the increase in DIC b was consistent with the increase in pCO 2 observed in a previous study (Kubo et al., 2017) in which the increase was due to the weakening of stratification in Tokyo Bay due to the cooling of surface water, a decrease in precipitation, and the weakening of the seasonal southward wind. These factors probably contributed to bringing high-fCO 2 water and resuspended organic sediments from the bottom to the surface, resulting in the CO 2 efflux to the atmosphere. Although we could not perform a similar analysis for the other two bays because of a lack of measurement data, we would expect a similar, considerable efflux tendency because the other bays share the same hydrographic and climate conditions.
Although the measurement data in this study was limited in the bays open to North Pacific Ocean, the relationship between the coastal urbanization and biotic carbonate parameters should be applicable to wide variety of regions because the regional differences in temperature and salinity are basically smaller than the variabilities due to the seasonality and the mixture of riverine water. Therefore, the comparison between Tokyo Bay and other two bays can suggest that the enhanced urbanization (increase in population and area of surrounding urban region) with wastewater treatment results in increased biotic CO 2 absorption. Because the development of coastal areas will likely continue for several decades, the biotic absorption of CO 2 in the bays is expected to be a mitigating factor for future climate change. We assumed that the area within 100 km of the global coastline (6.2 × 10 7 km 2 ) could absorb atmospheric CO 2 additionally at the same rate as the biotic flux in Tokyo Bay (0.77 mol m −2 year −1 ), and roughly estimated the potential for additional biotic absorption to be 0.057 Pg year −1 . This is on the same order of magnitude as the estimated estuarine CO 2 efflux (e.g., Chen et al., 2013). The estimation of worldwide abiotic CO 2 flux is difficult because the oceanic and riverine carbonate parameters differ at each location. In addition, the carbon export from urbanized areas is observed to increase along with the development (Barnes & Raymond, 2009;Lopes et al., 2020;Wang et al., 2017). This increase in carbon export could mitigate the increase in the CO 2 absorption.
The estimated consequences related to the CO 2 absorption should be evaluated for total benefit of the coastal land use. For example, the annual increase rate in fCO 2 in Tokyo Bay (5.83 μatm year −1 ) and Ise Bay (3.24 μatm year −1 ) might indicate the enhancement of the acidification in the bays although the long-term trend of the increase was obscure due to non-strict linearity (r 2 = 0.88 in both bays). On the other hand, the net biotic CO 2 absorption does not enhance the acidification because the absorbed CO 2 was converted to organic components and not decomposed near water surface. However, if the organic components were decomposed at bottom layer, the biotic flux would result in the enhancement of the hypoxia or HABs in the bay. For balancing the CO 2 absorption and the consequences like hypoxia and HABs, the regulation of phosphate in treated wastewater should be efficient considering the ratio of N:P (21-81:1) than the Redfield ratio (16:1) (Lui & Chen, 2012).
For accurate evaluation of the global CO 2 absorption and the consequences resulting from coastal urbanization, it is necessary to obtain more carbonate measurements from a variety of areas. For example, data from measurements at the lagoons in Ivory Coast and Guanabara Bay in Brazil suggested that strong CO 2 absorption and release are mixing in the tropical urbanized inner bay where the rate of wastewater treatment is low (Cotovicz et al., 2015;Koné et al., 2009). Both of these studies indicated that eutrophication and stratification were the factors regulating the undersaturation of CO 2 . Meanwhile, studies in the Chesapeake Bay reported that both the absorption and release of CO 2 were observed by numerical simulation model and abundant pH measurement, respectively (Herrmann et al., 2020;St-Laurent et al., 2020). Because the main bottleneck for further study is the difficulty of obtaining comprehensive measurements in coastal areas, the development of novel methods like the above Chesapeake Bay studies will facilitate filling gaps in the temporal and spatial distributions of available data. The review study of such measurements is expected to understand the quantitative relationship between the coastal urbanization and the carbon cycle.

Conclusions
In this study, we quantified the temporal and spatial variations of carbonate parameters in three urbanized inner and outer bays in Japan including riverine water and the surrounding oceanic water within a range of 100 km. Our results are the first in the world to fill the gap between studies of nearshore areas and the marginal seas. We found a notable atmospheric CO 2 influx due to biological activity in the areas about 20 km from land in summer, accounting for 27%, 12%, and 6% of the measured net CO 2 flux in Tokyo, Ise, and Osaka Bays, respectively. In addition, the biological effect in the highly developed Tokyo Bay significantly affected the seasonal variation of air-water CO 2 flux. The potential for atmospheric CO 2 absorption was mainly regulated by the water in the surrounding marginal seas. The biological effect seemed to be associated with the nutrient concentrations, which are related to the volumes of treated wastewater entering the inner bays. In addition, the degree to which the water area is enclosed likely influences the biological effect. Conditions that are more closed or open would enhance or mitigate the biological effect on the air-water CO 2 flux, respectively.
This study should contribute to future investigations into the carbon cycle in urbanized coastal areas, which will likely continue to expand for the next several decades. However, more detailed investigations in inner and outer bays are required for more precise evaluation of their contribution to the global carbon cycle. To expand the results of this study to the global scale will require further measurements in bays in a variety of regions. We think that the measurement of water temperature, salinity, fCO 2 in water and the quantification of the linear relationship between salinity and TA are the crucial and efficient for the analysis of the carbonate parameters and air-water CO 2 flux. The more VOS programs in bays should provide abundant dataset for the analysis. Furthermore, the field investigation of the above parameters in riverine water is required for quantifying the biotic and abiotic components like this study.