Modelling the impact of riverine DON removal by marine bacterioplankton on primary production in the Arctic Ocean

The planktonic and biogeochemical dynamics of the Arctic shelves exhibit a strong variability in response to Arctic warming. In this study, we employ a biogeochemical model coupled to a pan-Arctic ocean–sea ice model (MITgcm) to elucidate the processes regulating the primary production (PP) of phytoplankton, bacterioplankton (BP), and their interactions. The model explicitly simulates and quantifies the contribution of usable dissolved organic nitrogen (DON) drained by the major circum-Arctic rivers to PP and BP in a scenario of melting sea ice (1998–2011). Model simulations suggest that, on average between 1998 and 2011, the removal of usable riverine dissolved organic nitrogen (RDON) by bacterioplankton is responsible for a ∼ 26 % increase in the annual BP for the whole Arctic Ocean. With respect to total PP, the model simulates an increase of ∼ 8 % on an annual basis and of ∼ 18 % in summer. Recycled ammonium is responsible for the PP increase. The recycling of RDON by bacterioplankton promotes higher BP and PP, but there is no significant temporal trend in the BP : PP ratio within the ice-free shelves over the 1998–2011 period. This suggests no significant evolution in the balance between autotrophy and heterotrophy in the last decade, with a constant annual flux of RDON into the coastal ocean, although changes in RDON supply and further reduction in sea-ice cover could potentially alter this delicate balance.


Introduction
In response to the polar amplification of global climate change, air temperature in the lower atmosphere is increasing twice as fast in the Arctic as in temperate regions.By the end of the century, model projections suggest an average increase in the surface air temperature of 3.7 • C relative to 1981-2000(ACIA report, 2005)).In response to Arctic warming, plankton production and the biogeochemistry of the Arctic Ocean (AO) are rapidly evolving.Changes in phytoplankton communities (Li et al., 2009) as well as their phenology in spring (Kahru et al., 2011) and autumn (Ardyna et al., 2014) are being observed.Overall, the AO tends to be more productive (Bélanger et al., 2013) and is taking up more atmospheric carbon dioxide (1996-2007;Manizza et al., 2013).In the long term, model projections suggest an increase in spatially integrated primary production (PP) by the end of the twentyfirst century (Vancoppenolle et al., 2013).
The AO is the basin most influenced by continental freshwater.It receives 10 % of the freshwater that flows into the global ocean, but represents only 1 % of the global ocean volume (Opshal et al., 1999).Circum-Arctic rivers are potentially a significant source of inorganic nutrients and organic matter for shelf seas (Le Fouest et al., 2013;Tank et al., 2012).10 % of the global riverine inputs of organic carbon are conveyed into the AO (Rachold et al., 2004).This fraction is projected to increase in the near future due to the ac-V.Le Fouest et al.: Modelling the impact of riverine DON removal by marine bacterioplankton celerated thawing of permafrost (Frey et al., 2007).This pool of organic matter enters the carbon cycle, but little is known about its fate and pathways within the plankton ecosystem in Arctic waters prior to being exported into the Atlantic Ocean.
Bacterioplankton is a major biological component involved in the degradation and mineralization of dissolved organic matter in Arctic waters (Ortega-Retuerta et al., 2012a).It can significantly affect the fate and distribution of organic matter within the entire water column (Bendsten et al., 2002) as well as the microbial food web activity through the assimilation of remineralized nitrogen.However, the contribution of Arctic bacterioplankton to plankton production in the context of Arctic warming remains unknown.Despite the fact that the AO basin now acts as a sink for atmospheric carbon dioxide (1996-2007;Manizza et al., 2013), the balance between autotrophy and heterotrophy may change in the future based on observations of enhanced stratification of the water column (Li et al., 2009), increased sea temperature (Timmermans et al., 2014;Steele et al., 2008), which acts as a key driver of Arctic bacterioplankton metabolism (Piontek et al., 2014;Bendsten et al., 2002), and changes in the riverine inputs of nutrients due to an increase in freshwater discharge (Shiklomanov and Lammers, 2011).In near-shore AO waters, riverine inputs already sustain part of the bacterial activity (e.g.Vallières et al., 2008).
Using a relatively simple biogeochemical modelling approach, Tank et al. ( 2012) shed light on the potential impact of riverine nutrient inputs on the PP of the AO.In the present study, we propose building on the static view provided by the work of Tank et al. (2012) by explicitly modelling the effect of the interactions between riverine dissolved organic nitrogen (RDON) and bacterioplankton.The objective is to use a pan-Arctic ocean-sea ice coupled model to quantify the contribution of usable RDON processed by marine bacterioplankton to the production of both bacterioplankton and phytoplankton in a scenario of melting sea ice over the period 1998-2011.

