Carbon cycling on the East Siberian Arctic Shelf – a change in air-sea CO 2 flux induced by mineralization of terrestrial organic carbon

Measurements from the SWERUS-C3 and ISSS-08 Arctic expeditions were used to calibrate and validate a new physical-biogeochemical model developed to quantify key carbon cycling processes on the East Siberian Arctic Shelf (ESAS). The model was used in a series of experimental simulations with the specific aim to investigate the pathways of terrestrial 20 dissolved and particulate organic carbon (DOCter and POCter) supplied to the shelf. Rivers supply on average 8.5 Tg C yr dissolved inorganic carbon (DIC), and further 8.5 and 1.1 Tg C yr DOCter and POCter respectively. Based on observed and simulated DOC concentrations and stable isotope values (δCDOC) in shelf waters, we estimate that only some 20% of the riverine DOCter is labile. According to our model results, an additional supply of approximately 14 Tg C yr eroded labile POCter is however required to describe the observed stable isotope values of DIC (δCDIC). Degradation of riverine DOCter 25 and POCter results in a 1.8 Tg C yr reduction in the uptake of atmospheric CO2, while degradation of eroded POCter results in an additional 10 Tg C yr reduction. Our calculations indicate nevertheless that the ESAS is an overall small net sink for atmospheric CO2 (1.7 Tg C yr). The external carbon sources are largely compensated by a net export from the shelf to the Arctic Ocean (31 Tg C yr), and to a smaller degree by a permanent burial in the sediments (2.7 Tg C yr).

1 Introduction 30 The total soil carbon pool in Arctic permafrost regions amounts to approximately 1700 Pg C, which is about half of the estimated global pool (Tarnocai et al., 2009).This carbon pool exceeds the atmospheric carbon pool (~860 Pg C at 400 ppm) by a factor two. Permafrost thawing and coastline erosion are potentially large sources of terrestrial organic carbon (OCter = DOCter and/or POCter) to the East Siberian Arctic Sea (ESAS).One of the key questions concerning carbon cycling in in this areaand globallyis the degradability of OCter previously stored in frozen soils.A high mineralization of OCter could 35 significantly enhance CO2 outgassing, although that signal at least to some extent may be compensated by an increased plant growth (Schuur et al., 2009).Enhanced outgassing would then increase the atmospheric CO2 level and capacity to trap heat, which could in turn accelerate permafrost thawing and potentially generate a positive feedback from permafrost thawing to global warming.
The development of the permafrost carbon pool as well as the different processes behind the removal of carbon from 40 the pool and outgassing to the atmosphere have been discussed in detail by Schuur et al. (2008).Further, the Arctic carbon cycle and the sensibility to climate change for some of the key processes were reviewed by McGuire et al. (2009).While the Biogeosciences Discuss., doi:10.5194/bg-2017-115, 2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.fate of OCter previously stored in permafrost is not clear, model calculations indicate that possibly up to 1 Pg C yr -1 (corresponding to an increase in the atmospheric CO2 level by almost 0.5 ppm yr -1 ) might be outgassed to the atmosphere as a result of carbon mobilization by permafrost thawing (Schuur et al., 2008).Schaefer et al. (2011) estimated that up to 200 Pg C (corresponding to ~90 ppm atmospheric CO2) could be released from thawing permafrost by the year 2200.The total increase from present-day atmospheric CO2 levels has for comparison been estimated to end up in a range 100 to 1000 ppm by the year 5 2100 depending on the development of several different anthropogenic processes (Moss et al., 2010).
Organic carbon of marine origin (OCmar) is usually distinguished from OCter by stable isotope composition.The δ 13 C value of OCmar is typically found to be in the range -23‰ to -21‰ (Deutsch et al., 2012), although this value is heavily dependent on the δ 13 CDIC value of DIC assimilated by phytoplankton.For example, δ 13 C values in a range -32‰ to -29‰ have been reported for OCmar in ESAS waters, and these depleted values are believed to be the result of assimilation of DIC where 10 the δ 13 CDIC values were low because of OCter degradation (Alling et al., 2012).OCter is nevertheless usually more depleted in 13 C than OCmar; discharge-weighted δ 13 C values for DOCter and POCter in Lena are on average -27‰ and -29‰ respectively (Raymond et al., 2007;McClelland et al., 2016).The isotopic composition of DOC in coastal marine systems is partly determined by the relative influence of riverine carbon loads and marine production, and further by the fate (export, mineralization, burial) of OCter and OCmar in the system.15 In the upper layers (<500 m) of the Arctic Ocean, δ 13 CDIC values have been observed in a range of -0.6 to +2.2‰ (Bauch et al., 2015), while river water is typically much more depleted in 13 C (δ 13 CDIC ≈ -8‰ in Lena; Alling et al. 2012).In the ESAS there are clear gradients of e.g.surface water δ 13 CCO2 values along the coast, with lower values in the west than in the east (SWERUS-C3 cruise report, 2016)a response to a more or less strong terrestrial influence.In addition to the effects of conservative mixing of end-members, stable isotope values in shelf waters are strongly influenced by several other 20 processes.Ecosystem transformation processes such as e.g.carbon assimilation by primary producers, mineralization in the water column and sediments, and outgassing/uptake of CO2 leave characteristic isotopic fingerprints in the various carbon reservoirs of a marine system.
DOC can be transformed into DIC by photomineralization and microbial degradation.The bioavailability or degradability of DOCter mobilized through permafrost thawing differ depending on origin (Schuur et al., 2008).Bioavailability 25 is in addition not exclusively regulated by the lability of the organic substrate, but can further depend on concentration itself.
Even labile DOC may escape degradation by being too dilute to be consumed by microorganisms (Middelburg, 2015).Further, the microbes that degrade DOC may be restricted by other factors, such as e.g.nutrient and/or temperature limitation or predation (Thingstad et al., 1997;Zweifel, 1999).Kuliński et al. (2016) found that less than 20% of the riverine DOC in the subarctic Baltic Sea is mineralized, even 30 after months of incubation.This mainly refractory behavior of riverine DOC has earlier been reported by Hoikkala et al. (2015) and Dittmar and Kattner (2003).In contrast, Vonk et al. (2013) demonstrated that old DOC in stream waters draining Yedoma deposits is highly bioavailable (34% loss after 14 days incubation), while mixtures of Yedoma stream DOC with river and ocean waters generally showed somewhat lower bioavailability (18-32% loss after 14 days incubation).Alling et al. (2010) estimated an overall 30-50% removal for DOCter in ESAS waters . 35 In addition to the large DOCter reservoirs in shelf waters, a clear terrestrial influence is evident in coastal sediments of the Laptev and East Siberian Seas (LS and ESS) (Semiletov et al., 2005;2011;Vonk et al., 2014).This implies a large sedimentation of POCter supplied by rivers and/or coastline erosion.The riverine supply of POCter from Lena and Kolyma to LS and ESS has been estimated to amount to about 0.94 Tg C yr -1 combined (McClelland et al., 2016), while estimates of the OCter supply to the ESAS through coastline erosion vary by a factor twenty; from ~2 to 44 Tg C yr -1 (Grigoriev et al., 2004 40 and references therein;Semiletov et al., 2011;Vonk et al., 2012).
The main aim of this study is to estimate to what extent mineralization of riverine DOCter and eroded POCter respectively influence the exchange of CO2 between shelf waters and the atmosphere.To this end, a coupled physical- Discuss., doi:10.5194/bg-2017Discuss., doi:10.5194/bg- -115, 2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.biogeochemical model with a specific focus on carbon and stable carbon isotopes has been developed.The model is well suited for large scale budget calculations, and explicitly accounts for processes such as production and respiration of OCmar, supply and respiration of OCter, air-sea CO2 fluxes, exchange between sediments and water column, as well as transports between the shelf and the Arctic Ocean.By first adjusting DOC concentrations and isotope values in the model (as constrained by observations), we can then estimate the "residual" OCter mineralization required to reproduce e.g.observed DIC concentrations 5 and δ 13 CDIC values.This residual OCter mineralization must be primarily coupled to a supply of labile OCter from coastline erosion, since the supply of riverine POCter is comparatively small (cf.Section 2.3.2).We will further quantify the key carbon fluxes from land to sea, the net air-sea CO2 exchange, as well as the export of carbon from the shelf to the Arctic Ocean.

