Introduction

Pan-Arctic river discharge to the Arctic Ocean has increased over recent decades (Peterson et al. 2002; Haine et al. 2014) with climate model simulations forecasting this to continue and potentially accelerate during the 21st century (Ahmed et al. 2020; Wang et al. 2021). Hydrologic models informed with climate projections estimate increases of ~ 25 to 50% in freshwater discharge to the Laptev and East Siberian Shelf by 2100 (Wang et al. 2021). Ongoing terrestrial permafrost thaw and enhanced rates of coastal erosion in response to climate warming will additionally alter the source and loads of dissolved organic matter (DOM) exported from watersheds to coastal waters (Wang et al. 2021; Frey and McClelland 2008). An intensifying hydrologic cycle will therefore cause greater quantities of terrestrial dissolved organic carbon (terr-DOM) to reach coastal shelf waters, with changes in its overall composition due to redistribution of aged peat and permafrost derived materials (Mann et al. 2022; McGuire et al. 2009; Vonk and Gustafsson 2013).

The shallow Laptev Sea (50 m) receives the greatest quantities of freshwater across the ESAS (∼ 745 km3 year−1), predominantly supplied by the Lena River (566 km3 year−1; Cooper et al. 2008). Shelf waters of the Laptev Sea are characterized by high DOM concentration (up to 500 µM, Juhls et al. 2019) comprised mostly of terr-DOM, evidenced by stable isotope carbon composition and ratios consistent with terrestrial sources (Alling et al. 2010; Salvadó et al. 2016). Arctic shelf DOM also comprises DOM derived from coastal erosion and marine production, and is influenced by the contrasting fate (export, mineralization, flocculation, and burial) of different DOM fractions (e.g., non-humic vs. humic DOC) across the shelf. Overall losses of up to 10 % of the riverine DOC pool have been estimated across the Lena River estuary before export onto the ESAS shelf, likely due to processes such as photodegradation, flocculation and sedimentation (Alling et al. 2010; Gustafsson et al. 2000) which mainly remove humic DOM fractions in the lower salinity zone (between 4 and 6 psu, Forsgren et al. 1996; Eckert and Sholkovitz 1976; Stubbins et al. 2016). By contrast, moving offshore, where salinity is higher, DOC losses in the non-humic fraction of the DOM pool become increasingly more important, indicating a dominant role for bacteria in degrading terr-DOM in the outer shelf regions (Anderson et al. 2019; Alling et al. 2010; Amon and Benner 2003; Lobbes et al. 2000).

The fate and impact of terr-DOM on shelf biogeochemistry is, therefore, influenced by the capacity of heterotrophic bacteria to utilize terr-DOM to fulfil their carbon and nutrient requirements i.e., on the biological lability of terr-DOM, a key yet highly uncertain parameter (Alling et al. 2010; Holmes et al. 2008; Manizza et al. 2009). Despite this, the sensitivity of air–sea CO2 flux estimates to contemporary or future terr-DOM biological degradability has yet to been addressed. However, mounting evidence suggests future changes to coastal terr-DOM composition, due to permafrost thaw and reduced freshwater residence times (Mann et al. 2022), will increase the biological lability of terr-DOM promoting enhanced DOC losses (Catalan et al. 2016; Mann et al. 2015; Vonk et al. 2013) and supply of nutrients to support coastal productivity across Arctic shelves (Terhaar et al. 2019, 2021). Here, we use a complex marine biogeochemical model—the European Regional Seas Ecosystem Model (ERSEM, Butenschön et al. 2016) specifically augmented with a new formulation accounting for terr-DOM to assess the impact of terr-DOM concentration and its biological reactivity on the CO2 exchange between the ocean and the atmosphere.

Materials and methods

The biogeochemical model

ERSEM is a biomass and functional type based marine biogeochemical model describing carbon and nutrient cycling within the planktonic and benthic environment. Heterotrophic prokaryotes (DOM consumers) are modelled through a single functional type, hereafter generically called “bacteria”. DOM–bacteria interactions, including terr-DOM are summarized in Fig. 1.

Fig. 1
figure 1