The physical model
We used the MIT general circulation model (MITgcm) (Marshall et al., 1997) coupled with a sea-ice model.The model is configured on a "cubed-sphere" grid encompassing the Arctic domain with open boundaries at ≈ 55 • N in the Atlantic and Pacific sectors.Prescribed boundary conditions for potential temperature, salinity, flow, and sea-surface elevation are provided from previous integrations of a global configuration of the same model (Menemenlis et al., 2005).The grid has a variable horizontal resolution with an average mesh of ∼ 18 km.The mesh resolves major Arctic straits, including many of the channels of the Canadian Archipelago.The sea-ice and fluid dynamics equations are solved on the same horizontal mesh.The 50-level vertical grid is height based, varying from 10 m thick near the surface to ∼ 450 m at a depth of ∼ 6 km.Bathymetry is derived from the US National Geophysical Data Center (NGDC) two-minute global relief data set (ETOPO2), which uses the International Bathymetric Chart of the Arctic Ocean (IBCAO) product for Arctic bathymetry (Jakobsson et al., 2008).The ETOPO2 data are smoothed to the model's horizontal mesh and mapped to the ocean's vertical levels using a "lopped cell" strategy (Adcroft et al., 1997), which permits an accurate representation of the ocean bottom boundary.
The ocean model's hydrography is initialized with observations taken from the Polar Science Center Hydrographic Climatology (PHC) 3.0 database (Steele et al., 2001).Initial sea-ice distributions are taken from the pan-Arctic Ice-Ocean Modeling and Assimilation System data sets (Zhang and Rothrock, 2003).Atmospheric forcings (10 m surface winds, 2 m air temperature and humidity, and downward longwave and shortwave radiation) are taken from the 6-hourly data sets of the Japanese 25-year ReAnalysis (JRA-25; Onogi et al., 2007).Monthly mean estuarine fluxes of freshwater are based on the Arctic Runoff database (Lammers et al., 2001;Shiklomanov et al., 2000).The sea-ice component of the coupled model follows the viscous-plastic rheology formulation of Hibler (1979) with momentum equations solved implicitly on a C-grid (Arakawa and Lamb, 1977) using a procedure based on Zhang and Hibler (1997).Fluxes of momentum into ice due to the overlying atmospheric winds and momentum fluxes between sea ice and the ocean are calculated by solving for the momentum balance at each surface grid column (Hibler and Bryan, 1987).This model configuration was previously used to study the Arctic freshwater budget (Condron et al., 2009).Modelling studies of Condron et al. (2009) compared to observations by Serreze et al. (2006) concluded that this model configuration is able to realistically represent the freshwater budget of the AO, including the import and export of freshwater from the Bering and Fram straits and from the Canadian Archipelago.

The riverine DON (RDON) discharge
To realistically represent the RDON flux in the AO in our biogeochemical model, we follow the approach adopted by Manizza et al. (2009), which is based on seasonally explicit regression relationships.These relationships use covariations between water yield and dissolved organic carbon (DOC) concentrations in circum-Arctic rivers to define riverine DOC (RDOC) monthly averaged fluxes for 10 regions in the pan-Arctic domain.These regions are the Barents Sea, Kara Sea, Laptev Sea, East Siberian Sea, Chukchi Sea, Bering Strait, Beaufort Sea, Canadian Archipelago, Hudson Bay, and Hudson Strait using published watershed areas and seasonal water runoff (Lammers et al., 2001).The approach uses empirical relationships quantifying the covariation between discharge and RDOC to scale the Lam-mers et al. (2001) discharge estimates into estimates of RDOC export.Estimates of RDOC export for December-March, April-July, and August-November were divided into monthly bins according to measured distributions of RDOC export for those months in Arctic rivers.For each season, [RDOC]-discharge relationships were developed.North American and Eurasian rivers were considered separately.Data from the Yukon, Mackenzie, and Kuparuk rivers were used to define a runoff-[RDOC] relationship for drainage areas in North America, and data from the Ob', Yenisey, and Lena rivers were used to define a runoff-[RDOC] relationship for drainage areas in Eurasia.RDOC for the Yenisey, Ob', Lena, and Mackenzie were collected as part of the Pan-Arctic River Transport of Nutrients, Organic Matter, and Suspended Sediments (PARTNERS) project (McClelland et al., 2008).RDOC concentrations for the Kuparuk River were collected as part of the NSF Study of the Northern Alaska Coastal System (SNACS, http://www.arcus.org/arcss/snacs/index.php).In all cases, discharge data were acquired from ArcticRIMS (http://rims.unh.edu/).Recent sampling efforts on these rivers have provided exceptional seasonal coverage (McClelland et al., 2008) and the total annual discharge of RDOC in the model is 37.7 TgC yr −1 , which is consistent with the estimate of Raymond et al. (2007).To initialize the model, we used the three-dimensional RDOC field obtained from the 3-decade integration of the model by Manizza et al. (2009).After that time, RDOC distributions are relatively steady, because the flushing time for tracers through the surface waters of the basin is of the order of a decade.RDOC was converted into nitrogen currency (RDON) using a molar C : N ratio of 40 : 1 (Tank et al., 2012;Köhler et al., 2003).We assume that 15 % of the RDON entering the model at river grid cells is usable by bacterioplankton (e.g.Wickland et al., 2012).

The biogeochemical model
We couple to the MITgcm physical model a biogeochemical model that explicitly represents the plankton ecosystem dynamics.The biogeochemical model is improved from previous applications in sub-Arctic (Le Fouest et al., 2005, 2006) and Arctic waters (Le Fouest et al., 2011, 2013b).In the model, nitrogen is the currency and it includes 10 compartments (i.e.nine in the pelagic domain + RDON that couples the marine and terrestrial cycling of nitrogen), chosen according to the ecosystem structure observed in the AO.Phytoplankton is size-fractionated into large (> 5 µm) and small (< 5 µm) phytoplankton (LP and SP, respectively).These two compartments encompass the major phytoplankton groups relevant for plankton dynamics and biogeochemistry in the Arctic waters (e.g.Li et al., 2009;Coupel et al., 2012).The two zooplankton compartments represent large (LZ, mainly copepods) and small (SZ, protozooplankton) organisms.Dissolved inorganic nutrients are nitrate (NO 3 ) and ammonium (NH 4 ).Detrital (i.e.produced by the bio- ), large (> 5 µm) and small (< 5 µm) phytoplankton, large zooplankton, protozooplankton, bacterioplankton, detrital particulate and dissolved organic nitrogen (dPON and dDON, respectively), and usable riverine dissolved organic nitrogen (RDON).Green, red and blue arrows represent nutrient uptake, grazing, and nitrogen recycling, respectively.
geochemical model components) particulate and dissolved organic nitrogen (dPON and dDON, respectively) close the nitrogen cycle.Bacterioplankton (BACT) are explicitly represented following the model of Fasham et al. (1990).They grow on NH 4 , dDON, and on the usable fraction of RDON (see the Appendix for details).LP and SP growth depends on light, NO 3 and NH 4 availability according to Liebig's law of the minimum.LZ graze on LP and SZ, whereas SZ graze on SP and BACT.Fecal pellets and LP basal mortality fuel the dPON pool.The dDON pool is made of unassimilated nitrogen resulting from SZ grazing, SP, SZ and BACT basal mortality and dPON fragmentation.BACT release, LZ excretion and unassimilated nitrogen resulting from SZ grazing are the sources of NH 4 in the model.NH 4 is converted into NO 3 through the nitrification process.For phytoplankton, nitrogen is converted into carbon using the Redfield carbon to nitrogen (C : N) molar ratio (106 : 16; Redfield et al., 1963) and into Chl using variable C : Chl mass ratios computed according to a modified version of the phytoplankton photoacclimation model of Cloern et al. (1995).The plankton biogeochemical model (Fig. 1) is fully detailed in the Appendix.Differential equations are given in Table 1, whereas biological parameters are given in Table 2.
Nitrate data used for the model initialization are from the World Ocean Atlas 2005 (National Oceanographic Data Centre, 2006).LP and SP are assigned a constant field over the model grid (0.2 and 0.002 mmol N m −3 in the top eight layers and below, respectively; e.g.Sherr et al., 2003;Ducklow, 1999, Taniguchi, 1999).The same conditions are imposed for BACT (e.g.Sherr et al., 2003;Ducklow, 1999).
LZ and SZ are assigned a constant field over the model grid (0.1 and 0.001 mmol N m −3 in the top eight layers and below, respectively; e.g.Sherr et al., 2003, Taniguchi, 1999).
The same conditions are imposed a priori for dPON.A value of 1 mmol N m −3 of NH 4 (e.g.Kristiansen et al., 1994) and dDON (e.g.Charria et al., 2008) is imposed at each grid cell.Boundary conditions at the North Atlantic and North Pacific sectors are data from the World Ocean Atlas 2005 (NODC, 2006) for NO 3 , and null for the remaining nine biogeochemical tracers.Apart from RDON, there are no riverine inputs for the remaining nine biogeochemical tracers.