Study area 10
In this study we defined the ESAS area (or rather the model domain) as limited at 110ºE in the west and 180ºE in the east.To the north, the model domain is limited by the 100 meter isobaths.In our model setup we divided this area into three connected sub-basins (Fig. 1).LS is divided into one inner and one outer sub-basin (ILS and OLS); ILS is limited at 125ºE in the west, 140ºE in the east, and 75ºN in the north.ESS is limited at 140ºE in the west and 180ºE in the east.Bathymetric properties for these three basins are summarized in Table 1.Altogether the study area covers about 1300000 km 2 , the average depth is 34 m, 15 and the volume is 44300 km 3 .This section of the shelf is approximately 2000 km wide in the along-shore direction.In addition to the overarching background control exerted by boundary conditions, the great Siberian rivers that drain the ESAS catchment area have a defining influence on the hydrology of the shelf.With an annual freshwater supply of 490 km 3 , the Lena River dwarfs all other rivers and accounts for ~60% of the total runoff (Table S1, supplementary material).More than one third of the total runoff is supplied in June.The impact of the Lena River is visible throughout LS and further 25 north/north-east of the New Siberian Islands well inside ESS.The riverine signal is apparent not only in distributions of salinity and temperature, but also in e.g.DOC, DIC, and nutrient concentrations, as well as stable carbon isotope values in the water column and sediments (Semiletov et al., 2005;2016;Alling et al., 2010;2012;SWERUS-C3 cruise report, 2016).
Under present-day climate and weather conditions, ILS and OLS are usually altogether ice-free by late July/early August.The minimum ice extent of ESS occurs somewhat later; typically in early September.By mid-November, all three 30 basins are again mostly covered by ice (Fig. S1, supplementary material).

The model
To study carbon cycling in the ESAS we have developed a new physical-biogeochemical model.Each of the three model subbasins (Fig. 1) is described as horizontally homogeneous although with a high vertical resolution and a depth-dependent area distribution.The model parameterizations were adopted from a physical-biogeochemical Baltic Sea model (BALTSEM) that 35 has been described and validated in numerous publications: the model includes a hydrodynamic module (Gustafsson, 2000;2003) that has subsequently been expanded to include nutrient and plankton dynamics (Gustafsson et al., 2012;Savchuk et al., 2012), organic carbon and the carbonate system (Gustafsson et al., 2014a), and stable carbon isotopes (Gustafsson et al., 2015).
The model is described in some detail in Appendix A-C.All pelagic and benthic state variables included in the model are presented in Table B1 and B2  The model simulations in this study covered the period 1979-2014.All model runs were however first spun up over a period of 70 years to ensure that the output is independent from initial conditions.

Atmospheric forcing
Atmospheric forcing from January 1979 to August 2015 was obtained from the ERA-interim 0.75 x 0.75 gridded reanalysis 5 (Dee et al., 2011).The forcing parameters are: 10 m zonal/meridional wind velocity components, total cloud cover area fraction, 2 m air temperature, mean sea level pressure, total precipitation, and relative humidity.A description of the calculations of relative humidity is available in Text S1 (supplementary material).72.75ºN),(129.75ºE,75.75ºN)and (159.75ºE,75ºN) grid cells were chosen to represent model boxes 1, 2 and 3 respectively (Fig. 1).The ERA-interim forcing conditions are updated four times a day. 10