Model schematic describing the interactions between bacteria and organic carbon. On the right side is described the dynamics of autochthonous DOC as modelled in ERSEM (Butenschön et al. 2016). On the left side is described the addition of terr-OC modelled based on Anderson et al. (2019)

ERSEM is fully described in Butenschön et al. (2016), which includes a detailed description of model assumptions and mathematical equations. The current implementation is based on the Butenschön et al. (2016) description with 4 phytoplankton and 3 zooplankton functional groups. Here we limit our description to the equations accounting for terr-DOM which have been newly developed in this work partially based on Anderson et al. (2019) and to bacterial respiration which is a key variable for the purpose of this study. terr-DOM is described through two state variables: T1 and T2 (both given in mg C m−3). T1 represents a coloured, light-absorbing ‘humic’ fraction which is susceptible to both microbial and photochemical degradation. T2 by contrast, accounts for the non-light absorbing ‘non-humic’ fraction and is susceptible only to degradation by bacteria. Following Anderson et al. (2019), the T1 fraction is also assumed to be prone to flocculation processes. terr-DOM is assumed to have a fixed C:N and C:P molar ratio equal to 22 and 950, respectively (Cauwet and Sidorov 1996; Kutscher et al. 2017; Sanders et al. 2022). Model parameters, other than those explicitly described here, are taken from Butenschön et al. (2016).

The evolution of terr-DOM can be summarized as:

$$\frac{{\text{dT}}_{1,2}}{\text{dt}}=\text{uptake}-\text{photodeg}-\text{flocc}-\text{diffusion}-\frac{1}{\text{R}}\left({\text{T}}_{1,2}-{\text{T}}_{\text{obs}}\right),$$
(1)

where \({\text{T}}_{1,2}\) is the terr-DOM fraction concentration, \(\text{uptake}\) is the bacterial consumption (see below) and \(\text{flocc}, \text{photodeg}\; \text{and}\; \text{diffusion}\) are the terr-DOM changes due to flocculation, photodegradation and turbulent diffusion, respectively. R is a relaxation constant (10 days). The purpose of R is to enable the introduction of a source of terr-DOM (\({\text{T}}_{\text{obs}},\) driven by the seasonality of the Lena outflow) in our 1D model setup where horizontal advection is not resolved.

The term uptake (mg C m−3 day−1) is equal to:

$$\text{uptake}= \frac{{{\alpha }}_{1,2} \cdot \text{T}_{1,2}}{\text{Total substrate}} \cdot \text{min}\left(\text{Pot},\text{Total substrate}\right)$$
(2)

where \(\text{Pot}\) (mg C m−3 day−1) is the potential bacterial uptake (i.e. the uptake under non-limiting carbon conditions) and is given by:

$$\text{Pot}=\text{r} \cdot \text{B} \cdot {\text{T}}^{10} \cdot \text{O}2\text{lim}$$
(3)

where \(\text{r}\) (day−1) is the max specific uptake rate of bacteria, \(\text{B}\) is the bacterial concentration (mg C m−3) and \({\text{T}}^{10}\) and \(\text{O}2\text{lim}\) are adimensional functions describing temperature and O2 dependency, respectively (Butenschön et al. 2016). \({{\alpha }}_{1,2}\) are non-dimensional parameters describing the degree of lability of T1 and T2 with respect to the degradability of the labile DOC fraction (\({\sigma},\) 1 day−1, Butenschön et al. 2016; Polimene et al. 2006).

\(\text{Total substrate}\) (mg C m−3 day−1) is the sum of all the carbon sources available to bacteria multiplied by the relative lability coefficient \({{\alpha }}_{\text{n}}\) and \({\sigma }\):

$$\text{Total substrate} = T_{1} \cdot \alpha _{1} \cdot \sigma + T_{2} \cdot \alpha _{2} \cdot \sigma + \sum\limits_{1}^{n} {DOM_{n} } \cdot \alpha _{n} \cdot \sigma + \sum\limits_{1}^{m} {POM_{m} } \cdot \alpha _{m} \cdot \sigma ,$$
(4)