Coupled model integrations
The model is spun up by repeating the January 1980-December 1989 decade twice.It is thereafter initialized with the physical and biogeochemical fields obtained from 31 December 1989 to run the 1990-2011 time period.Two simulations are then carried out: without usable RDON removal by bacterioplankton (our control run, hereafter CTRL run) and with usable RDON removal by bacterioplankton (hereafter RIV run).The difference between the two simulations provides information on the potential impact of the interactions between bacterioplankton and usable RDON on bacterioplankton production (BP), nutrient regeneration, and ultimately primary production (PP) in the Arctic basin.

Primary production
Shelf seas influenced least by riverine inputs of RDON show comparable simulated annual rates of total PP in the CTRL and RIV runs (Fig. 2).In the Barents Sea, simulated PP averaged over 1998-2011 reaches up to ∼ 80 gC m −2 yr −1 , in line with previous remote sensing estimates (up to 70-80 gC m −2 yr −1 on average over 1998-2010;Bélanger et al., 2013).In the central Chukchi Sea, simulated PP lies between 50 and 80 gC m −2 yr −1 , in agreement with the observed range (15-80 gC m −2 yr −1 on average over 1998-2007;Bélanger et al., 2013).The largest differences in total PP between the two runs are found in the river-influenced Eurasian seas (East Siberian Sea, Laptev Sea, and Kara Sea; Fig. 2).In the CTRL run, maximum simulated PP rates reach ∼ 30 gC m −2 yr −1 , which is more than 3-fold lower than satellite-derived and in situ estimates that can exceed 100 gC m −2 yr −1 (Bélanger et al., 2013;Codispoti et al., 2013;Sakshaug, 2004).In contrast, PP rates simulated in the RIV run (80-90 gC m −2 yr −1 ) show a better agreement with observations.
The increase in the 1998-2011 averaged annual PP in the RIV run relative to the CTRL run is due to the increase in NH 4 -supported PP (Fig. 3d, e and f).In contrast, overall, new PP remains unaffected by the bacterial use of RDON (Fig. 3a, b and c).In the Kara Sea, Laptev  Direct estimates of NH 4 -supported PP based on measurements are rare in Arctic coastal waters.Nevertheless, they can be crudely derived by subtracting new PP from total PP estimated by Sakshaug (2004).In the Eurasian and North American shelves, NH 4 -supported PP in the CTRL run is < 10 gC m −2 yr −1 (Fig. 3d) overall.This is low relative to the rates derived from Sakshaug's data, which would range from ∼ 25 to ∼ 40 gC m −2 yr −1 .By contrast, in the RIV run, rates simulated in offshore shelf waters are ∼ 10-30 gC m −2 yr −1 .However, closer to the coast, local rates reach 40-50 (Laptev Sea) and 70-80 gC m −2 yr −1 (Kara Sea; Fig. 3e).
Averaged over 1998-2011, the total PP simulated by the model and integrated over the whole AO is 662 ± 91 TgC yr −1 in the CTRL run and 717 ± 95 TgC yr −1 in the RIV run.These values are within the range of previously reported rates based on remote sensing or in situ data (385-1008 TgC yr −1 , Bélanger et al., 2013;Codispoti et al., 2013;Hill et al., 2013;Arrigo and van Dijken, 2011).Between the two model runs, the annual total PP increased by ∼ 8 %, on average, between 1998 and 2011.In September-October, when the simulated sea-ice concentration reaches its seasonal minimum, the annual total PP increase is more than twice this value (∼ 18 %, on average).