River loads and atmospheric data
River runoff data were collected from two sources; i) discharge data from before 1999 were extracted from the R-ArcticNET database (http://www.r-arcticnet.sr.unh.edu/v4.0/index.html),and ii) discharge data from 1999 to present as well as river loads of nutrients, total alkalinity (AT), inorganic and organic carbon were extracted from the PARTNERS and Arctic-GRO databases (http://arcticgreatrivers.org/data.html).15 Average river loads of DOC and POC have been estimated to 5.6 and 0.81 Tg C yr -1 respectively from Lena, and 0.82 and 0.12 Tg C yr -1 respectively from Kolyma (Holmes et al., 2012;McClelland et al., 2016).Lena and Kolyma represent 75% and 63% of the total runoff to LS and ESS respectively (Table S1, supplementary material).In this study we assumed that the other rivers entering LS and ESS share properties in terms of carbon and nutrient concentrations with Lena and Kolyma respectively.This means that the simulated river loads of dissolved and particulate constituents to LS and ESS are 33% and 20 60% larger respectively than the loads from Lena and Kolyma alone.The stable isotope values for riverine carbon were set to -8‰ for DIC (Alling et al., 2012), -27‰ for DOC (Raymond et al., 2007), and -29‰ for POC (McClelland et al., 2016), and were for simplicity assumed to be constant.
The stable isotope values for coastal ice complex and top soil permafrost have been estimated to approximately -26.3 and -28.2‰ respectively (Vonk et al., 2012).In the sensitivity experiments below where we calibrated the supply of eroded 25 OCter (Section 2.4), it was assumed that the average isotope value for eroded OCter is -27‰ (and constant).

Atmospheric CO2 levels and isotope values were based on measurements by the Earth System Research Laboratory
Global Monitoring Division at the Barrow station in Alaska (http://www.esrl.noaa.gov/gmd/dv/data/index.php?site=brw).

Boundary conditions and initial profiles
Boundary conditions for all state variables at the four model boundaries (Fig. S2, supplementary material) were based on data 30 downloaded from the World Ocean Database (WOD; https://www.nodc.noaa.gov/OC5/WOD/datageo.html) together with data from the SWERUS-C3 expedition.

Calibration and validation data
Model calibration and validation data were based on measurements from the ISSS-08 and SWERUS-C3 expeditions, and in addition WOD data.Data from ISSS-08 and SWERUS-C3 include measurements of the carbonate system and stable carbon 35 isotopes, and also temperature, salinity, and nutrient concentrations.We downloaded additional measurements of salinity, temperature, and nutrients from WOD. Carbon and carbon isotope data were however not available in WOD.The observed profiles are mainly from July and August, but in a few cases also September.Sampling procedures and data analyses have Biogeosciences Discuss., doi:10.5194/bg-2017-115, 2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.been described in detail by Alling et al. (2010;2012) and in the SWERUS-C3 cruise report (2016).Sea ice data were provided by the National Snow and Ice Data Centre (Cavalieri et al., 1996).

Model sensitivity experiments
The model was first used to calibrate the labile fraction of DOCter supplied by rivers and atmospheric depositions, and then to estimate the degradation of OCter supplied to the sea by coastline erosion.In an additional experiment we adjusted the C:N:P 5 ratios for phytoplankton assimilation.The model optimization procedure was performed in several consecutive steps: i) In a first simulation we assumed that the DOCter supplied by rivers and atmospheric depositions was 100% refractory (Sim. 1, cf.Table 2), and in a second run that the DOCter was on the contrary 100% labile (Sim.2).In these two experiments we assumed further that the riverine POCter was 100% labile, and that the supply of eroded OCter was zero.The purpose of these two experiments was to determine the ranges of DOC and DIC concentrations and δ 13 CDOC and δ 13 CDIC values 10 that are possible to obtainand particularly if it is at all possible to reproduce observed DOC, DIC, and their isotope valuesby DOCter mineralization alone.
ii) Then, in a series of simulations we optimized the labile DOCter fraction and further the supply of labile eroded OCter.For simplicity, it was assumed that the eroded OCter was all in particulate form (POCter) although previous publications suggest that some 20% of the eroded carbon is dissolved (Sánchez-García et al., 2011 and references therein).The eroded 15 POCter was further assumed to be 100% labileor in other words, we only considered the labile part of eroded POCter.One fraction of the simulated POCter sinks to the sea floor where it is partly mineralized and partly buried.The other fraction is broken down into DOCter in the water column where it can be mineralized or transported to adjacent sub-basins and the Arctic Ocean.Thus, eroded POCter contributes to some extent also to the pelagic DOCter concentrations and δ 13 CDOC values in the model.For that reason the model was run in an iterative approach were both the labile DOCter fraction and the supply of labile 20 POCter were gradually adjusted until observed DOC and DIC concentrations as well as their isotope values were reproduced as far as possible.The best model run from this step is referred to as Sim. 4. For reference, we also kept a model run with the same labile DOCter fraction as in Sim. 4, but with zero supply of eroded POCter.This model run will be referred to as Sim. 3 (cf.Table 2).
iii) In the experiments above, Redfield ratios were used for the C, N, and P assimilation by autotrophs (i.e., 106:16:1 25 molar).In a final experiment we used the best model simulation as described above (Sim.4) but adjusted the C:N:P ratios in an attempt to further improve the model results.This final model simulation is referred to as Sim. 5.

Results
In this section we first validate the physical (Section 3.1.1)and then the biogeochemical parameters (Section 3.1.2).Budget calculations are then used to quantify carbon transformation processes and transports in the ESAS (Section 3.2-3).30

Salinity and temperature profiles
Here we compare simulated average salinity and temperature profiles from the period July-September 2000-2014 to observed values (Fig. 2).Many of the observations from primarily ILS (and to some extent OLS) that we used to validate the model output are directly within the Lena Plume, characterized by a very low salinity and high temperature in surface waters.This 35 water mass could not be reproduced within the 1D-modelling approach with horizontally integrated layers.Overall we found a smaller variability in the model output than in the observed data.The variability of the observed data is the result of horizontal gradients (mainly) and temporal changes; the variability in the model is due to temporal variations.Discuss., doi:10.5194/bg-2017Discuss., doi:10.5194/bg- -115, 2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.

