Biological lability of terrestrial DOM increases CO2 outgassing across Arctic shelves

Arctic shelf seas receive greater quantities of river runoff than any other ocean region and are experiencing increased freshwater loads and associated terrestrial matter inputs since recent decades. Amplified terrestrial permafrost thaw and coastal erosion is exposing previously frozen organic matter, enhancing its mobilization and release to nearshore regions. Changing terrestrial dissolved organic matter (terr-DOM) loads and composition may alter shelf primary productivity and respiration, ultimately affecting net regional CO2 air–sea fluxes. However, the future evolution of Arctic Ocean climate feedbacks are highly dependent upon the biological degradability of terr-DOM in coastal waters, a factor often omitted in modelling studies. Here, we assess the sensitivity of CO2 air–sea fluxes from East Siberian Arctic Shelf (ESAS) waters to changing terr-DOM supply and degradability using a biogeochemical model explicitly accounting for bacteria dynamics and shifting terr-DOM composition. We find increasing terr-DOM loads and degradability trigger a series of biogeochemical and ecological processes shifting ESAS waters from a net sink to a net source of CO2, even after accounting for strengthening coastal productivity by additional land-derived nutrients. Our results suggest that future projected inputs of labile terr-DOM from peat and permafrost thaw may strongly increase the CO2 efflux from the Arctic shelf sea, causing currently unquantified positive feedback to climate change.

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 km 3 year −1 ), predominantly supplied by the Lena River (566 km 3 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 nonhumic 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 CO 2 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; and supply of nutrients to support coastal productivity across Arctic shelves (Terhaar et al. 2019(Terhaar et al. , 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 CO 2 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.
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: T 1 and T 2 (both given in mg C m −3 ). T 1 represents a coloured, lightabsorbing 'humic' fraction which is susceptible to both microbial and photochemical degradation. T 2 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 T 1 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: where T 1,2 is the terr-DOM fraction concentration, uptake is the bacterial consumption (see below) and flocc, photodeg and 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 ( T 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: where 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: where r (day −1 ) is the max specific uptake rate of bacteria, B is the bacterial concentration (mg C m −3 ) and T 10 and O2lim are adimensional functions describing temperature and O 2 dependency, respectively (Butenschön et al. 2016). 1,2 are non-dimensional parameters describing the degree of lability of T 1 and T 2 with respect to the degradability of the labile DOC fraction ( , 1 day −1 , Butenschön et al. 2016;Polimene et al. 2006). 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 n and : (2) uptake = 1,2 ⋅ T 1,2 Total substrate ⋅ min(Pot, Total substrate) (3) Pot = r ⋅ B ⋅ T 10 ⋅ O2lim Fig. 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) where DOM and POM are the dissolved and particulate organic matter produced by the marine planktonic system respectively, and n and 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 Total substrate < Pot upt , 1,2 is a first order degradation rate of T 1,2 .
The term photodeg (mg C m −3 day −1 ) is only applied to T1 and is given by: where is the reference photodegradation rate (assumed to be 0.03 day −1 , Mann et al. 2012), I the incident irradiance (W m −2 ) and I 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 T 1 ) is described by the term flocc (mg C m −3 day −1 ) where S is salinity in psu, is the maximum flocculation rate [2 × 10 −6 day −1 (mg C) −1 , Anderson et al. 2019], S x is the ln of salinity at which flocculation is maximum (~ 2 psu) and is the parameter of the bellshaped curve. We have used = 1.35, implying that the max flocculation rate decreases by one order of magnitude at 35 psu. (4) Bacterial respiration ( RESP, mg C m −3 day −1 ) is given by the sum of two terms describing activity and rest respiration, respectively: where Eff (unitless) is the maximal bacterial growth efficiency (

Model physical set up
We coupled the model described above with a onedimensional 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).
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. T 1 has been added to the light absorbing substrates assuming a specific absorption coefficient of 0.0002 (m 2 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 CO 2 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 ( T obs ) which was assimilated by the model into the T 1 and T 2 state variables (50% each) as described in Eq. 1. Present day terr-DOM fluxes (SI, Fig. S3 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) Core experiment with parameters and initial conditions as above.
(2) Experiment S1 assuming terr-DOM was totally refractory to bacteria (i.e. 1 and 2 set to zero) (3) Experiment S2 with doubled nutrient initial conditions (nitrate, phosphate and silica). (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 CO 2 fluxes rather than providing quantitative estimates.

Results and discussion
Air-sea CO 2 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 CO 2 balance close to zero in the absence of allochthonous DOM due to the close interdependence and balance between CO 2 consumption (driven by primary productivity) and production driven by community respiration. Positive, seasonally averaged air to sea CO 2 fluxes (i.e. CO 2 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 CO 2 uptake (~ 4 mmol CO 2 m −2 day −1 ) in good agreement with recent estimates for the Arctic basin (4 ± 4 mmol CO 2 m −2 day −1 , Yasunaka et al. 2016). Air-sea CO 2 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 CO 2 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 CO 2 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 CO 2 emissions.
Increased CO 2 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 CO 2 production through respiration (del Giorgio and Cole 1998).
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 CO 2 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 Fig. 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 Fig. 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 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 CO 2 production from consumption and causing sharp increases in projected CO 2 outgassing (Fig. 2). Notably, for each given terr-DOM concentration, decreasing terr-DOM lifetime strongly enhanced the CO 2 outgassing (up to 20 times in the + 100% scenario with 0.3 years terr-DOM lifetime).
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 landderived 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 (PO 4 , NO 3 and SiO 2 ) 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 CO 2 air-sea fluxes in both experiments (Fig. 6). However, the response to nutrient enrichments differed between experiments. Ocean CO 2 uptake increased under all scenarios tested in exp. S2 (Fig. 6a), whereas increased CO 2 uptake was only simulated at low terr-DOM concentrations and high terr-DOM lifetime in exp. S3 (Fig. 6b).
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 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 effects on CO 2 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 SiO 2 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 CO 2 Ocean uptake.
Our results question the capacity of the Arctic Ocean to serve as a net sink for atmospheric CO 2 in agreement with prior studies  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 CO 2 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 CO 2 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;. 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 seaice cover will also modify coastal light availability and timing. As our model assumptions implies that Fig. 6 Differences in CO 2 average fluxes between the core experiment (Fig. 1) and a Experiment S2 and b Experiment S3. Negative differences (blue colour) indicate increased CO 2 ocean uptake, positive values (red colour) increased CO 2 emissions 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 CO 2 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 CO 2 efflux from the Arctic shelf causing, currently unquantified, positive feedback to global warming.