where DOM and POM are the dissolved and particulate organic matter produced by the marine planktonic system respectively, and \({{\alpha }}_{\text{n}}\) and \({{\alpha }}_{\text{m}}\) the relative lability coefficients (Butenschön et al. 2016). The number of DOM and POM fractions (n and m) described in ERSEM map onto the different functional groups considered (Butenschön et al. 2016).

Note that if \(\text{Total substrate} < {\text{Pot}}_{{upt}}, {{\alpha }}_{1,2}\) is a first order degradation rate of T1,2.

The term \(\text{photodeg}\) (mg C m−3 day−1) is only applied to T1 and is given by:

$$\text{photodeg}= {\phi } \cdot \frac{\text{I}}{{\text{I}}_{\text{ref}}} \cdot \text{T}_1$$
(5)

where \({\phi }\) is the reference photodegradation rate (assumed to be 0.03 day−1, Mann et al. 2012), \(\text{I}\) the incident irradiance (W m−2) and \({\text{I}}_{\text{ref}}\) the reference irradiance assumed to be 130 W m−2. According to Anderson et al. (2019), byproducts of photodegradation are partially directed into the inorganic carbon and nutrient pools, while the remaining fraction (0.2) is redirected into the T2 component. It should be stressed that we assumed photodegradation to be a function of the total irradiance (Ultraviolet or UV light is not modelled) and this could lead to an overestimation of this process as UV absorption by water is significantly higher than at larger wavelengths. However, since photodegradation parameters were not altered during our experiments, this does not affect the quality of our results. Loss due to flocculation (which only applies to T1) is described by the term \(\text{flocc}\) (mg C m−3 day−1)

$$\text {flocc} = \vartheta \cdot e^{{ - \frac{{({\text{ln }}S - S_{x} )}}{{2 \cdot \beta ^{2} }}}} \cdot T1^{2}$$
(6)

where \({S}\) is salinity in psu, \({\vartheta }\) is the maximum flocculation rate [2 × 10−6 day−1 (mg C)−1, Anderson et al. 2019], \({\text{S}}_{\text{x}}\) is the ln of salinity at which flocculation is maximum (~ 2 psu) and \({\beta}\) is the parameter of the bell-shaped curve. We have used \({\beta}=1.35,\) implying that the max flocculation rate decreases by one order of magnitude at 35 psu.

Bacterial respiration (\(\text{RESP},\) mg C m−3 day−1) is given by the sum of two terms describing activity and rest respiration, respectively:

$$\text{RESP}= \text{uptake} \cdot \left[1-\text{Eff}-\text{EffO}2 \cdot \left(1-\text{O}2\text{lim}\right)\right]+{\text{R}}_{\text{r}} \cdot \text{B} \cdot {\text{T}}^{10}$$
(7)

where \(\text{Eff}\) (unitless) is the maximal bacterial growth efficiency (0.6, Butenschön et al. 2016), \(\text{EffO}2\) the bacterial growth efficiency under oxygen limitation (0.2, Butenschön et al. 2016) and \({\text{R}}_{\text{r}}\) the daily carbon specific rest respiration rate (0.1 day−1, Butenschön et al. 2016).

Model physical set up

We coupled the model described above with a one-dimensional hydrodynamic model, the General Ocean Turbulence Model (Burchard et al. 1999) and implemented it at the shelf-ocean break of the Siberian Laptev Sea (Lat 75, Lon 130, depth ~ 38 m). This site is well documented to receive significant quantities of terr-DOM from the Lena River and coastal erosion (Bauch et al. 2013; Juhls et al. 2019).

The model was forced with tidal data (Egbert and Svetlana 2002) and reanalysis meteorological data, including net short wave radiation (hourly ERA-5, generated using Copernicus Climate Change Service information, Hersbach et al. 2018) and initialized with vertically uniform values of nutrients (mmol m−3 of NO3, PO4 and SiO2) estimated from literature (NO3 = 4, PO4 = 0.8 and SiO2 = 8, Kattner et al. 1999; Sorokin and Sorokin 1996; Bauch and Cherniavskaia 2018).