Biogeochemical parameters
Observed deep water DOC concentrations and δ 13 CDOC values were reasonably well reproduced by the model assuming a labile DOCter fraction in a range 0-30% (not shown), but far from observed values assuming instead a 100% labile fraction (Sim.2).
It was not possible to reproduce observed δ 13 CDIC values by adjusting only the labile fraction of DOCter; a large supply of eroded labile POCter was consequently required to obtain reasonable δ 13 CDIC values.The best model results were obtained by 5 assuming a 20% labile DOCter fraction combined with a supply of 14 Tg C yr -1 eroded labile POCter (the Sim. 4 experiment, Table 2).By adjusting the C:N:P ratios for phytoplankton assimilation and instead of Redfield ratios using molar C:N:P ratios of 106:8:1, it was possible to further improve the model results for a few parameters (Sim.5).In Fig. 3 2) and compared to observed values.In Fig. 4, average simulated July-August profiles of simulated AOU, 10 nitrate (NO3 -), phosphate (PO4 3-), and CO2 fugacity (fCO2) are compared to observations.As mentioned above, the model is not capable of reproducing the Lena plume, and this was particularly evident in above-pycnocline DIC, DOC, δ 13 CDIC, and δ 13 CDOC values in ILS and OLS, though not in ESS that is not as strongly influenced by the Lena discharge (Fig. 3).In the simulations with no coastline erosion (Sim.1-3), δ 13 CDIC values were considerably overestimated below the pycnocline, indicating that the OC respiration was too low.The only exception was the experiment 15 with 100% labile DOCter (Sim.2); in this case δ 13 CDIC values were reasonably well reproduced in OLS.However, a 100% labile DOCter fraction is not feasible as shown by simulated DOC and δ 13 CDOC profiles.In the Sim. 4 and Sim. 5 experiments with calibrated POCter supply from coastline erosion, AOU and δ 13 CDIC values were reasonably well reproduced below the pycnocline.In the Sim. 5 experiment where C:N:P ratios for phytoplankton were adjusted the discrepancy between simulated and measured surface water PO4 3-concentrations was least pronounced (Fig. 4).20 We have no observations of δ 13 CDOC values in ILS.The model results nevertheless suggested a very strong influence of terrestrial DOC (i.e., values close to -27‰ which is the preset value for DOCter), particularly above the pycnocline whereas the below-pycnocline values were somewhat larger because of a relatively larger contribution from DOCmar.In OLS we found particularly large discrepancies between simulated and measured DOC and δ 13 CDOC values above the pycnocline.However, the few DOC measurements in the western OLS (that are not influenced by the Lena plume) were well reproduced.25 Several signs pointed towards an underestimated productivity in the simulated ESS basin: Too high fCO2 and too low δ 13 CDIC values in the surface water but the opposite signal in the deep water.Deep water AOU, NO3 -, and PO4 3-concentrations in ESS were further underestimated in all model simulations.In this sub-basin we generally found that model results improved by assuming non-Redfield ratios for the C, N, and P uptake by autotrophs (the Sim. 5 experiment).