Bacterioplankton activity
The PP increase is tightly linked to a higher bacterioplankton activity that promotes RDON recycling into nutrients usable by both phytoplankton and bacterioplankton.The bacterioplankton biomass (BB), integrated between the sea surface and 50 m and averaged over April-June (spring) and July-September (summer), is shown in Fig. 4. As for PP, the Barents and Chukchi seas show comparable levels of BB in CTRL and RIV runs.In the Chukchi Sea, the BB simulated in spring (< 100-250 mgC m −2 ; Fig. 4a and b) overlaps with the measured range (222-358 mgC m −2 ; Kirchman et al., 2009).It is similar in summer, when simulated (∼ 100 mgC m −2 to > 800 mgC m −2 ; Fig. 4d and e) and measured BB levels (250-507 mgC m −2 ; Kirchman et al., 2009;Steward et al., 1996) are higher than in spring.In the Barents Sea, the simulated BB increases from < 100 mgC m −2 in spring to < 250 mgC m −2 in summer, falling within the measured range (from ∼ 80 mgC m −2 in spring to ∼ 400 mgC m −2 in summer, on average; Sturluson et al., 2008).In the highly river-influenced shelf seas, the two runs show notable differences in their simulated BB (Fig. 4c and f).In the central part of the Kara Sea, influenced by the Ob' and Yenisey river plumes, BB measured in late summer along a south-north transect from the Yama Peninsula to the Novaya Zemlya island is reported to range from ∼ 0.1 to 7 mgC m −3 (Sazhin et al., 2010).For the same time period and along a comparable transect, simulated values of BB are < 2 mgC m −2 in the CTRL run.However, in the RIV run, BB increases up to ∼ 6-7 mgC m −3 to match the values measured by Sazhin et al. (2010).
The depth-integrated (0-50 m) bacterioplankton production (BP) simulated in both the CTRL and RIV runs in summer in the Chukchi Sea (< 280 mgC m −2 d −1 ) is consistent with measurements reported for the same season (5-301 mgC m −2 d −1 ; Kirchman et al., 2009;Rich et al., 1997; Figure 3. Mean annual new primary production (gC m −2 ; upper panels) and NH 4 -supported primary production (gC m −2 ; lower panels) over 1998-2011 simulated in the CTRL run (left panels a and d) and the RIV run (middle panels b and e).Right panels (c and f) provide the absolute difference (gC m −2 ; RIV run -CTRL run).Steward et al., 1996).In the Beaufort Sea, influenced by the Mackenzie River plume, simulated BP is lower than ∼ 6 mgC m −2 d −1 in the CTRL run, which is far below measurements made within the area (25-68 mgC m −2 d −1 ; Ortega- Retuerta et al., 2012a;Vallières et al., 2008).By contrast, in the RIV run, simulated BP (< 30 mgC m −2 d −1 ) approaches the lower range of observations.Similarly, BP simulated in the CTRL run for the Kara Sea (< 30 mgC m −2 d −1 ) does not exceed the first mid-range of measurements given by Meon and Amon (2004; 12-79 mgC m −2 d −1 ).In the RIV run, the simulated BP (∼ 4-90 mgC m −2 d −1 ) overlaps the measured range (12-79 mgC m −2 d −1 ; Meon and Amon, 2004) to reach up to 120 mgC m −2 d −1 locally.This result is consistent with enrichment experiments conducted with surface oceanic water sampled in the Beaufort Sea that showed a 43 % increase in BP when Mackenzie River water was included in samples (see Ortega-Retuerta et al., 2012a).
Averaged over 1998-2011, the total BP simulated by the model and integrated over the whole AO is, on average, 26 % higher in the RIV run (68 ± 9 TgC yr −1 ) than in the CTRL run (54 ± 8 TgC yr −1 ).Bacterioplankton recycle RDON into nutrients that can be used by both phytoplankton and bacterioplankton, hence promoting their growth.In addition, bacterioplankton and small phytoplankton are grazed by microzooplankton that, in turn, are grazed by mesozooplankton.More organic matter is channelled towards the upper trophic levels, a flow that also contributes to fuelling the dDON and NH 4 pools through recycling.By enabling the removal of RDON by bacterioplankton in the biogeochemical model, the biomass of microzooplankton and mesozooplankton, averaged over 1998-2011, increased by ∼ 16.1 and 43.6 %, respectively.