To reproduce the seasonal water column structure evolution, temperature and salinity profiles were constrained between artificially created summer (August) and winter (January) profiles (Figs. S2 and S3). Summer temperature and salinity profiles range from 4 degrees and 10 psu, respectively at the surface to − 1.5° and 34 psu at the bottom (Bauch et al. 2013) while winter values were constant through the water column and equal to the summer bottom values (i.e. − 1.5 and 34 psu, respectively). These profiles were linearly interpolated to generate seasonal time series which were assimilated by the model. Resulting temperature and salinity values, averaged over the light season, were consistent with the data reported in Semiletov et al. (2016) (SI, Figs. S1 and S2).

Light extinction through the water column is dependent on the concentrations of organic particulates (Butenschön et al. 2016) and the coloured fraction of terr-DOM (T1). Light absorption by inorganic particles (which are not explicitly modelled) has been simulated by using the mass specific light absorption for inorganic suspended solid matter reported in Butenschön et al. (2016) and a concentration of 800 mg m−3. The latter was estimated from the total suspended matter measured in the proximity of the model location (~ 1 g m−3, Wegner et al. 2003) assuming a 20% contribution of living and detrital organic particles. T1 has been added to the light absorbing substrates assuming a specific absorption coefficient of 0.0002 (m2 mg C−1). This value is at the lower range of those reported in Matsuoka et al. (2014) for the ESAS. At the lower boundary of the water column, a simple remineralization closure is applied (Polimene et al. 2014). According to this, sinking organic particles are exported from the water column and re-injected into the bottom waters as dissolved nutrients and inorganic carbon.

Ice coverage reduction due to global warming is expected to affect primary production (PP) in the Arctic Seas (Lewis et al. 2020) with potential consequence on CO2 fluxes. To account for this, we used a conservative approach by not considering ice coverage in the model setup. This implies model phytoplankton experience the longest productive season theoretical possible i.e. coinciding with the light season. The light season was here defined based on modelled non-zero values of daily photosynthetically active radiation simulated from March to October. Simulations were run for 9 years (2010–2018) to check for model stability (i.e. absence of drift in state variables). After 2 years of simulation, a stable repeating cycle for the main biogeochemical variables was achieved and the year 2012 was taken for the sensitivity exercise. The physical setup of the model has been kept constant in all the scenarios tested.

Terr-DOM scenarios

The starting terr-DOM concentration (present day scenario) was estimated from late summer DOC observations (on average ~ 200 µmol L−1, Juhls et al. 2019) and assuming a 60% terr-DOM content (Kattner et al. 1999). Spring and Winter terr-DOM concentrations were then estimated from the seasonal evolution of the Lena River DOC discharge (Cauwet and Sidorov 1996; Juhls et al. 2020). terr-DOM loads during May to June and over winter months (January–March) are approximately double and half (respectively) that of late summer concentrations (August). These values were linearly interpolated to construct a terr-DOM seasonal time series (\({\text{T}}_{\text{obs}}\)) which was assimilated by the model into the T1 and T2 state variables (50% each) as described in Eq. 1. Present day terr-DOM fluxes (SI, Fig. S3) were increased by 25, 50, 75 and 100% covering the range of increases estimated for Arctic Shelf Seas (20–46%, Frey and Smith 2005; Frey et al. 2007) and model scenarios (up to 100%, Terhaar et al. 2019).

Terr-DOM biological lability is expressed as life time i.e. the time by which a terr-DOM pool [X] is degraded to a value equal to [X]/e (Hansell 2013). A range of life time parameters were used based upon decomposition rates reported in literature: 0.7 years (river and estuarine water incubations; Holmes et al. 2008; Vonk et al. 2013), 3.3 years (field observations; Alling et al. 2010) and 10 years (modelling estimate; Manizza et al. 2009). To examine more extreme future changes in terr-DOM lability (e.g. enhanced permafrost thaw; Mann et al. 2022), we widened the life times examined to include 0.3 and 20 years. The inverse of these life times (degradation rates) have been used as input values for the model parameters \({{\alpha }}_{1}\) and \({{\alpha }}_{2}\). Total terr-DOM is assumed to be composed by the same amount of T1 and T2 (i.e. 50–50%) while, according to Anderson et al. (2019), T2 is considered to be ~  3 times more labile than T1. The values of \({{\alpha }}_{1}\) and \({{\alpha }}_{2}\) used are: \({{\alpha }}_{1}=[0.012, 0.006, 0.0012, 0.0003, 0.00015].\) and \({{\alpha }}_{2}=[0.012, 0.006, 0.0012, 0.0003, 0.00015].\) The average of each pair of these values corresponds to the life times described above (from 0.3 to 20 years).