Budget calculations 30
Key carbon fluxes from the five experiments are listed in Table 2.The net ecosystem production (NEP) in Table 2 is defined according to Woodwell and Whittaker (1968): NEP = gross primary production -respiration by autotrophs and heterotrophs.
NEP is equal to the internal pelagic DIC sink, which is exactly equal to the internal pelagic total organic carbon (TOC = DOC + POC) source.Sediment mineralization is a sink for sediment carbon that is exactly equal to a DIC source to the water column (this means that neither NEP nor sediment mineralization influences the total carbon (TC = DIC + TOC) budget).The carbon 35 budget is closed by a net accumulation of carbon in the water column and active (i.e., not permanently buried) sediment layer (net accumulation = river load + coastline erosion + atmospheric deposition + net air-sea CO2 flux -net export -permanent burial).
Comparing the Sim. 1 and Sim. 2 experiments (i.e., 100% refractory vs. 100% labile DOCter), we found that the system in both cases was an overall net sink for atmospheric CO2.In Sim. 2, the net CO2 uptake was reduced by more than 5 40 Tg C yr -1 compared Sim. 1 (9.2 vs. 3.9 Tg C yr -1 , cf.ILS was a net CO2 source to the atmosphere every year in all five simulations (not shown).OLS and ESS were some 5 years CO2 sources and some years CO2 sinks, but on average they were CO2 sinks in all experiments.In the three experiments that did not include coastline erosion (Sim.1-3), the system was an overall sink for atmospheric CO2.In the Sim. 4 experiment where we included the calibrated supply of eroded labile POCter (14 Tg C yr -1 ), the CO2 absorption decreased by 10 Tg C yr -1 compared to the Sim. 3 experiment (Sim.4 is apart from coastline erosion identical to Sim. 3).In the Sim. 4 experiment the system was on average an overall CO2 source (1.8 Tg C yr -1 ).In the Sim. 5 experiment with adjusted phytoplankton C:N:P 10 ratios, NEP more than doubled compared to Sim. 4, and further the net CO2 uptake was enhanced by 3.5 Tg C yr -1in this simulation the ESAS was thus a net CO2 sink of 1.7 Tg C yr -1 .The adjusted C:N:P further induced considerable changes in sediment mineralization, carbon export, and burial fluxes increase as well (Table 2).Main pelagic and benthic carbon pools and fluxes integrated for the entire ESAS as calculated in the Sim. 5 experiment are presented in Fig. 5.

Pathways of OCter supplied to the shelf 15
In Table 3, we present budget calculations for OCter in the ESAS according to the Sim. 5 experiment (where we overall were able to best reproduce observed biogeochemical parameters; Fig. 3-4).The calibrated load of POCter from coastline erosion exceeded the riverine TOCter load by almost 45% (14 vs. 9.6 Tg C yr -1 ).Overall, the two main sink terms for OCter were transport out of the system and mineralization in the sediment.Pelagic mineralization and permanent burial contributed further to the OCter removal.We estimated the residence time of OCtercalculated as total OCter pool (pelagic + benthic) divided by 20 sink terms (net export + burial + pelagic and benthic mineralization)to be in a range 2.0-4.6 years depending on sub-basin, and 3.4 years for ESAS overall.Table 3 further includes the net air-sea CO2 flux for each sub-basin as well as the change in this flux that was the result of OCter degradation.Our calculations implied that without degradation of OCter in the ESAS, there would be an additional 11.9 Tg C yr -1 absorption of atmospheric CO2.

Discussion 25
The resemblance between simulated and measured profiles is not sufficiently close to conclusively determine the air-sea CO2 exchange between shelf waters and atmosphere, but, what we can say is that under present river loads, coastline erosion, and productivity, the ESAS appears to be an overall small sink for atmospheric CO2.Degradation of organic carbon (terrestrial and marine) is in other words more than compensated by phytoplankton CO2 assimilation, carbon burial, and export to the Arctic Ocean.Earlier results from a coupled physical-biogeochemical column model for LS (Wåhlström et al., 2012) indicated an 30 overall undersaturation of surface water fCO2, resulting in a net absorption of atmospheric CO2.Our results imply that LS (ILS + OLS) would indeed have been a net sink for CO2 were it not for the supply of eroded labile POCter.Including coastline erosion in the calculations, ILS and OLS are together a net CO2 source while overall the ESAS is still a CO2 sink.Our results agree qualitatively with several earlier studies that have identified inner LS and western ESS as net sources of CO2 to the atmosphere, and the outer/eastern parts of ESS as a net sink for atmospheric CO2 (e.g.Semiletov et al., 2007;2013; Anderson 35 et al., 2011;Pipko et al., 2011).
Air-sea CO2 fluxes in the ESAS are heavily influenced by respiration of POCter supplied by coastline erosion, while on the other hand mineralization of DOCter supplied from rivers seems to be of minor importance.However, an increased supply of labile eroded POCter does not immediately translate into CO2 outgassing: According to our calculations, approximately 71% of the supply is compensated by CO2 outgassing (or rather reduced CO2 absorption), 23% is compensated Biogeosciences Discuss., doi:10.5194/bg-2017-115,2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.by export of both DOCter and DIC, and the remaining 6% by permanent burial.If on the other hand the river loads increase, 74% of the supply is compensated by export, 25% by air-sea exchange, and less than 1% by burial (assuming again that riverine TOCter is composed of 13% particulate and 87% dissolved OCter, and further that the particulate fraction is 100% labile and the dissolved fraction 20% labile).Thus, although a detailed future permafrost melt scenario is beyond the scope of this study, our simple sensitivity experiments indicate that coastline erosion is probably going to be a more important process behind 5 future CO2 emissions than changes in riverine DOC concentrations.
In the sensitivity experiment that produced the best results (Sim.5), we calculated an overall net uptake of 1.7 Tg C yr -1 for the entire ESAS (Table 2).However, as evident from the comparison between the Sim.average DIC excess of 2.1, 1.6, and 6.9 Tg C in ILS, OLS, and ESS respectively.In total, the DIC excess in the ESAS is thus 10.6 Tg C according to our model calculations (matching more or less our calculated reduction in CO2 absorption due to mineralization of OCter).
The overall global oceanic CO2 uptake has been estimated to some 2.6 Pg C yr -1 (corresponding to 25% of the anthropogenic emissions), while the accumulation in the atmosphere is approximately 4.5 Pg C yr -1 (Le Quéré et al., 2016).In 20 order to increase the atmospheric CO2 level by 1 ppm, about 2.2 Pg C must be added to the atmospheric carbon reservoir.In this context, fluxes of and order of 10 Tg C yr -1 are obviously rather marginal.The overall influence on atmospheric CO2 from OCter degradation in the ESAS is at least under present river loads negligible in a larger perspective.However, rising temperatures in the future are expected not only to accelerate permafrost thawing, but further to prolong the ice-free season.This could on one hand increase productivity and CO2 assimilation, but on the other hand enhance also the coastline erosion.25 As discussed above, our model sensitivity study identifies the supply of eroded labile POCter as the main process behind CO2 release to the atmosphere.Further, the ESAS covers only some 0.25% of the total global surface area.Thus, similar to other coastal and shelf seas in the world, carbon fluxes in the ESAS are generally disproportionally large in relation to the areal expanse of the region.

Model performance 30
Overall, the model performance is reasonable given the methodological limitations.One benefit of horizontally integrated models is the ability to do long-term and large-scale budget calculations, and to quantify key processes as well as model sensitivity to changes in e.g.external loads and climatic forcing.A couple of shortcomings/limitations of the 1D-modelling approach are expected and clearly evident.The horizontally integratedbasin-widesimulated values are here compared to measurements at several individual stations in a region characterized by strong horizontal gradients.35 The very conspicuous Lena signal apparent from observations in LS is in the model distributed over the entire horizontal areas of the respective sub-basin (ILS and OLS), and consequently much less pronounced.In addition, many of the ILS and OLS measurements lie more or less directly within the Lena path, which means that they may not represent the entire basins well.A few profiles from western OLS indeed indicate a much less stratified and less terrestrially influenced water column (Fig. 2-3).40 Another methodological issue is that although the measured profiles cover a large horizontal area, they are snapshots of a marine environment that is constantly changing and influenced by several processes at the same time.Large seasonal and Biogeosciences Discuss., doi:10.5194/bg-2017-115,2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.year-to-year variability is expected due to variations in sea ice coverage, the timing of ice retreat, river water pulse, and plankton blooms, changes in wind strength and diapycnal mixing, changes in advective transports of water and ice in and out of the system, etc.We do not expect the model to reproduce conditions during a few days a certain year, but rather to reasonably well reproduce seasonal patterns.
On the other hand, one of the major advantages using a multi-year model is that the snapshots obtained from 5 measurement expeditions can be put into a context.This is particularly important in comparatively remote and harsh regions where measurements are more or less limited to a few months per year.In seasonally ice-covered areas, factors like the timing of ice retreat and main river pulse of melt water in spring control to a very large degree both biotic and abiotic parameters on the shelf.For example, Janout et al. (2016) showed that the timing of summer phytoplankton blooms in LS can differ by approximately two months depending on ice retreat.10 In Fig. 6 we show the large seasonal variability of a number of simulated parameters.The DOC concentration in ILS surface water, for example, increases from around 200 to almost 500 µmol kg -1 when the river pulse is released over the shelf.
The DIC concentration in contrast simultaneously declines from some 2200 to about 1200 µmol kg -1 , and salinity declines from 29 to 13.The short-term impact of the river pulse then declines the farther we move from the Lena delta.DOC and DIC concentrations as well as their isotopic compositions are further largely influenced by the timing of phytoplankton blooms in 15 the model; DIC decreases as a direct response to phytoplankton CO2 assimilation, while DOC increases because of phytoplankton exudation and zooplankton excretion.During phytoplankton blooms, δ 13 CDIC and δ 13 CDOC values both increase; δ 13 CDIC values because of the preferential 12 C uptake by phytoplankton, and δ 13 CDOC values because of the production of DOCmar that is less depleted in 13 C than DOCter.If the timing of ice-retreat, river pulse, and phytoplankton blooms in the model deviates from events in the real ESAS, inconsistencies between modelled and measured profiles are inevitable.20 The model simulates ice of several different thicknesses.In Fig. 7, simulated 2 and 20 cm ice is compared to observed ice concentrations in the years 2010-2014.The ice-free season is represented rather well by the model in both ILS and OLS (especially considering the 1D-modelling approach), although the observed ice retreat is more gradual than in the model.This means also that the CO2 accumulated under the ice in winter (cf.Fig. 6d) is exhaled to the atmosphere more gradually in the real LS than in the modelled counterpart.In ESS, the model overestimates the ice retreat, which suggests that the modelled 25 CO2 fluxes in summer (in both directions) might be overestimated in this sub-basin.

POCter and DOCter sources and sinks
Mineralization of organic material is calculated as a temperature dependent rate in the model.In reality, organic material is mineralized by a combination of bacterial degradation and phototransformations; labile DOC can be photomineralized into DIC, refractory DOC can be phototransformed into labile DOC, or labile DOC can be bleached by sunlight and thus become 30 less available for bacterial degradation (Deutsch et al., 2012).The calibrated DOC mineralization in the model can be regarded as the net effect or effective mineralization resulting from all of these processes.Our calibrated 20% labile DOCter fraction is lower than the 30-50% estimate by Alling et al. (2010), but on the other hand somewhat higher than the labile fraction observed for DOCter in subarctic rivers (Hoikkala et al., 2015;Kuliński et al., 2016), and within the range found by Vonk et al. (2013) for Yedoma stream DOC mixed with river and ocean waters (18-32%).35 Removal of DOCter on the shelf depends not only on the assumed labile and refractory fractions, but also on the residence time on the shelf.Further, the spatial distribution is not only determined by advection between the sub-basins, but could also be significantly altered by flocculation and sedimentation.Flocculation of DOCter is probably largely occurring in the vicinity of river mouths, where fresh water meets sea water (Mulholland, 1981).Aggregation of DOCter into particles could potentially increase the bioavailability, because particles can be ingested by filter feeders and also colonized by bacteria -40 condensation could turn OC into a more viable energy source for bacteria.Based on our model calculations, DOCter removal by flocculation and sedimentation is probably not an important process in the overall carbon cycling on the shelf.Implementing In such an open-ended problem with several unknown parameters and process rates it is possible to obtain similar results through different combinations of rates and fluxes.For that reason, we made the assumptions that a) the riverine POCter is 100% labile, b) the eroded carbon is 100% particulate, and further c) all of the eroded POCter is labile.Thus, in our calculations we have not included the supply and burial of refractory POCter, neither from rivers nor erosion (although a fraction of the labile POCter is buried, as evident from Table 2).Refractory POCter does not contribute to either 10 DIC or DOC concentrations, stable isotope values, AOU, or air-sea CO2 fluxes, and for that reason it is not meaningful to incorporate this carbon fraction in the model at this point.
We could have divided the riverine POCter into labile and refractory fractions; the 100% bioavailability assumed in the sensitivity experiments is of course unrealistic.However, if we had reduced the supply or labile riverine POCter, then it would have been necessary to increase the calibrated supply of eroded POCter by the same amount in order to achieve the same 15 results in terms of DIC and DOC concentrations and their isotopic compositions.