The bacterioplankton production versus primary production ratio (BP : PP)
The BP : PP ratio is computed over the AO shelf, delimited here by the 200 m isobaths, for ice-free waters (i.e. with less than 30 % ice cover).On average for the 1998-2011 period, the simulated BP : PP ratio is 0.19 ± 0.02 in the CTRL run and 0.21 ± 0.01 in the RIV run.These values lie within the range observed in open (0.02; Kirchman et al., 2009) and coastal (0.37-0.43;Ortega-Retuerta et al., 2012a; Garneau et al., 2008) waters.When looking at the temporal evolution of BP : PP in the RIV run (Fig. 5), the model simulates a significant increase in PP (r = 0.57, p < 0.05) and BP (r = 0.63, p < 0.05) between 1998 and 2011, with a production maximum simulated in 2007, the year showing the higher seaice minimum.However, there is no evidence of a significant temporal trend of BP : PP (r = −0.09,p > 0.05) over 1998-2011.This result suggests that, with a constant annual flux of RDON into the coastal AO, the significant increase in simulated BP in the model is not high enough to promote a higher contribution of heterotrophy with respect to autotrophy within the upper water column.

Discussion
The coupled model suggests that NH 4 produced from the remineralization of RDON by the microbial food web contributed ∼ 8 % to annual pan-Arctic PP over the 1998-2011 period.This is about twice the value given in the study by Tank et al. (2012) that, in addition to RDON, accounted for the contribution of riverine inorganic nutrients as well as of the photochemical transformation of RDON into NH 4 .In our coupled model, the uptake of RDON by marine bacterioplankton and its subsequent recycling into reduced nitrogen is the sink term that shapes, with ocean transport, the spatial and temporal distribution of RDON.The photoammonification process is not parameterized but, if so, it would fuel the stock of NH 4 available for phytoplankton and bacterioplankton use, particularly in summer (e.g.Le Fouest et al., 2013;Xie et al., 2012).The RDON contribution to plankton production simulated by the coupled model can thus be considered as a minimum estimate.
From the total input of RDON, only a fraction is directly usable by the plankton (e.g.Wickland et al., 2012).The fraction that enters the coupled model by the 10 river source points is set to 15 % of the total RDON input according to a study by Wickland et al. (2012), which suggests that about 15 % of the total RDON pool can be degraded within less than 1 month.This value was chosen based on annual averages calculated from measurements or from model outputs for the Mackenzie River, Yukon River, Kolyma River, Lena River, Yenisey River, and Ob' River (e.g.Wickland et al., 2012).Note, however, that the average values given in Wickland et al. (2012) vary between seasons and rivers.They are lowest in the Lena River (8 %) and highest in the Ob' River (19 %).Maximum values as high as 24 % of usable RDON are reported for the Ob' River.Sensitivity analyses with different parameterizations of the usable RDON fraction set amongst rivers and seasons would hence be informative on the amplitude of the PP and BP response to spatial and temporal variations of the usable RDON flux.To be robust, they should be combined with sensitivity analyses of the freshwater discharge to better constrain the RDON flux.In the Mackenzie River, strong inter-annual variations in terms of peak of discharge and maximum spring flow were observed in the last four decades (Yang et al., 2015).Nevertheless, the use a constant fraction of usable RDON as preformed in the present study provides a first-order estimation of its contribution to BP and PP that is consistent with the current state of knowledge about the RDON inputs.In addition to the usable RDON flux into coastal oceans, autochthonous sources of DONl (usable RDON + dDON) are important in fuelling BP.Despite improved BP estimates simulated in the RIV run, the rates remain within the lower range of the observations.It can result from unresolved sources of DONl within the model such as ice-edge and under-ice phytoplankton blooms (Arrigo et al., 2012;Perrette et al., 2011), and from miss-ing biological processes like sloppy mesozooplankton feeding and viral lysis.
In the biogeochemical model, the usable RDON, dDON, and NH 4 produced by the plankton components are taken up by bacterioplankton to build up biomass.The synthesis of cell proteins requires at least carbon and nitrogen.Bacterioplankton obtain all their carbon and some of their nitrogen from DONl (usable RDON + dDON).The simulated NH 4 uptake supplements their nitrogen requirements.The growth function is formulated using the Fasham et al. (1990) model.It assumes, in a balanced growth situation, where N and C assimilation occurs simultaneously and where bacterioplankton have fixed stoichiometry, that the ratio of NH 4 uptake to DONl uptake is constant (0.6; see Appendix A) to ensure that the biomass of the required C : N ratio is produced from DONl with a given C : N ratio.If there is not enough NH 4 available, the uptake rate of both DONl and NH 4 decreases, allowing both N and energy limitation.In Arctic waters, the inhibition of DOC uptake by bacterioplankton under inorganic nitrogen limitation was shown by Thingstad et al. (2008).However, as DONl in the model is made a proxy of DOC, the C : N ratio of the substrate is assumed constant.As a consequence, any explicit stoichiometric treatment of the simulated bacterioplankton metabolism is precluded as well as any stoichiometric coupling between DOC and inorganic nutrients (e.g.Thingstad et al., 2008).In addition, the implicit treatment of DOC in the model implies that all of the DOC required for growth is in N-containing forms.Hence it assumes that bacterioplankton cannot be N-limited in substrate.However, N-limitation of bacterioplankton production was observed in summer in surface waters of the Beaufort Sea (Ortega-Retuerta et al., 2012b).This pattern contrasts with the organic carbon limitation observed in the Yenisei and Mackenzie river plumes and the adjacent Kara and Beaufort seas (Meon and Amon, 2004;Vallières et al., 2008), hence highlighting the difficulty in drawing a general pattern on the AO scale.Nevertheless, making the C : N ratio of substrates of terrigeneous and marine origin vary in a realistic way in biogeochemical models would further be required.Single explicit pools of DOC and DON represented as two different state variables, as well as a distinction between readily usable molecules (turnover within days) and more complex ones (turnover within a month) would also make the model more realistic.The parameterization of variable C : N ratios is not trivial as it requires large in situ data sets (see Letscher et al., 2015) and, in Arctic river-influenced shelf seas, a good knowledge of the characteristics of the terrigenous dissolved organic matter flowing into the coastal ocean (e.g.Mann et al., 2012).Appropriate values for the maximum uptake rates and half-saturation constants may not be easily obtained from existing data in the Arctic.As a result, the coupled model that is used in the present study is an interesting compromise relative to more complex (in terms of number of biological equations and parameters) models of www.biogeosciences.net/12/3385/2015/Biogeosciences, 12, 3385-3402, 2015 bacterioplankton growth applied to shelf waters (e.g.Auger et al., 2011;Anderson and Williams, 1998).
In the model, bacterioplankton compete with phytoplankton for the NH 4 remineralized from the usable RDON and dDON pools.This competition for a nutrient resource acts as a bottom-up control of the simulated phytoplankton and bacterioplankton production and, finally, of the BP : PP ratio.For bacterioplankton, the maximum growth rate is temperature normalized.At 10 • C, which corresponds to a sea surface temperature within the upper range of observations over the shelves in summer, it takes a value of 3 d −1 .Hence the NH 4 uptake efficiency (α = maximum growth rate/half-saturation constant for uptake) can reach 30 m 3 mmol N −1 d −1 , which is about 2 times higher than for small phytoplankton (14 m 3 mmol N −1 d −1 ).In their study, Lignell et al. (2013) report values of α that are about 10 times higher for bacterioplankton than for small phytoplankton, hence a larger difference compared to the model.Nevertheless, the difference in α in the model is comparable to that estimated for Isefjord at the entrance of the Baltic Sea (see Lignell et al., 2013).Although the model parameterization can be improved in that respect, the model is in fair agreement with the theoretical and empirical results showing that smaller cells are more efficient in nutrient uptake than larger ones (Lignell et al., 2013).In contrast to bacterioplankton, phytoplankton uptake of inorganic nutrients is also limited by light.In the model, the diffuse attenuation of the incident light caused by the pool of coloured dissolved organic matter (0.05 m −1 ) is set as constant in the model.This results in the light attenuation in the water column being the same in river plumes as in open and clearer waters.However, river plumes transfer to the coastal marine environment large amounts of optically active coloured dissolved organic matter of terrigenous origin that strongly attenuate the incident light propagation with depth.As the model does not account for the stronger light attenuation in river plumes, it may overestimate the simulated phytoplankton growth on NH 4 recycled from RDON by bacterioplankton and underestimate the BP in river plumes.As a consequence, the spatial and temporal evolution of the simulated BP : PP ratio can be impacted on shelves.In addition, the ability of Arctic phytoplankton to assimilate low molecular weight DON compounds (50 % of total nitrogen assimilated annually; see Simpson et al., 2013) is likely to also play an important role in the phytoplankton-bacterioplankton competition on shelves.A more accurate representation of the simulated underwater light field and uptake of nutrients in river plumes in the coupled model will certainly improve its ability to simulate the competition for nutrients between phytoplankton and bacterioplankton, and hence predict the temporal evolution of the BP : PP ratio within Arctic waters.