The model was run under four configurations each combining the terr-DOM concentrations and labilities described above (i.e. up to 25 runs, SI Table S1):

  1. (1)

    Core experiment with parameters and initial conditions as above.

  2. (2)

    Experiment S1 assuming terr-DOM was totally refractory to bacteria (i.e. \({{\alpha }}_{1}\) and \({{\alpha }}_{2}\) set to zero)

  3. (3)

    Experiment S2 with doubled nutrient initial conditions (nitrate, phosphate and silica).

  4. (4)

    Experiment S3 with double N and P content in terr-DOM (i.e. with double N:C and P:C ratios).

Model performance in the core experiment was assessed by comparing simulated PP with satellite and field estimates presented in literature (SI, Fig. S4).

Model caveats and limitations

The presented modelling approach has caveats and limitations mainly due to the simulation of the physical features of the system. The effect of climate changes on the physical structure of the Arctic shelf is likely to be complex (Timmermans and Marshall 2020; Årthun et al. 2021) and our theoretical model setup was not meant to reproduce such complexity. Consequently, the simulation of the biogeochemical and ecological dynamics we have highlighted is quantitatively affected by the uncertainty associated with our simplified physical setup. However, we stress that our goal was to investigate potential trends in CO2 fluxes rather than providing quantitative estimates.

Results and discussion

Air-sea CO2 fluxes as a function of terr-DOM concentration and life time were examined using our model (Fig. 2). In an idealized water column as here, we could expect a CO2 balance close to zero in the absence of allochthonous DOM due to the close interdependence and balance between CO2 consumption (driven by primary productivity) and production driven by community respiration. Positive, seasonally averaged air to sea CO2 fluxes (i.e. CO2 uptake) are possible representing produced carbon that is buried in sediments and/or respired over longer time scales. Addition of terr-DOM in such a system should alter this balance potentially leading to dissolved inorganic carbon accumulation in the upper part of the water column.

Assuming initial conservative life time estimates of 10–20 years and present day terr-DOM supply, our model estimates a small annual CO2 uptake (~ 4 mmol CO2 m−2 day−1) in good agreement with recent estimates for the Arctic basin (4 ± 4 mmol CO2 m−2 day−1, Yasunaka et al. 2016). Air-sea CO2 fluxes were sensitive to both increased terr-DOM supply and/ or reductions in average life times (Fig. 2). Increased terr-DOM supply alone only shifted shelf waters toward a net CO2 source under the doubled scenario (+ 100% terr-DOM) when life times were higher (3.3, 10 and 20 years). However, air–sea fluxes were highly responsive to terr-DOM life times, with shelf waters shifting to net CO2 sources under all terr-DOM supply scenarios- including under present day terr-DOM loads, when shorter life times of 0.7 and 0.3 years were examined. Concurrent changes in both terr-DOM load and life times, caused more rapid shifts toward net CO2 emissions.

Fig. 2
figure 2

CO2 air-to-sea flux as function of terr-DOM concentration and life time. Positive values indicate CO2 sinks, and negative values net sources. Values have been averaged over the light season

Increased CO2 air sea emissions were proportional to increased bacterial terr-DOM uptake, which in turn varied across the terr-DOM scenarios tested (Fig. 3a). Our model indicates that < 10% of total DOC uptake is subsidized by terr-DOM when life times > 3.3 years, whereas terr-DOM contributed up to 60% of the DOC assimilated when life times were < 1 year (Fig. 3a). Increased future terr-DOM subsidies to marine bacteria may therefore be expected to drive greater quantities of CO2 production through respiration (del Giorgio and Cole 1998).

Fig. 3
figure 3