Phytoplankton cell stoichiometry
In the final experiment (Sim.5), we lowered the N demand by phytoplankton (compared to Redfield ratios) in order to improve the model performance in terms of e.g.PO4 3-, AOU, and δ 13 CDIC values.Numerous publications demonstrate the very flexible C:N:P cell stoichiometry for phytoplankton, particularly during nutrient stress, low temperatures, and poor light conditions 20 (e.g. Flynn, 2003;2010;Kreus et al., 2015;Spilling et al., 2015).In the present model setup we are however restricted to constant ratios.The 106:8:1 molar ratio implemented in Sim. 5 seems to have resulted in a somewhat overestimated OC mineralization in OLS deep waters, whereas in ESS, productivity appears to still be underestimated as evident from many state variables (Section 3.1.2).Sakshaug (2004) estimated the total particulate primary production (i.e., excluding DOC production) in LS and ESS 25 to be in a range 25-40 g C m -2 yr -1 .Using the areas in Table 1, those estimates would correspond to 11-17 and 22-35 Tg C yr - 1 in LS and ESS respectively.Simulated phytoplankton assimilate CO2 to produce POCmar, and there is in addition a 25% excessive CO2 assimilation used to produce extracellular DOCmar.According to our calculations in the Sim. 5 experiment, the POCmar production or particulate primary production is on average 8.1 Tg C yr -1 in LS (2.6 and 5.5 Tg C yr -1 in ILS and OLS respectively) and 16.9 Tg C yr -1 in ESSi.e., some 50-75% of the rates suggested by Sakshaug (2004).In all, the simulated 30 POCmar production is thus on average 25 Tg C yr -1 , and there is in addition a 6.2 Tg C yr -1 production of extracellular DOCmar (cf.Fig. 5).
Further, based on calculated fCO2 values, Anderson et al. (2011) estimated that primary production in the northeastern parts of ESS (approximately 50% of the basin) is in a range 4 ± 10 Tg C yr -1 .Estimates based on fCO2 corresponds to our calculated NEP (Section 3.2) which on average is 6.0 Tg C yr -1 in LS (1.5 and 4.4 Tg C yr -1 in ILS and OLS respectively) and 35 13.1 Tg C yr -1 in ESS.
The CO2 assimilation by phytoplankton is a comparatively large source of uncertainty in the model; by changing from the Redfield ratios in Sim. 4 to the adjusted ratios in Sim. 5, the ESAS turned from a small CO2 source to a small CO2 sink (Table 2).In both experiments however, the overall picture is that LS is on average a net CO2 source and ESS a net sink on annual basis.This uncertainty is a result of the rather simplistic plankton modelling approach that we use.However, in this 40 area a more complex model formulation with flexible C:N:P ratios would be difficult to calibrate and validate because of the poor coverage of observations.Biogeosciences Discuss., doi:10.5194/bg-2017-115,2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.