Conclusions
A pan-Arctic physical-biogeochemical model was used to quantify the contribution of usable dissolved organic nitrogen drained by the major pan-Arctic rivers to marine bacterioplankton and phytoplankton production in a scenario of melting sea ice (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011).By accounting for the removal of RDON by bacterioplankton in the coupled model, the ability to predict PP and BP in river-influenced shelves is improved.The key points of the study are that 1. on average between 1998 and 2011, the removal of usable RDON by bacterioplankton is responsible for an increase of ∼ 26 % in the annual BP, and an increase of ∼ 8 % in the total annual PP; 2. recycled ammonium is responsible for the total PP increase; total summertime PP is increased by ∼ 18 %, on average, over 1998-2011; and that 3. the processing of usable RDON by bacterioplankton promotes a higher annual BP and PP, but there is no significant temporal trend in the BP : PP ratio over 1998-2011 on the ice-free shelves; this suggests no significant evolution in the balance between autotrophy and heterotrophy in the last decade, with a constant annual flux of RDON into the coastal ocean.
The effect of the predicted warming on the Arctic watersheds is linked to a potential regional increase in RDON inputs into the AO shelf by 32-53 % before the end of the century (Frey et al., 2007).Combined with the accelerated sea-ice decline (Comiso et al., 2008) and an increase in seawater temperature on Arctic shelves (Timmermans et al., 2014), this new biogeochemical and physical setting might exacerbate the competing effect for resources between autotrophs and heterotrophs as sea ice recedes in summer.As a consequence, the metabolic state of the AO shelves could be altered.Nevertheless, to obtain robust predictions of the response of the microbial food web functioning and mass fluxes, coupled models would require improvements in parameterized landocean fluxes in terms of spatial and temporal variability of freshwater discharge and nutrient fluxes.In their study combining in situ data sets and modelling, Holmes et al. (2011) show that annual fluxes of RDOC in the Lena River, estimated between 1999 and 2008, can vary by about a factor of 2. Such variations accentuate the significance of considering the short-term and inter-annual variability of the continental fluxes into the coastal ocean when deriving temporal trends in plankton production and investigating potential changes in trends related to the Arctic warming.Finally, model predictions of future trajectories of PP (e.g.Vancoppenolle et al., 2013) would probably benefit from considering riverine nutrient fluxes as important drivers of PP on Arctic shelves in future decades.The set of differential equations that include the mechanistic formulations cited below is given in Table 1.The biological parameters related to the mathematical equations are detailed in Table 2.