a Percentage of bacterial DOC uptake subsidized by terr-DOM as function of terr-DOM concentration (x axis) and lifetime (colours) and b PP as a function of terr-DOM concentration (x-axis) and lifetime (colours). Fluxes have been depth-integrated and averaged over the light season

Future increases in bacterial respiration may be partially balanced by concomitant increases in PP, due to the fertilization effect of nutrients associated with terr-DOM. However, although incorporating the impacts of nutrient enrichment from terr-DOM (Fig. S5), our model simulated decreased PP under all scenarios of increasing terr-DOM concentrations (Fig. 3b) causing a positive feedback to CO2 outgassing (Fig. 2).

Patterns of decreasing PP were due to the interplay of two factors: light limitation and increased grazing pressure on primary producers. When terr-DOM was assumed to have longer life times (> 10 years), declining PP was primarily caused by reduced light penetration caused by greater terr-DOM light absorption in the water column. Model simulations repeated assuming terr-DOM was completely refractory (acting as a “light screen” alone) demonstrated PP decreases between 4 and 16% with 25 to 100% increases in terr-DOM supply (experiment S1; Fig. 4). By contrast, when terr-DOM lifetimes were shorter, increased grazing pressure on primary producers became progressively more important in all terr-DOM scenarios (Figs. 5 and S6). When terr-DOM was more biologically available, bacterial biomass became higher during winter and early spring months (i.e. before the phytoplankton bloom, Fig. S7), inducing a shift to a more heterotrophic dominated food chain with relatively high zooplankton concentrations (Fig. 5 and Fig. S6). Increased zooplankton control phytoplankton through grazing, reducing PP. As such, the combined effects of light reduction and increased zooplankton pressure resulted in more pronounced declines in PP (up to 60% PP reduction in the + 100% terr-DOM scenario), further decoupling CO2 production from consumption and causing sharp increases in projected CO2 outgassing (Fig. 2). Notably, for each given terr-DOM concentration, decreasing terr-DOM lifetime strongly enhanced the CO2 outgassing (up to 20 times in the + 100% scenario with 0.3 years terr-DOM lifetime).

Fig. 4
figure 4

Decrease in average, depth-integrated primary production (a) and photosynthetically active radiation simulated at 5-m depth (b) with respect to the present day scenario in the experiment S1. Uncertainties have been estimated by considering the Coefficients of Variation (CV = (std/mean) × 100) in the + 100% scenario normalized by the CVs in the present day

Fig. 5
figure 5

Relative increase in zooplankton biomass as function of terr-DOM life times in the + 100% scenario with respect to present day. Uncertainties have been estimated by considering the coefficients of variation (CV = (std/mean) × 100) in the + 100% scenario normalized by the CVs in the present day

Recent modeling studies suggest nitrogen subsidies from terr-DOM may support up to 51% of contemporary Arctic shelf PP and infer that terrestrial nutrient supplies from land will affect the future evolution of the Arctic Ocean PP (Terhaar et al. 2021). To examine the relative importance of land-derived nutrient supply, we compare PP simulated under present day scenarios assuming high and low terr-DOM life times (Fig. 3). Model simulation incorporating high terr-DOM life time (20 years) examine PP when land-derived nutrients are limited in availability for use by primary producers, whereas low terr-DOM life time (0.3 years) assess PP when nutrients are readily accessible to phytoplankton. In good agreement with Terhaar et al. (2021), simulated PP under short terr-DOM life times (0.3 year−1) was approximately double of that simulated under 20 years terr-DOM lifetimes (Fig. 2b), suggesting terrestrially-derived nutrients remineralised by bacteria and their grazers have potential to sustain up to ~ 50% of simulated PP under the present day terr-DOM concentrations. However, our model results further suggest that with increasing terr-DOM concentrations, this “fertilization” effect is progressively offset by light limitation imposed by terr-DOM absorption and/or a stimulated zooplankton grazing on phytoplankton populations.

