The ocean's overall role as a sink for anthropogenic carbon (Canthro) has been determined unequivocally via global air-sea carbon dioxide (CO2) flux observations (1), global oceanic anthropogenic carbon inventory observations (2)(3)(4)(5)(6), and ocean physicalbiogeochemical model simulations (7,8). As such the ocean has played an important role in shouldering CO2 storage and mitigating global warming but also has suffered from the acidification stress resulting from the absorbed CO2 (9,10). However, it has not yet been revealed how much Canthro accumulated in the coastal ocean and where it comes from using observational data within coastal oceans.
Though small in surface area, the coastal regions, including estuaries, tidal wetlands, and the continental shelf and upper slope, are important locations for global carbon budget with high variability and uncertainty due to the heterogeneity within and among coastal systems (11)(12)(13). After two decades of intensive field and modeling research, it is clear now that the continental shelf is presently a sink for atmospheric CO2 while estuaries are a CO2 source driven by both the river and wetland carbon export (14,15). Previous models and conceptual frameworks have suggested that the continental shelf is an increasing CO2 sink, either switching from a source (16) or being enhanced from a weak sink in the preindustrial time to a stronger sink at present (17), implying the continental shelf as a potentially critical location to take up anthropogenic CO2 and to export it to the open ocean (18,19). Recent syntheses of observational data and model results have determined that the global atmospheric CO2 uptake rate is about 0.25 PgC yr -1 (11,14,20) for the continental shelf (to 200 m). However, so far, we have no evidence to conclude on how much of this contemporaneous CO2 uptake flux is stored as Canthro, though conceptual model, mass-balance analysis, and numerical models often ignore coastal gradients and heterogeneities and suggest a substantial amount of the total flux is anthropogenic CO2 uptake by the coastal ocean (18).
An even less known but immensely important question considers whether there is a lateral transport of anthropogenic carbon from the terrestrial system to the open ocean along the land-ocean aquatic continuum (LOAC). A recent analysis by Regnier et al. (20) has suggested that, globally, net Canthro flux from the estuaries-wetlands to the continental shelf over the industrial period is essentially zero (0.05±0.10 PgC yr -1 ). This is a downward revision of the LOAC Canthro flux of about 0.1 PgC yr -1 previously proposed by Regnier et al. (2013) as it is believed to be compensated for by the more recent organic carbon loss in coastal wetlands (also see Dai et al. 2022(15)). However, the above conclusions are largely derived from indirect mass balance assessments and no direct observational evidence has been provided to support or negate them. In contrast, another model-based global river carbon budget analysis showed the carbon export from river to ocean (~0.9 PgC yr -1 , in which DIC is ~0.6 PgC yr -1 ) has been stable, meaning no net Canthro, during the 20 th century (21). Furthermore, a recent global synthesis of riverine bicarbonate (HCO3 -) export is also inconclusive regarding whether anthropogenic activities have increased river HCO3 -flux on a global scale (22), though other syntheses have suggested that part of the terrestrial organic carbon export could be converted to dissolved inorganic carbon (DIC) and together with the CO2 uptake from the atmosphere, they have increased the shelf to open ocean export of DIC (11,14). In summary, the input and accumulation rates, variabilities, and controls of the Canthro in coastal oceans are still unclear. An assessment based on direct observations of anthropogenic carbon accumulation in downstream coastal oceans may provide a new angle to address this unresolved question which is of broad interest to ocean and global carbon cycle and climate change communities.
To distinguish Canthro from the natural background in the coastal ocean is challenging as only the total carbon changes can be directly measured and as the ratio of signal (anthropogenic trend) to noise (natural variability) may be small in coastal regions due to high natural variability and limited period of observational data (8,23). Hence coastal Canthro estimation has so far mainly relied on either biogeochemical model simulations (17)(18)(19) or a combination of conceptual model understanding and mass-balance analysis (19,20). While attempts of estimating Canthro have been made on large marginal sea basins such as the South China Sea (24,25) and the Mediterranean Sea (26), these isolated deep basins are treated like the open ocean basins, involving no data from low salinity shelf and nearshore waters. In addition, Canthro content was also assessed in the upwelling dominated eastern boundary ocean margins (27,28), yet it was extrapolated from the decadal repeated carbonate observations of the open ocean to the coastal waters (29,30).
While Canthro was also estimated and its impact on acidification was assessed in nearshore systems, such studies were generally limited to surface (31,32) or seasonally hypoxic shallow subsurface water (33). Thus, high-quality synoptic coastal cruise data have barely been used to estimate regional Canthro directly so far. Although using the high-quality and long-term carbonate parameter observational data within shallow coastal oceans is the preferred path to move forward, time-series investigations of accurate and precise carbonate parameters of sufficiently long periods (> two decades) are rare in coastal ocean. To our knowledge, the North American East Coast may be the only ocean margin, among those seriously impacted by terrestrial carbon exports, which meets the criteria with a time series of high-quality carbonate chemistry investigations that can be traced back more than 20 years to the mid-1990s (34)(35)(36)(37).
Here, we take advantage of such a dataset from 1996 to 2018 to estimate the Canthro change (ΔCanthro) and explore its controlling processes in the U.S. Mid-Atlantic Bight (MAB) (fig. S1). Our results show that ΔCanthro is nearly zero in the low salinity (S<32) nearshore waters, meaning no increase in the Canthro export from estuaries-wetlands. We further conclude that only a small fraction of anthropogenic CO2 uptake by the continental shelf has been stored within the shelf while the majority is exported to the open ocean, aligning with recent model results. These conclusions may have major implications for Canthro re-distribution among the Earth systems. Our conclusions may stimulate a community-wide level effort to quantify coastal Canthro storage rates as more and more high-quality and long-term carbonate parameter data become available in the highly heterogeneous global coastal oceans in the next two decades, improve the constraint of boundary conditions in regional and global ocean carbon models, and potentially provide the necessary baseline for assessing the marine carbon dioxide removal (mCDR) intervenes.
We first identified a consistent spatial distribution pattern and an evolving temporal change of DIC from 1996 to 2018 across five research cruises in the MAB (Fig. 1, fig. S1, table S1). The spatial variations of DIC content are largely determined by their source water properties (fig. S2 -S3). On the shelf, DIC concentrations are low in the low salinity nearshore water and the low salinity and cold shelf water from the alongshore Labrador Current (LC) (38,39) (Fig. 1A -1C and fig. S2). DIC concentrations increase moving offshore, especially in waters proximity to the Gulf Stream (Fig. 1A and 1B and fig. S3). On the continental slope, DIC is low in the top 200m of the high-temperature and salinity Gulf Stream water (GS) (Fig. 1 S2). However, due to the large natural variations, additional analysis is needed to isolate how much DIC changes is contributed from ∆Canthro.
To identify the anthropogenic effect on the secular trends of DIC changes, we estimated ΔCanthro and its increasing rate in the subsurface using a time series analysis of regional extended multiple linear regression (regional eMLR) (42,43) built from salinity, temperature, oxygen, and silicate concentration (Methods, table S3 and table S4). The results are further validated by the ensemble eMLR approach(30) (table S5 and table S6).
While similar in principle to the eMLR approaches developed and used for the open ocean basins, regional eMLR includes a cross term to deal with the more dynamic coastal oceans (42,43). On the surface (depth < 20 m), we estimated ΔCanthro using the constant air-sea disequilibrium (ΔpCO2 = seawater pCO2 (pCO2sw) -atmospheric pCO2 (pCO2air)) method to avoid the pitfalls of the eMLR approach(3).
We found a strong ΔCanthro vertical variation (Fig. 2A and 2B) in different slope water masses and distinct increasing trends with time (Fig. 2C). The highest ΔCanthro (26 ± 2 µmol kg -1 ) over the period of 22 years (or annual increasing rate of ~1.1± 0.1 µmol kg -1 yr -1 ; table S6) is located in the top 200 m, where the water shows GS features with high salinity, temperature, and TA. Below the GS water, the lowest ΔCanthro of 9 ± 3 µmol kg -1 (or annual increasing rate of 0.3 ± 0.1 µmol kg -1 yr -1 ) is found in the AAIW+ water between 200 m to 500 m. Below 500 m, ΔCanthro in the uNADW water is 16 ± 3 µmol kg - 1 or an annual increasing rate of 0.5 ± 0.1 µmol kg -1 yr -1 (table S5 and S6). This sandwich structure of ΔCanthro in the slope water reflects the history and efficiency of the anthropogenic carbon uptake of the source waters, depending on the water-mass age (i.e., the time of last contact with the atmosphere) and water chemical properties. We illustrate and explain the mechanism on a ΔCanthro -pCO2sw space (Fig. 3) with the assumption that the air-sea CO2 gradient, ΔpCO2, in each source water has remained constant over the recent few decades(3). Here, each group of contours represents the ΔCanthro corresponding to a historical atmospheric CO2 increase in different water masses.
On the slope, the GS water-mass age is within one year because it is in the surface and upper water column (44,45). Therefore, under full air-sea equilibrium, we estimated that the GS water with high salinity and TA should have taken up about 26 µmol kg -1 of ΔCanthro between 1996 and 2018 when the pCO2air increased from 363 to 408 ppm (Fig. 3, red circles and solid arrow), which is compatible with the eMLR-derived ΔCanthro of ~26 ± 2 µmol kg -1 (table S5).
The deep uNADW, originated from the Labrador Sea by winter-time convection (46), takes about 18-22 years for the surface Labrador Sea water to subduct to 1000 m depth and move to the MAB region along isopycnals (46). Therefore, the uNADW ΔCanthro accumulation signal from 1996 to 2018 actually reflected the air-sea exchange state of the Labrador Sea surface seawater around two decades before the respective survey times.
Theoretically, ΔCanthro accumulation in the uNADW should be about 16 µmol kg -1 for the 20-year period when pCO2air increased from 330 ppm in 1976 to 364 ppm in 1998 (Fig. 3, the blue circles and solid arrow). However, not all uNADW in this region is newly convected (46), and thus the actual ΔCanthro could be lower. Our eMLR-derived ΔCanthro (16 ± 3 µmol kg -1 , table S5) from observational data closely matches this theoretical prediction within uncertainties. Finally, the AAIW+ water originated from the Antarctic has a water mass age of at least 50 years (47). As pCO2air increase was small (~ 7 ppm or less) during a 22-year time span in the mid or early part of the 20th century (48), the ΔCanthro accumulation rate is very low and indistinguishable from the uncertainty of estimation.
For the shelf water (water depth < 200 m), we used the same approaches (regional eMLR approach in the subsurface and constant ΔpCO2 method in the surface) to evaluate ΔCanthro. We also tested other methods of estimating ΔCanthro in the surface layer, though variance is shown in ΔCAnthro content in low salinity water, they all generated similar distribution pattern and conclusions (Supplementary Material section 4). Vertically, we found that the ΔCanthro is higher in the surface (σ < 24) than in the subsurface (Fig. 4A -4D). Horizontally, the most distinguished feature is a decreasing ΔCanthro from offshore to nearshore. Thus, both cross-shelf mixing and air-sea CO2 exchange are important factors influencing the distribution of Canthro in continental shelf waters.
The relationship between salinity and ΔCanthro (Fig. 4E) can be used to separate water mass mixing and ΔCanthro storage and to explore the mechanisms behind the complex distribution pattern of shelf water ΔCanthro. Once a water parcel leaves the surface and is isolated from the atmosphere, its ΔCanthro is mainly affected by the water's physical mixing and can be considered quasi-conservative. Thus, in the subsurface, a triangle shape distribution of ΔCanthro against salinity (red triangle in Fig. 4E) indicates a three end-members mixing of the low-salinity nearshore water from estuaries, the mid-salinity southward alongshore LC water, and the high-salinity northward GS water. Interestingly, the two coastal end-members have distinctly different ΔCanthro (fig. S7). The alongshore LC water has ΔCanthro of 17 ± 4 µmol kg -1 (Fig. 4E, table S5), matching our theoretical ΔCanthro estimation of ~19 µmol kg -1 from 1996 to 2018 at S = 33, T =10°C (Fig. 3, blue square). However, in stark contrast, the nearshore low salinity end-member originated from local estuaries and wetlands (S < 30, T = 25°C) has a near-zero ΔCanthro (6 ± 7 µmol kg -1 , Fig. 4E and table S5) between 1996 and 2018, which greatly deviates from the ΔCanthro estimation of 20 µmol kg -1 under air-sea gas equilibrium assumption.
A near-zero ΔCanthro in the nearshore subsurface water from estuaries can first be attributed to two characteristics of its source water when it is on the surface: 1) high sea surface pCO2 and low buffer capacity, and 2) short water residence time. We again illustrate and explain the control mechanism on the Canthro -pCO2sw space (Fig. 3), focusing on surface water with high pCO2. The increase of DIC caused by the same pCO2air increase becomes smaller when the water pCO2sw is higher and water's buffer capacity is lower (Fig. 3, represented by the changing slopes of the contours). For example, high seawater pCO2 (500 ppm) in the low salinity nearshore water can only allow a small amount of ΔCanthro accumulation (~12 µmol kg -1 ) corresponding to the 45 ppm atmospheric CO2 increase between 1996 and 2018 (Fig. 3, red squares and dark red solid arrow). In addition, the short water residence time also limits ΔCanthro storage in the nearshore waters. Near the estuarine mouth, the water residence time is usually short (10 -100 days) (49)(50)(51), which causes a faster DIC exchange with the end-member waters (here we used a mixture of estuarine water and the offshore water). We simulated the effect of residence time using a simple 1-D box model, which assumes constant physical, chemical, and biological parameters but allows the pCO2air to increase under a suite of different mixing scenarios from the marine-dominant case (10% estuarine and 90% marine) to the estuarine-dominant case (90% estuarine and 10% marine) (see Method, fig. S8). The simulation results illustrated that ΔCanthro in 10-day and 100-day residence time scenarios can accumulate only about 20% and 40%, respectively, of the ΔCanthro values with the infinitely long residence time scenario derived from the estuarinedominant case (Fig. 3, yellow and orange arrows, fig. S8).
In the surface water, the ΔCanthro, though also decreasing with the salinity decrease, is higher than the subsurface ΔCanthro in the low salinity end, which indicates that other than water-mass mixing, an extra source contributes to the surface water Canthro (Fig. 4E). If
Our results provide the first solid evidence in constraining the anthropogenic carbon flux in the continental shelf realm based on the long-term oceanic observational dataset. This observational conclusion further clarifies and fortifies the existing continental shelf carbon budget analyses, previously only reliant on models, mass balance analysis, and conceptual frameworks.
First, the derived near-zero ΔCanthro in the low salinity nearshore water (Fig. 4) indicates that the contribution of the ΔCanthro from estuaries and wetlands to the shelf is rather limited. The ΔCanthro could be negative when it is extrapolated to an even lower salinity (Fig. 4 and Fig. 5). This conclusion supports the mass balance analysis results of Regnier et al. (20) that net Canthro flux from estuaries and wetlands to the continental shelf is essentially zero and can be negative (0.05±0.10 PgC yr -1 ) as the riverine Canthro flux (0.10±0.05 PgC yr -1 ) is counter balanced by the reduced wetland area and its associated carbon loss and reduced CO2 uptake from the atmosphere within the past few decades. This notion is supported by the report that the North America has a wetland loss rate 3 times of its growth rate(52). This conclusion is also consistent with another global river analysis which suggests no obvious change of riverine carbon export during the 20 th century (21). The possibility of the estuaries being a site of reduced/removing Canthro uptake requires additional independent assessments but, if proven to be true, it will have important ramifications for global carbon cycle studies under a changing climate, in particular in the context of current interest in considering coastal wetlands and estuaries zone as viable storage reservoirs for carbon dioxide removal (CDR) (53).
Another insight from our study worthy of further exploration is the fate of Canthro from the local atmospheric CO2 uptake. Here we first quantify that the volume-integrated total anthropogenic carbon inventory change rate is about 0.055 Tg C yr -1 in the MAB shelf, which is only about 3.4% of total air-sea CO2 uptake (1.6 Tg C yr -1 ) (see Methods) (Fig. 5). Assuming that the air-sea Canthro flux is 30% -70% of the total air-sea flux (15,18,20,54,55), then we can conclude that only 5% -11% the air-sea Canthro flux is accumulated in the water column and that the majority of Canthro (89% -95%) must be exported to the open ocean. Even if we take the perturbed particulate organic carbon (POC) burial estimations into consideration (10% -22%, twice of the Canthro in the water column) (17,20,56), the majority of Canthro uptake (68% -85%) would still be exported. This observation-based conclusion agrees with the model result of Bourgeois et al (18), though our export percentage is even higher (Fig. 5). We further simulated the Canthro accumulated in the water and the Canthro uptake from the atmosphere using the same 1-D box model assuming pCO2air increase is the only driver under a suite of different mixing scenarios (See Method). The box model simulation also suggested that only 4% -28% of Canthro uptake through air-sea flux is accumulated within the shelf due to its limited volume and relatively short water residence time. The majority of this local gasexchange-induced Canthro (72% -96%) is exported to the open ocean in the MAB (Fig. 5).
Therefore, we further concluded that while the continental shelf is not efficient in accumulating Canthro, it is highly efficient in exporting its locally absorbed anthropogenic CO2 from the atmosphere to the open ocean and thus may play a critical role in the ocean acidification pathway.
Our observation-based findings address the questions of how and how much ΔCanthro is accumulated in the coastal region and where it comes and goes, and substantiated recent assessments from numerical models, mass balance exercises, and conceptual frameworks.
Our conclusions will also have implications for setting up or verifying boundary conditions in regional and global carbon cycle models and budgets (57). Though some previous studies have shown that Canthro accumulations in coastal regions can be influenced by overturning circulations (58), water ventilation (59,60), upwelling(61) and water buffer capacity (62), our study specifically addresses the Canthro accumulations mechanism in an open continental shelf with substantial influences from the terrestrial and estuarine sources and the open ocean. However, the coastal region is diverse and dynamic with local features. Consequently, accurately quantifying the contribution and storage percentage of ΔCanthro in different types of continental shelves globally (FCS, Fig. 5) remains challenging, given these complexities, but are highly needed. Our study only examined one example of using high-quality summer coastal cruises to conduct independent and high-resolution Canthro estimation and inventory changes. The coastal cruises usually have higher temporal and spatial sampling frequency, which can reveal the region-specific Canthro features and detect Canthro evolution. The method established and validated in this study can potentially be used in different coastal regions and enhance the understanding of past, present, and future coastal Canthro changes. Especially relevant and timely, the application of this method can benefit the understanding of coastal Canthro dynamics, which is critical to evaluating and verifying the efficacy of the various marine CDR strategies in the following decades.
We also advocate that future research should quantify and understand the contributions of anomalous CO2 fluxes (that is non-steady state changes) from estuaries to the continental shelf (FEC, Fig. 5) (see Supplemental Material section 1 and 2) due to the changes in the wetland area, nutrient concentration, increased chemical weathering under anthropogenic pressure and climate change stresses (12,22,63), and terrestrial organic and inorganic carbon inputs forced by hydrological cycle (e.g., wet vs. dry years, extreme events, damming and etc.( 64)) and coastal currents, which may substantially alter the Canthro accumulation in rates and status on a regional scale. For example, anthropogenic activities have increased the net fluxes/burial of mud-associated organic carbon in some coastal environment but probably lead to a net decrease in coastal wetlands (52,64). Even though most Canthro estimations have assumed that Canthro does not interact with organic carbon cycle (65), there are exceptions in this assumption, such as in the consideration of ocean acidification (OA) caused by the Canthro, one also considers how OA modifies the natural carbon cycle and what are their synergistic impacts (31,(66)(67)(68). This will need further clarification on whether the low ΔCanthro from the estuarine has been transformed into other forms of carbon and stored within the estuaries before it is exported to the continental shelf. A further collection of long-term and high-quality observational data and improvement of realistic and high-resolution global and regional coastal ocean models, and development of methodology in quantifying the component of non-steady state carbon and its interaction with climate change trend are clearly needed to validate and improve our conclusions and implications. With more high-quality and comprehensive data available in coastal oceans in the coming decade, we will be able to piece together a more complete picture of the contributions of the coastal ocean in anthropogenic carbon storage, ocean acidification, the global carbon cycle.
We used data from five research cruises in the Mid-Atlantic Bight (MAB) conducted in the summer of 1996 (Ocean Margin Program, OMP Cruise, http://po.msrc.sunysb.edu/omp/) (36), 2007 (Gulf of Mexico Ecosystem and Carbon Cycle Cruise 1, GOMECC1, https://www.aoml.noaa.gov/ocd/gcc/GOMECC1/), 2012(Gulf of Mexico Ecosystem and Carbon Cycle Cruise 2, GOMECC2, https://www.aoml.noaa.gov/ocd/gcc/GOMECC2/), 2015 (East Coast Ocean Acidification Cruise 1, ECOA1), and 2018 (East Coast Ocean Acidification Cruise 2, ECOA2, https://www.aoml.noaa.gov/ocd/gcc/ECOA/) (fig. S1 and table S1). Except for the 1996 cruise, we are part of the team that collected the data under the auspice of the Ocean Acidification Program (OAP), National Oceanic and Atmospheric Administration (NOAA).
Samples for dissolved inorganic carbon (DIC), total alkalinity (TA), oxygen (O2), and dissolved inorganic nutrients (nitrate and nitrite (N), phosphate (P), silicate (Si)) were sampled from Niskin bottles. Salinity, temperature, and pressure data were from conductivity, temperature, and depth (CTD) sensors. In this work, DIC samples were measured by coulometric titration and TA was measured by acidimetric titration using the open-cell titration method. DIC and TA were checked or calibrated daily against Certificated Reference materials from Dr. Andrew Dickson's laboratory at Scripps Institution of Oceanography. The overall uncertainties of DIC and TA measurements were within 0.1%. The precision and accuracy details on all physical and chemical parameters sampled can be found in each cruise report and have been archived (https://www.ncei.noaa.gov/access/ocean-carbon-acidification-data-system-portal/).
We assessed the cruise-to-cruise measurement biases by comparing crossover deep stations (table S2) and by comparing the hydrographic parameters from the same water mass (fig. S5, S6). Since previous studies have proven that the GOMECC and ECOA cruise are consistent and comparable to the GLODAP data (69,70), we especially compared the deep stations from the OMP cruises with the GLODAP stations (depth >1000 m) within 1° × 1° grid and 200 m depth differences. The hydrographic data differences between OMP and GLODAP are similar to the comparison of deep stations between Coastal Ocean Data Analysis Production in North America (CODAP-NA) and GLODAP data. In the deep-water masses (AAIW+ and uNADW), we found no statistical differences except silicate concentrations. Sensitivity analysis shows that the differences in these hydrographic parameters have little influence on ΔCanthro and thus do not undermine our conclusions (see Method Section 2.1).
A critical concept of interpreting present ocean interior DIC storage is to divide it into natural and anthropogenic carbon components (71). The natural component is the ocean DIC pool (Cnatural) that exists in preindustrial time while the anthropogenic component (Canthro) represents the perturbation by human activities (65). For a DIC content difference (ΔDIC) between two times (t1 and t2), it includes the changes caused by the anthropogenic carbon increase (ΔCanthro) alone and natural variability (ΔCnatural) (43). Thus ΔCanthro is defined as:
Both ΔCanthro and ΔCnatural can be further separated into steady-state and non-steady-state components (72). The steady-state components represent the change of DIC in the absence of climate variability(73) while the non-steady-state components emerges from varying climate yet atmospheric CO2 is unchanged (65). Following this definition, the natural steady-state variation is assumed to be zero when a sufficiently long time is given (73) whereas the natural non-steady-state variation contributes the most of ΔC natural(t 2 -t 1 ) .
The steady-state ΔCanthro is the direct ocean response of atmospheric CO2 increase only while the non-steady-state is anthropogenic CO2 modified by other the climate changes (See details in Supplemental Material session 1).
Most methods of Canthro estimation only constrain the steady state component (72).
However, ΔCanthro from extended multiple linear regression (eMLR) method cannot be easily defined. This eMLR method is used in the subsurface and may capture some non-
steady-state signals based on a study of model-generated synthetic data.(73) Notably, ΔCanthro from eMLR depends on some important assumptions, e.g. constant coefficients of natural components, small hydrographic variability and etc., (see the discussion by analytical analysis on decomposing the coefficients of eMLR in Supplemental Material Section 1 and 2). We validate these assumptions with observation data during our application of eMLR (see details in Method 2.1). The surface ΔCanthro calculation is more challenging in the coastal regions as strong seasonality and air-sea gas exchange make the establishment of eMLR in the surface mixed layer difficult and the results not reliable. The observed ∆pCO2 climatology change can only capture the overall DIC pool changes(74). The widely adopted constant DIC disequilibrium (∆TCO2) method and constant ΔpCO2 method(3) can estimate the steady-state component but not include nonsteady-state component (see detailed discussions in Method 2.2 and Supplemental Material Section 3 and 4). All Canthro methods have caveats, especially if they are applied in the coastal ocean. We expect that non-steady state component of ΔCanthro would play a more important role in the coastal region relative to the open ocean, which can be influenced by the changes of nutrient concentration, biological process shifts, terrestrial organic and inorganic carbon inputs, chemical weathering under anthropogenic stress. We discussed the influences of these factors on subsurface and surface ΔCanthro estimation in Supplemental Material
Section 2 and Section 4, respectively. Therefore, we recommend careful evaluations of the methodology before using any Canthro estimations in coastal regions.
We determined the subsurface ΔCanthro increase rate using a time series analysis on the ΔCanthro calculated from the eMLR method with interaction term in the coastal region (regional eMLR) (42). The eMLR method has been broadly applied to the repeated transects in the open ocean (29,30,(75)(76)(77). However, such a method has not been used to constrain ΔCanthro in coastal areas due to limited long-term observations and complex natural processes that may mask part of the anthropogenic signals and lead to large uncertainties in the ΔCanthro estimates. Recent papers in the coastal regions have successfully built regional multiple linear regression (MLR) with interaction terms that can reconstruct the carbonate parameters with a high coefficient of determination (R 2 ) (78)(79)(80)(81)(82), which provides a solid foundation to apply.
The first step of the regional eMLR method is to build regional MLR models separately for five different cruises from 1996 to 2018. Here, we built MLR models following the steps in Li et al (70). To avoid the collinearity among the predictor variables and their interaction terms, all terms that lead to variance inflation factors (VIF) > 10 were excluded from the model. The combinations of T, S, O2 and Si were chosen as the predictor parameters to build the MLR function because they deliver the smallest VIF values. We excluded water depth < 20 m since the summertime mixed layer depth in the U.S. East Coast is less than 20 m(83). The lower limit was set at 2000 m because the residuals began to shift with a trend in the deep layer (Fig. S9) and model performance deteriorated. To perform the eMLR method among different cruises, all models in different cruises must have the same predictor terms in the model for all cruises and O2×Si is the only interaction term. Therefore, the final model equation for all cruises is,
X (𝑡 𝑖 ) represents the anomalies (measured values subtract the mean) of input matrix of hydrographic data (i.e., salinity, temperature, O2, and Si) and C (𝑡 𝑖 ) is the matrix of intercepts and coefficients (table S3),
1,(𝑡 𝑖 ) 𝛽 2,(𝑡 𝑖 ) 𝛽 3,(𝑡 𝑖 ) 𝛽 4,(𝑡 𝑖 ) 𝛾 (𝑡 𝑖 ) ] 𝑇
The final regression equation is,
where Xm values are the mean of each data group (table S4).
To perform eMLR, we set 1996 as the reference year (t1). Thus, ΔCanthro indicates the amount of anthropogenic carbon accumulated between any year (ti) and the reference year (t1), There are a few assumptions in the application of the eMLR method, e.g., no temporal changes of relationship between predictive and predicted parameters in the natural components, small differences in the residuals between any two MLRs, and limited effects of hydrographic parameter variation on ΔCanthro estimation. We examined these assumptions of eMLR and re-confirmed their validity in our study areas by conducting uncertainty analyses (see more discussion in Supplemental Material, sections 1 and 2).
We especially tested the influence of hydrographic parameter changes on anthropogenic carbon by setting 2018 (instead of 1996) as the reference year and calculated ΔCanthro increase from 1996 to 2018 (Fig. S10). The ΔCanthro increase in different water masses shows similar results as those from setting 1996 as the reference year (Fig. 2).
We quantitatively estimated the uncertainties in the ΔCanthro that arose from the MLR and eMLR models. The uncertainties in the MLR models come from the statistical model, the measurement uncertainty of response carbonate parameters and inputs (73,(84)(85)(86). We combined these errors as the square root of the sum to evaluate the overall uncertainties (Eest) (70).
Here, Eresponse represents DIC measurement uncertainties in the whole data set (0.1%).
The Einput is the uncertainty of the input parameters propagated in the regression, which can be calculated as √∑ (𝑈 𝑖 𝛽 𝑖 ) 2 𝑛 𝑖=1
. Ui is the i-th input uncertainties of the predictor properties. We set the Ui values 0.005, 0.005°C, 1% and 2% for S, T, O2, and Si respectively (87,88). The 𝛽 𝑖 terms are the i-th regression coefficients in the MLR models.
EMLR refers to the model root mean square error (table S3). In addition, we conducted a sensitivity analysis on how much the change of hydrographic data from time one to time two affects the ΔCanthro. We used the mean differences of hydrographic parameters in AAIW+ and uNADW water among the five cruises as input sensitivity tests (Figure S5).
Besides, the sea surface T and S change in the MAB is about 0.2°C and 0.2 S(89, 90). Our results suggest that changes of 0.2°C, 0.2 S, 10 µmol kg -1 of O2, and 3 µmol kg -1 of Si correspond to a ΔCanthro change of 1.6 ± 0.5 µmol kg -1 , 2.4 ± 0.4 µmol kg -1 , 3.0 ± 0.4 µmol kg -1 and 2.5± 1.9 µmol kg -1 , respectively.
We also tested the model robustness by substituting Si with N or P in the combination of predictors in the eMLR (86). We evaluated the mean ΔCanthro (ensemble MLR method (30) )
and compared the ΔCanthro from different models to estimate the uncertainty ranges caused by the choice of different parameters (table S5 and Fig. S11). The ensemble eMLR also derived similar ΔCanthro increase rates as shown in eMLR approach (table S6).
Finally, to validate the overall accuracy of our algorithm, we applied our regional eMLR method to a dataset with known Canth additions to test if it can capture the low Canthro especially in the low salinity waters (Fig. S12).
The surface ΔCanthro calculation is more challenging in the coastal regions as strong seasonality and air-sea gas exchange make the establishment of eMLR in the surface mixed layer difficult. In the previous studies, several approaches are used for calculating surface ΔCanthro: 1) Extend the eMLR coefficients derived from the subsurface layer to the surface(91); 2) Constant DIC disequilibrium (ΔTCO2) method( 31), and 3) Constant
ΔpCO2 method(3). In this paper, we adopted the constant ΔpCO2 method because it takes the seawater initial pCO2 and buffer capacity into consideration (31), especially in the nearshore regions where pCO2 is high and mainly driven by biological fluxes (92,93). We also examined and compared with the surface ΔCanthro derived from other methods (in the session and in Supplementary Material section 3 and 4).
For convenience, we start with known pCO2(t2), DIC(t2), and TA(t2) and estimate the CO2 state of in the past (t1). For the constant ΔpCO2 method, we assumed that the gradient between pCO2sw and pCO2air is constant during the study period (Δ𝑝CO 2(𝑡 1 ) = Δ𝑝CO 2(𝑡 2 ) ) under instantaneous equilibrium assumption. Thus, the estimated seawater pCO2 at time t1 is,
Here, pCO2(t1), sw, eq is the surface seawater pCO2 in equilibrium with the atmosphere at time t1. The estimated DIC at time t1 is calculated as,
Then, ΔCanthro in the surface layer is calculated as,
where DIC (𝑡 2 ),𝑜𝑏𝑠 is the observed DIC at time t2. The constant ΔpCO2 assumption may not be validate for the coastal regions as suggested by previous studies (16,17) that the coastal regions, e.g., MAB, are a growing CO2 sink. However, the 22-year time span is not long enough to detect the surface anthropogenic pCO2sw signal changes based on the time of trend emergence analysis (Supplemental Materials, Section 3).We are aware that the constant ΔpCO2 method is not perfect for mixed layer ΔCanthro estimation and highly influenced by the seasonality, yet it provides a credible range of ΔCanthro in the surface layer and the uncertainty is assessed to be acceptable (Supplementary Material section 4).
In addition, we compared the surface ΔCanthro estimations to two other approaches: 1)
extending eMLR coefficients from subsurface to surface, and 2) constant ΔTCO2 method.
At high salinity end (S>35), all methods gave similar ΔCanthro as the constant ΔpCO2 method. At low salinity (S<32), we got smaller ΔCanthro values when 1) extending eMLR coefficients from the subsurface (Fig. S13) and larger ΔCanthro values when 2) applying constant ΔTCO2 method (Fig. S13). The results of the constant ΔpCO2 method fall within the range of the ΔTCO2 method and that by extending eMLR coefficients to surface waters. Results from ΔTCO2 delivered relative conservative results as it does not take seawater pCO2 into consideration (31,32). Note that if we use the constant TCO2 approach, the air-sea pCO2 gradient and flux will become greater with time, whereas airsea pCO2 gradient and flux remain the same in the constant pCO2 (see Methods).
Nonetheless, all methods show a decreased ΔCanthro from the open ocean towards the nearshore. Further details and estimations of surface ΔCanthro from different methods can be found in the Supplemental Material section 4.
To quantify the Canthro flux across the land-ocean interface, we calculated the coastal Canthro accumulation inventory and air-sea gas total flux using observation data. The airsea CO2 flux of the MAB is calculated as follows (94):
where k = 0.251 × U10 × (Sc/660) -0.5 is the gas transfer velocity calculated from wind speed at 10 m (U10) from ERA5 monthly averaged data on the sea surface (0.25° × 0.25° gridded, https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5) and the Schmidt number (Sc), K0 is the CO2 gas solubility (95). pCO2air is calculated using CO2 S7).
We set the simulation time step as one day. For each time step, DIC concentration in the coastal box was calculated as, DIC 𝑡+1 = DIC 𝑡 + ∆DIC 𝑡,(𝑝ℎ𝑦) + ∆DIC 𝑡,(𝑏𝑖𝑜) + ∆DIC 𝑡,(𝑔𝑎𝑠) (10) where ∆DIC 𝑡,(𝑝ℎ𝑦) , ∆DIC 𝑡(𝑏𝑖𝑜) , and ∆DIC 𝑡,(𝑔𝑎𝑠) represent the DIC changes induced by physical mixing of riverine and marine end-members, net photosynthesis and respiration, and gas exchange, respectively. ∆DIC 𝑖,(𝑏𝑖𝑜) and ∆DIC 𝑖,(𝑔𝑎𝑠) was calculated as, ∆DIC 𝑡,(𝑏𝑖𝑜) + ∆DIC 𝑡,(𝑔𝑎𝑠) = (NCP + 𝐹CO2, t) / MLD (11) where the term FCO2, t is calculated following Equation ( 9
The DICmix,t is calculated as
We set a series of different fractions of estuarine and marine end-member combinations (from 𝑓 𝑚𝑎𝑟𝑖𝑛𝑒 = 0.1 to 𝑓 𝑚𝑎𝑟𝑖𝑛𝑒 = 0.9 ) in the simulation (Fig. S8). S7) as the estuarine end-member (𝐷𝐼𝐶 𝑒𝑠𝑡𝑢𝑎𝑟𝑖𝑛𝑒 ) and presumed it to be unchanged and not influenced by the atmospheric CO2 increase during the time of our concern due to low water buffer capacity and a small CO2 outgassing. With the new DICt+1, a new pCO2, t+1 was calculated with constant TA, T, S, and nutrients in CO2SYS
(MATLAB version( 97)). The new pCO2, t+1 is substituted into the simulation iteration to calculate ∆DIC 𝑡,(𝑔𝑎𝑠) in the next time step.
In addition, we ran another two sets of simulations with different estuarine pCO2 values:
an extremely high sea surface pCO2 scenario (700 ppm) representing small estuaries like Hudson River and a relatively low pCO2 scenario (400 ppm) representing the large estuaries like Chesapeake Bay in the MAB (Fig. S14 and Table S7). All simulation results indicated that fast mixing processes and short residence time restrains sea surface pCO2 increase in the coastal box at a rate slower than both the atmospheric CO2 and open ocean box pCO2 increase (Fig. S8 and Fig. S14).
Commented [LX1]: Our references are 18 more than required. (Main text reference 1:98, 99-end for the supplements) Reference list should be a reasonable accounting of the critical literature supporting the conclusions of the paper. Ensure inclusion of a diverse set of authors. Do not include excess references. The suggested upper limit for research articles is 80 and for review articles 150. Al Commented [LX2R1]: Well It is suggested and I go over the paper again and found nothing can be taken out In the maintext, we have 67 citations. 22 in methods. Fig. 1. Dissolved inorganic carbon (DIC) spatial distributions and decadal changes. The depth profiles of DIC in the most southern transect (A) Chesapeake Bay (CB) and the most northern (B) Massachusetts (MA) transect were exactly repeated in all five cruises from 1996 to 2018. DIC content increases from nearshore to offshore and increases with depth. The temperature and salinity diagram (C) reveals the major water masses with distinct DIC content. On the shelf, the nearshore The DIC content in the water column increases with time (D), especially in the off-slope top 200 m of Gulf Stream (GS) water (fig. S1). S6). (2022) and the yellow box indicates our research area, the continental shelf water, while the light blue (left box) and dark blue (right box) indicate the estuaries-wetland box and open ocean box, respectively. The numbers in black are present carbon flux from Regnier et al., (2022). The unit is PgC yr -1 . The red number or flux (F) indicates anthropogenic perturbations. FEC represents the anthropogenic flux from estuaries to continental shelf; FAC from atmosphere to continental shelf; FCO from continental shelf to ocean. FCS is the Canthro accumulated in the continental shelf box. (b) The Canthro flux (F) and export percentage of the continental shelf box derived from this research in the MAB (red numbers) and from other models or synthesis in global scale (blue numbers).
Here we further divided FCS into the Canthro accumulated in the water column (FCS_W) and buried in sediments (FCS_S).
* Observational results from this study. ** The Canthro air-sea gas flux is assumed to be ~ 30% -70% of the total air-sea flux (1.6± 0.3 Tg C yr -1 ) (15,18,20,54,55) #