A1 Phytoplankton
The growth rate (µ LP,SP , d −1 ) of large and small phytoplankton (LP and SP, respectively) depends on both light and nitrogen availability.It is computed according to Liebig's law of the minimum between the nutrient-based and light-based growth rates (µ LP, SP N and µ LP, SP light , respectively): The nutrient-based growth rate is computed as follows: where µ LP, SP max is the maximum growth rate and lim LP, SP N the total nutrient limitation term (dimensionless) computed according to the substitutable model of O'Neill et al. (1989): where lim LP, SP NO 3 and lim LP, SP NH 4 are the nitrate (NO 3 ) and ammonium (NH 4 ) uptake fractions (dimensionless), respectively.K LP, SP NH 4 and K LP, SP NO 3 are the half-saturation constants for NH 4 and NO 3 uptake, respectively.NH 4 is set to be the preferred inorganic nitrogen source (Dorch, 1990) with a higher affinity for SP (Tremblay et al., 2000).This is expressed in the model by half-saturation constants for NH 4 uptake (K LP, SP NH 4 ) lower than for NO 3 that, when used with the substitutable model, allow for an inhibitory effect of NH 4 on NO 3 uptake as often observed (Dorch, 1990).It implies that NO 3 uptake by LP and SP is inhibited by NH 4 at concentrations 2-fold and 10-fold lower than NO 3 concentrations, respectively where C : Chl is the carbon to Chl ratio (g g −1 ) and α LP, SP is the initial slope (mg C (mg Chl) −1 (Ein m −2 d −1 ) −1 ) of the photosynthesis-irradiance curve.Photoacclimation translates the adaptative response through varying C : Chl ratios in response to light and nutrient availability (e.g.Cloern et al., 1995;Geider et al., 1997;MacIntyre et al., 2002).
Varying C : Chl ratios are computed using a modified version of the empirical relationship of Cloern et al. (1995) successfully applied to Hudson Bay in the Arctic (Sibert et al., 2011).The ratios can vary up to 4-to 6-fold based on the general photoacclimation rule given by MacIntyre et al. (2002) and on Arctic nano-and pico-phytoplankton data (DuRand et al., 2002;Sherr et al., 2003) as follows: where K LP, SP E is the half-saturation parameter driving the curvature of the C : Chl versus light relationship.E z (Ein m −2 d −1 ) is the downwelling photosynthetically active radiation (PAR) propagating within the water column according to the Beer-Lambert law: where PAR0 is the PAR just below the sea surface.The diffuse attenuation of PAR with depth (z) is due to the simulated Chl (kchl) (m −1 ; Morel, 1988), water molecules (kw; 0.04 m −1 ; Morel, 1988), and non-chlorophyllous matter (knonchl).knonchl is set to 0.05 m −1 .kchl is calculated according to Morel et al. (1988)

A2 Zooplankton
Mathematical formulations and parameters related to large zooplankton (LZ) dynamics were chosen to primarily reflect mesozooplankton.Grazing (G LZ , d −1 ) is described by an Ivlev function: LZ graze upon LP and protozooplankton (SZ) with a preyspecific grazing rate assumed to be proportional to the relative biomass of the prey (Campbell et al., 2009) defined for LP as follows: Losses in LZ biomass are due to NH 4 release, fecal pellet production (non-assimilated nitrogen ingested), and mortality.Mortality encompasses senescence and predation (Eiane et al., 2002).It is described by a density-dependent quadratic function.It implicitly represents cannibalism as well as predation by macrozooplankton (Forest et al., 2012;Berline et al., 2008) and limits the occurrence of oscillations generated in such non-linear systems (Edwards and Bees, 2001).The constant of mortality is set to 0.2 (mmol N m −3 ) −1 to simulate realistic mortality rates (e.g.Ohman et al., 2004).SZ grazing (G SZ ) upon SP and bacterioplankton (BACT) is formulated by a sigmoidal Holling-type-III function: where G max SZ and K G are the maximum grazing rate (d −1 ) and the half-saturation constant for grazing (mmol N m −3 ), respectively.The grazing function provides a threshold-like limit for low SP biomass that enhances the biological system stability (e.g.Steele and Henderson, 1992).In polar waters, there is evidence that protozooplankton only exert a control on small phytoplankton biomass beyond a threshold (Lancelot et al., 1997).As for LZ, SZ graze upon both SP and BACT, with a prey-specific grazing rate (d −1 ) assumed to be proportional to the relative biomass of the prey defined for SP as follows: We set the fraction of food ingested and converted into biomass to 30 % (Straile, 1997).Lehrter et al. (1999) report that > 30 % of the total nitrogen release by SZ could be in the dissolved organic form.In the model, assuming that 40 % is released as detrital DON (dDON), the 60 % remaining are lost as NH 4 .Other SZ loss terms are grazing by LZ and mortality.Similarly to LZ, SZ mortality is expressed by a density-dependent quadratic function to encompass grazing amongst SZ.