We re-ran our model scenarios (as Fig. 1) assuming: (1) a doubling in initial nutrient concentration inputs (experiment S2) and, (2) doubled nutrient content of the terr-DOM pool (e.g. doubled N:C and P:C ratios; experiment S3). Experiment S2 provides insights into nutrient supplies not associated to terr-DOM as, for example, through the additional supply of inorganic nutrients (PO4, NO3 and SiO2) advected from the Atlantic and Pacific Oceans (Lewis et al. 2020). Experiment S3, examines the effect of increased nutrients (organic P and N) within terr-DOM (i.e. assuming more organic nutrients per unit of terr-DOM). Increased nutrient loads only marginally impacted the CO2 air–sea fluxes in both experiments (Fig. 6). However, the response to nutrient enrichments differed between experiments. Ocean CO2 uptake increased under all scenarios tested in exp. S2 (Fig. 6a), whereas increased CO2 uptake was only simulated at low terr-DOM concentrations and high terr-DOM lifetime in exp. S3 (Fig. 6b).

Fig. 6
figure 6

Differences in CO2 average fluxes between the core experiment (Fig. 1) and a Experiment S2 and b Experiment S3. Negative differences (blue colour) indicate increased CO2 ocean uptake, positive values (red colour) increased CO2 emissions

In both nutrient enrichment experiments, increases in PP (Fig. S7) were partially balanced by increased heterotrophic respiration (bacterial respiration alone increases up to 18%), resulting in limited overall effects on CO2 uptake (Figs. 6 and S8). Large differences between the two experiments were simulated when terr-DOM concentrations and degradability were high; this is due to the dominance of diatoms under these conditions i.e. when smaller phytoplankton groups are top-down controlled by grazers which are enhanced by increased bacteria (Fig. S6). Consequently, in exp S3 which provided no SiO2 enrichment, increasing N and P associated to terr-DOM are less effective in fertilizing phytoplankton since PP is mainly limited by silica. Together, these results indicate that future increases in nutrients supply or even PP will not necessarily translate to proportional enhancements in net CO2 Ocean uptake.

Our results question the capacity of the Arctic Ocean to serve as a net sink for atmospheric CO2 in agreement with prior studies (Semiletov et al. 2016) and may help to explain recent satellite observations indicating decreasing PP in the Laptev Sea over recent decades (Demidov et al. 2020). Contrary to prior studies (e.g. Bates and Mathis 2009), we suggest that Arctic shelf waters are susceptible to shifting from net sinks to net atmospheric sources of CO2 under future changes in terr-DOM supply and origin. Additional nutrients from land-derived sources, or elsewhere, only marginally offset these trend due to complex concomitant changes in biogeochemical and ecological processes. Our results highlights the need for future studies to incorporate complex biogeochemical models explicitly accounting for bacterial dynamics to simulate emergent properties affecting CO2 production and consumption in coastal areas.

Terrestrial permafrost thaw and enhanced coastal erosion rates across the Arctic can mobilize peat and permafrost-derived terr-DOM to coastal waters, modifying nutrient and carbon loads and their relative biological degradability (Mann et al. 2015; Vonk et al. 2013). Small subsidies of permafrost-derived terr-DOM to riverine and estuarine waters (∼ 1%) have been shown to result in marked increases in biodegradability rates (20–60% dependent on DOM pool), and thus reduced mean life times (Mann et al. 2022). Changes to terr-DOM sources combined with an intensification in Arctic hydrologic cycles (e.g. 25–50% freshwater increase to ESAS shelf by 2100, Wang et al. 2021), is therefore likely to deliver greater quantities of more bioavailable terr-DOM to Arctic shelves over coming decades. Declining coastal sea-ice cover will also modify coastal light availability and timing. As our model assumptions implies that ice is absent during the light season, our model estimates already effectively incorporate future ice-free conditions and thus model simulated PP is likely overestimated with resulting CO2 fluxes tending to be conservative. To reliably quantify the future evolution of Arctic Ocean climate feedback complex biogeochemical models coupled to a three-dimensional hydrodynamic framework will be needed, on which to run climate-scenario simulations. Our study suggests that future inputs of terr-DOM from peat and permafrost thaw may induce net CO2 efflux from the Arctic shelf causing, currently unquantified, positive feedback to global warming.