Sediment properties
The active sediment pool of OCter exceeds the pelagic pool almost by a factor two according to our Sim.5 experiment (Table 3, Fig. 5).As discussed above, it is assumed that the eroded POCter is 100% bioavailable (Section 4.2).If we instead had assumed that some fraction of this carbon were impossible to degrade, then the burial loss would have increased while both the outgassing and net export would have been reduced correspondingly (and we would need a larger total supply than 14 Tg 5 C yr -1 ).Further, the δ 13 C value for organic carbon in the active sediment layer would decrease because of the relative increase in OCter compared to OCmar.In the Sim. 5 experiment, the isotopic composition of organic carbon stored in the active sediment layer is around -27‰ and -25‰ in ILS and ESS respectively (Fig. 8a); roughly corresponding to observed values (Semiletov et al., 2005).In Fig. 8b we show in addition the simulated average sedimentary distribution of terrestrial/marine OC in the three sub-basins.ILS sediments are according to our calculations largely dominated by OCter (as inferred from both isotope 10 values and terrestrial vs. marine OC fractions), while sediments in both OLS and ESS have larger contributions from marine sources.

Summary and concluding remarks
A coupled physical-biogeochemical model has been constructed with the intention to describe and quantify carbon cycling in the ESAS.A specific aim was to determine crucial transport and transformation processes of OCter supplied through river loads 15 and coastline erosion.Observations gathered during two Arctic expeditions (ISSS-08 and SWERUS-C3) have together with WOD data been used to calibrate and validate the model.A reasonable agreement between simulated and observed parameters has been achieved, and the model can thus be used as a tool for e.g.budget calculations and sensitivity studies.The main findings can be summarized as follows:  DOCter supplied by rivers appears to be largely refractory; we estimated that only some 20% is mineralized in the 20 shelf waters.This estimate agrees qualitatively with earlier publications.
 A comparatively large supply of eroded labile POCter (14 Tg C yr -1 ) is required to reproduce DOC and DIC concentrations as well as their stable isotope values.Our simulated supply of eroded POCter is approximately in the middle of earlier published estimates (~2-44 Tg C yr -1 ).
 In all, a total carbon supply of almost 33 Tg C yr -1 is delivered to the shelf by river loads together with coastline 25 erosion and atmospheric depositions.There is further a 1.7 Tg C yr -1 net uptake of atmospheric CO2 according to our model results.These carbon sources are balanced largely by a net export to the Arctic Ocean (31 Tg C yr -1 ), and to a lesser degree by burial (2.7 Tg C yr -1 ).
 The main sink terms for OCter are benthic mineralization and a net transport off the shelf.There are additional contributions from pelagic mineralization as well as burial in the sediments.DIC produced by OCter mineralization 30 contributes to both CO2 outgassing (or reduced uptake) and carbon export to the Arctic Ocean.
 Our results imply that without OCter mineralization, the net uptake of atmospheric CO2 would have been approximately 11.9 Tg C yr -1 larger.Thus, although the ESAS is an overall net sink for atmospheric CO2 according to our calculations, this sink could have been considerably more efficient.
 We calculate a summertime DIC excess of 10.6 Tg C, and this excess corresponds to the simulated decline in the 35 uptake of atmospheric CO2 that is a consequence of OCter mineralization on the shelf.quite heavily temperature limited, and cyanobacteria blooms are for that reason not expected to occur in the model simulations under present-day water temperatures.One functional group for heterotrophs represents all organisms that feed on phytoplankton and detritus.
Dead organic matter is included in the forms of dissolved and particulate organic carbon, nitrogen and phosphorus as well as biogenic silica (diatom shells).Detrital carbon is separated into two state variables, distinguishing POCter from POCmar 30 (Gustafsson et al., 2014a).The same distinction is made for DOC; DOCter and DOCmar pools are further separated into refractory and labile fractions.In all, four state variables are thus used to account for DOC in the model.DOCter is imported by rivers and atmospheric deposition, and produced through solubilization of POCter, whereas DOCmar is produced by phytoplankton exudations, zooplankton excretions, and solubilization of POCmar.
The depth dependent sediment pools of organic nitrogen, phosphorus, biogenic silica, as well as OCter and OCmar are 35 replenished by sedimentation of particulate organic matter.Resuspension, near-bottom transports and resettlement of organic matter is parameterized as a successive downward transport.Temperature dependent mineralization of organic matter occurs both in the water column and in the sediments.A fraction of the settled organic matter is permanently buried, while the majority is mineralized and returned to the water column as inorganic carbon and nutrients.Discuss., doi:10.5194/bg-2017Discuss., doi:10.5194/bg- -115, 2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.

Biogeosciences
An oxygen dependent fraction of the phosphate produced by benthic mineralization is released into the overlying water while the remaining fraction is retained in the sediments.Further, a fraction of the nitrate produced by mineralization and nitrification in the sediments is released to the water column whereas an oxygen dependent remainder is denitrified.If the oxygen concentration were to approach zero, first nitrate and then sulfate would be used instead of oxygen to oxidize the organic matter (denitrification and sulfate reduction).Under anaerobic conditions all mineralized nitrogen is released to the 5 water column as ammonium.Further, sulfate reduction would result in a production of hydrogen sulfide in anoxic waters (although we would not actually expect anoxic conditions in the ESAS).
Mineralization of organic carbon produces CO2 which equilibrates with bicarbonate and carbonate, and further adjusts fCO2 and the air-sea CO2 fluxes.Most biogeochemical processes result in either a production or consumption of AT; internal AT sources and sinks are calculated according to Wolf-Gladrow et al. (2007).It has been shown for the Baltic Sea and other 10 coastal-and shelf seas that internal AT generation in anoxic sediments can be a considerable source of AT to the water column (cf.Gustafsson et al., 2014b and references therein).Internal AT generation can increase the CO2 buffering capacity and in extension enhance the adsorption of atmospheric CO2 (or moderate CO2 evasions).It is however not known whether or not such processes are significant in the ESAS.