A3 Bacterioplankton
Bacterioplankton is explicitly simulated following the model of Fasham et al. (1990).They grow on NH 4 , dDON and usable RDON.Usable RDON is considered as 15 % of total RDON (e.g.Wickland et al., 2012) and is converted into N currency (RDON) using a C : N ratio of 40 (Tank et al., 2012;Köhler et al., 2003).Dissolved organic matter (DOM) is a complex bacterial substrate representing a source of nitrogen (DON) and carbon (DOC).As nitrogen is the sole currency of the model, the simulated DONl is made a proxy of DOC for bacterioplankton uptake.This means that bacterioplankton in the model obtain all of their carbon and some of their nitrogen from the usable fraction of RDON and from detrital DON (dDON).This assumes that all of the DOC required for growth is in N-containing forms.By contrast, the simulated ammonium uptake supplements the bacterioplankton N requirements for growth.The DONl uptake rate (Ubact DONl , d −1 ) is represented by a Michaelis-Menten model: Ubact DONl = Ubact max BACT DONl where K BACT DONl is the half-saturation constant for DONl uptake (mmol N m −3 ), and S the total nitrogenous substrate (mmol N m −3 ) defined as S = (NH 4 , 0.6DONl). (A20) According to the study by Bendtsen et al. (2002) in the Greenland Sea, a Q 10 function was introduced using a Q 10 factor of 3 (Kirchman et al., 2005): where T is the simulated seawater temperature.A temperature normalized maximum uptake rate (Ubact max ) of 1 d −1 was used to simulate maximum growth rates in line with those measured in polar waters (e.g.Nedwell and Rutter, 1994).A growth efficiency of 20 % (Ortega-Retuerta et al., 2012a;Meon and Amon, 2004) was imposed.The fractioning of the two DONl pools (i.e.usable RDON and dDON) is set as follows: Similarly, the uptake rate of NH 4 (Ubact NH 4 , d −1 ) is represented as follows: where K BACT NH 4 is the half-saturation constant for NH 4 uptake (mmol N m −3 ).This formulation ensures that the uptake of NH 4 will be 0.6-fold the uptake of DONl, as required by the balanced growth model (e.g.Fasham et al., 1990).BACT V. Le Fouest et al.: Modelling the impact of riverine DON removal by marine bacterioplankton cannot grow on NO 3 in the model.NO 3 uptake is an energetically expensive process, so that bacterioplankton usually accounts for more ammonium than nitrate uptake (Lipschultz, 1995).Furthermore, although a substantial nitrate uptake by bacterioplankton was reported at high latitudes, it occurred under very specific conditions such as in high-nitrate lowchlorophyll waters (Kirchman and Wheeler, 1998) or in lowchlorophyll waters dominated by cyanobacteria (Fouilland et al., 2007).Such conditions are not achieved in the model.Senescence is in the NH 4 form, and it represents 5 % of the biomass.

A4 Detritus
The pool of detrital particulate organic nitrogen (dPON) is fuelled by the production of fecal pellets by LZ, and by LZ and LP mortality.The sedimentation loss term (d −1 ) is expressed as a quadratic function allowing for the implicit increasing aggregation of particles with increasing dPON concentration (see Guidi et al., 2008): where sed_dpon is the sedimentation constant (m d −1 (mmol N m −3 ) −1 ).The second loss term is the dPON fragmentation into dDON (e.g.Grossart and Ploug, 2001).
The dDON pool results from dPON fragmentation, SP, SZ, and BACT mortality and SZ release.It is explicitly remineralized into NH 4 by BACT.

A5 Nutrients
NH 4 resulting from the remineralization by BACT and from the LZ and SZ release fuels regenerated primary production and BACT production.In turn, NH 4 undergoes nitrification (d −1 ) into NO 3 as follows: where nitrif max is the maximum nitrification rate (d −1 ), and K N nitrif and K light nitrif the half-saturation constants for NH 4 (mmol N m −3 ) and light (Ein m −2 d −1 ) use, respectively.The latter is defined as a fraction of surface PAR (E 0 ) as follows: Mean values taken from the literature (Guerrero and Jones, 1996;Olson, 1981aOlson, , 1981b) ) are used to set parameters.

Figure 4 .
Figure 4. Seasonal climatology of the 0-50 m integrated bacterial biomass (mmolN m −2 ) for spring (upper panels) and summer (lower panels) over the 1998-2011 period simulated in the CTRL run (left panels a and d) and in the RIV run (middle panels b and e).Right panels (c and f) provide the absolute difference (gC m −2 ; RIV run -CTRL run).

Figure 5 .
Figure 5.Time course of primary production (PP, TgC yr −1 ) (top panel), bacterioplankton production (BP, TgC yr −1 ) (middle panel), and the BP : PP ratio in the ice-free shelves (see text for details) of the Arctic Ocean domain (> 66.5 • N) simulated in the RIV run.The dashed straight lines represent the linear trend computed over the 1998-2011 period.
. The equation used to compute the light-based growth rate is µ LP,SP light = µ LP,SP max lim LP,SP light , (A6) where lim LP, SP light is the light limitation term (dimensionless) expressed as lim LP,SP light = 1 − e