Appendix C. Stable carbon isotopes 15
Stable carbon isotopes are taken into account by adding one extra state variable for each carbon containing state variable (DIC, all plankton groups, terrestrial and marine detrital and dissolved organic carbon, terrestrial and marine benthic organic carbon; cf.Table A1-A2, supplementary material).Thus, 12 C and 13 C concentrations are modelled separately, and the isotope signatures for e.g.CO2, DIC, DOC, and POC can for that reason be calculated.The carbon fractionation related to air-sea CO2 fluxes, biogeochemical processes (production and respiration), and external loads were described by Gustafsson et al. (2015), and will 20 here only be briefly recapped.
The isotopic composition (in terms of 12 C and 13 C) of a carbon sample can be described by its ratio (R = 13 C/ 12 C) and/or δ values.The δ value represents a comparison between the R value of a sample and that of a standard (e.g.Zeebe and Wolf-Gladrow, 2001): 25 δ 13 Csample = (Rsample/Rstandard -1)*1000 (‰) A negative δ value means that the sample is depleted in 13 C or enriched in 12 C depending on the view taken compared to the standard (and vice versa).
A preferential assimilation of 12 C (i.e., isotope fractionation) by phytoplankton has the effect that the more CO2 that 30 is assimilated, the more the δ value of DIC in the productive layer increases.Similarly, 12 C travels more readily between air and water than 13 C. Uptake of atmospheric CO2 will for that reason decrease the δ value of surface water DIC, whereas CO2 outgassing on the other hand serves to increase the δ value.Mineralization of OC will affect the isotopic composition of the DIC pool slightly differently depending on OC origin as OCter usually is somewhat more depleted in 13 C than OCmar.
Mineralization of OCter consequently produces a more strongly 13 C depleted signal than mineralization of OCmar.The effect of 35 mineralization on δ 13 C values can be partially masked by compensating processes such as a preferential uptake of 12 C by phytoplankton or a preferential release of 12 C to the atmosphere.
The simulated δ 13 C value of OCmar depends on the δ value of assimilated CO2 (and temperature), which has the advantage that it is not necessary to assume a fixed isotope value for OCmar in the model.As the model conserves mass of both 12 C and 13 C, a Rayleigh distillation effect (e.g.Zeebe and Wolf-Gladrow, 2001) occurs in the productive layeri.e., when 40 primary producers preferentially assimilate 12 C, the DIC pool is gradually becoming enriched in 13  Tables Table 1: Properties of the three sub-basins (cf.Fig. 1); surface area (km 2 ), volume (km 3 ), average depth (m), and maximum depth (m).

Notation
At the western boundary, there is an exchange of North Atlantic and ESAS waters between the Kara Sea and LS.At the Eastern boundary, there is an exchange between North Pacific and ESAS waters between the Chukchi Sea and ESS.The entire study area has an open boundary toward the Arctic Ocean in the north.While the North Atlantic water is characterized by high salinity and temperature, the North Pacific water is in comparison low-saline and coldresulting in strong horizontal 20 gradients.
3 and Sim. 4 experiments, the net CO2 uptake would be approximately 10 Tg C yr -1 higher without the mineralization of eroded POCter.Mineralization of riverine DOCter and POCter on the other hand, together results in a 1.8 Tg C yr -1 decline in CO2 absorption.Based on 10 observations in the summer of 2008, Anderson et al. (2009) found an excess DIC of approximately 10 Tg C in the ESAS.The DIC excess was computed as the difference between observed DIC profiles and reference DIC profiles (where reference DIC profiles were calculated from observed AT profiles and the atmospheric CO2 partial pressure).Here we shall for comparison calculate excess DIC profiles based on simulated July and August DIC and AT profiles (in the year 2014), and an atmospheric pCO2 value of 390 µatm.Integrating the excess DIC profiles throughout the water column of the sub-basins results in an 15 Biogeosciences Discuss., doi:10.5194/bg-2017-115,2017 Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.a large removal of DOCter in the model has the effect that it becomes impossible to reproduce observed DOC and δ 13 CDOC values.When it comes to pathways for OCter supplied to the shelf, there are several parameters that are not known in detail: i) labile versus refractory fractions of riverine DOCter and POCter respectively, ii) the loads of POCter and DOCter related to coastline erosion (both total loads and distribution between the sub-basins), and iii) the labile versus refractory fractions of 5 eroded POCter/DOCter respectively.
2000-2014 OCter fluxes, net air-sea CO2 flux, and change in the flux induced by OCter degradation (Tg C yr -1 ), pools (Tg C), and estimated residence time (yr) (cf.Section 3.3) for the ILS, OLS, and ESS sub-basins as well as the total ESAS according to the Sim. 5 experiment (cf.

Figure 2 :
Figure 2: Observed (circles) salinity and temperature profiles compared to the corresponding simulated mean ± 1 std values (model data cover the period July-September 2000-2014).

Figure 5 . 5 Biogeosciences
Figure 5. Schematic overview of simulated carbon pools (Tg C) and fluxes (Tg C yr -1 ) in the ESAS (2000-2014 mean values).Net carbon sources include river loads, coastline erosion, atmospheric deposition, and air-sea CO2 flux.Net sinks include permanent burial and export to the Arctic Ocean.The figure further includes phytoplankton POC production and extracellular DOC release, POC sedimentation, OC mineralization in the sediments, and the accompanying DIC efflux to the water column.5

Figure B1 :
Figure B1: Schematic illustration of biogeochemical state variables and processes in the model.

Table 2
) because of DOCter mineralization that resulted in elevated fCO2 in the water column.The decreased CO2 uptake in Sim. 2 was compensated by a corresponding decrease in carbon export to BiogeosciencesDiscuss., doi:10.5194/bg-2017Discuss., doi:10.5194/bg--115,2017Manuscript under review for journal Biogeosciences Discussion started: 3 April 2017 c Author(s) 2017.CC-BY 3.0 License.