The Gulf of St. Lawrence Biogeochemical Model: A Modelling Tool for Fisheries and Ocean Management

The goal of this paper is to give a detailed description of the coupled physical-biogeochemical model of the Gulf of St. Lawrence that includes dissolved oxygen and carbonate system components, as well as a detailed analysis of the riverine contribution for different nitrogen and carbonate system components. A particular attention was paid to the representation of the microbial loop in order to maintain the appropriate level of the different biogeochemical components within the system over long term simulations. The skill of the model is demonstrated using in situ data, satellite data and estimated fluxes from different studies based on observational data. The model reproduces the main features of the system such as the phytoplankton bloom, hypoxic areas, pH and calcium carbonate saturation states. The model also reproduces well the estimated transport of nitrate from one region to the other. We revisited previous estimates of the riverine nutrient contribution to surface nitrate in the Lower St. Lawrence Estuary using the model. We also explain the mechanisms that lead to high ammonium concentrations, low dissolved oxygen, and undersaturated calcium carbonate conditions on the Magdalen Shallows.


INTRODUCTION
The Gulf of St. Lawrence (GSL) is a large semi-enclosed sea (covering an area of about 240000 km 2 ) located in eastern Canada that connects to the northwest Atlantic through the relatively shallow Strait of Belle Isle (50-80 m; Bailey, 1958) and the deep (480 m) Cabot Strait (Figure 1). Waters flowing in through the Strait of Belle Isle are mostly of Arctic origin [see Figure 1 of Lavoie et al. (2019)], while waters entering at Cabot Strait are a mixture of cold water from the Labrador Current and warmer North Atlantic Central Water. The GSL is characterized by deep channels (Laurentian, Esquiman, and Anticosti) bordered by relatively narrow shelves, except for the Magdalen Shallows, a vast shelf that covers most of the southern GSL with an average depth of 60-65 m. The GSL is also characterized by an important inflow of freshwater. The St. Lawrence hydrographic system is one of the largest in the world and the third largest in North America, after the Mississippi and Mackenzie rivers 1 , with 80% of the runoff flowing into the Lower St. Lawrence Estuary (LSLE, Lavoie et al., 2017). This runoff generates an estuarine-like circulation, with fresher water moving out of the system at the surface, and deep water in the channels moving up to the channels head. The Laurentian Channel remains deep up to its head, where it rises rapidly towards shallow sills (up to 50 m) near the Saguenay Fjord. This particular topography leads to upwelling of deep water under the action of tides, with subsequent mixing of the deep water with surface water (Ingram, 1983;Saucier et al., 2003). These deep waters are rich in nutrients, and with the additional contribution of riverine nutrient and organic matter fluxes (Hudon et al., 2017), it leads to high primary production in the LSLE (Therriault and Lacroix, 1976;Greisman and Ingram, 1977;Therriault and Levasseur, 1985;Cyr et al., 2015) and northwest GSL (nwGSL). The latter is separated from the LSLE by natural topographic features to the west (corresponding to the Pointe-des-Monts transect in Figure 1) and by Anticosti Island to the east. It is characterized by a cyclonic gyre and a strong current along the southern shore generated by the LSLE outflow (the Gaspé current). Primary production is generally low in the center of the gyre while high rates can be observed in the Gaspé Current (e.g., Sévigny et al., 1979;de Lafontaine et al., 1991).
The shallow southern GSL, represented in large part by the Magdalen Shallows, receives part of the nutrients carried by the relatively fresh surface waters, although the importance of this nutrient source for local primary production is uncertain (Sinclair et al., 1976;de Lafontaine et al., 1991;Levasseur et al., 1992;Cyr et al., 2015). The northeast GSL, on the other hand, is less subject to the direct influence of freshwater runoff. The northeast is less stratified than the other regions of the GSL and under the direct influence of water flowing in from the Labrador Shelf and Atlantic Ocean, which have lower nutrient concentrations compared to the LSLE (Starr et al., 2004). Differences in nutrient sources lead to heterogeneous conditions in primary production in the GSL (de Lafontaine et al., 1991 and references therein;Lavoie et al., 2008). The geographic location of the GSL, in combination with its complex topography and circulation, also lead to the aggregation of large quantities of krill (Simard and Lavoie, 1999;Maps et al., 2015;Lavoie et al., 2016) and copepods (Maps et al., 2011;Brennan et al., 2021) that are sought as food by endangered whales (blue whales (Lesage et al., 2018), North Atlantic right whales (Simard et al., 2019) that visit the GSL for many months each year.
The deep waters of the GSL also present particularly low concentrations of dissolved oxygen (DO), pH (∼7.5, total scale) and calcium carbonate saturation states (below unity for both aragonite and calcite), especially in the LSLE (Gilbert et al., 2005;Mucci et al., 2011;DFO, 2020). At depths below about 200 m, DO concentrations decrease from about 180 mmol m −3 at Cabot Strait to 60 mmol m −3 at the head of the Laurentian Channel (Gilbert et al., 2005). This decrease in DO is accompanied by an increase in DIC (Mucci et al., 2011). Nitrate and silicate concentrations, on the other hand, increase by about 2 mmol N m −3 and 15 mmol Si m −3 respectively (Coote and Yeats, 1979), as the water moves up towards the head of the Laurentian Channel. Two main factors are responsible for these gradients: (1) the proportion of Labrador Current and North Atlantic Central waters in the mixture that enters the Laurentian Channel, since these water masses have different chemical properties, and (2) the sinking and remineralization of organic matter produced in the surface layer by primary producers and supplied by the rivers (Gilbert et al., 2005;Claret et al., 2018;Lavoie et al., 2019). These processes lead to loss of oxygen and acidification in the deep water of the LSLE where ventilation is weak (Gilbert et al., 2005;Mucci et al., 2011). Low pH (<7.7, total scale) and calcium carbonate saturation states below unity (even for calcite over a few months in summer or fall) are also found in the shallower southern GSL (Menu-Courey et al., 2019;DFO, 2020).
Substantial interannual variability in physical, biological and chemical conditions is also present in the system. For example, the sea-ice cover in the GSL, whose demise controls the initiation of the phytoplankton bloom in many regions (Le Fouest et al., 2005), is highly variable from one year to the next; record high and low volumes were observed within a few years in the last decades . Large anomalies in river runoff and wind forcing, which affect circulation patterns and strength in the GSL (Saucier et al., 2009;Lavoie et al., 2016;Lavoie et al., 2017) have been observed locally , and in the northwest Atlantic (Hu et al., 2011;Peterson et al., 2017). These large interannual variations in forcing result in large primary and secondary production variability (Plourde et al., 2014), and have an impact on the properties of the water masses that enter the GSL.
The GSL is also impacted by climate variability and change. Long term trends in physical and chemical conditions have been recorded in the GSL, such as warming of surface and sub-surface waters Loder et al., 2013), deoxygenation (Gilbert et al., 2005;Claret et al., 2018;Jutras et al., 2020), and acidification (Mucci et al., 2011) of the deep Laurentian Channel waters. These trends are similar to those obtained with Earth System Models near the mouth of the Laurentian Channel and are projected to continue in the future (Lavoie et al., 2019). These changes have an impact on the health and distribution of the various marine organisms that live in these waters (Bianucci et al., 2016;Brennan et al., 2016;Stortini et al., 2017;Greenan et al., 2019), of which many are of commercial importance. Commercial fisheries totaled $3.0 billion in the DFO's Atlantic regions (Maritimes, Gulf, Quebec and Newfoundland and Labrador) in 2016 (Lobster Fishing Areas 27 -

-Integrated Fisheries Management
Plan 2 ). Within the GSL, American lobster and snow crab are the most important fisheries, especially in the southern Gulf, with shrimp and groundfish also being important in the deep channels of the northern GSL (DFO, 2011a,b;Gaudet and Leger, 2011). The GSL is also of great importance for Atlantic Salmon as it has many salmon rivers and salmon smolts spend a few months in the GSL in summer before migrating to the Labrador Sea (e.g., Ohashi and Sheng, 2018).
To better understand the functioning of the ecosystem and to project the evolution of primary production, deoxygenation and acidification within the GSL, the use of a high-resolution regional model is required.
Different NPZD (nutrient-phytoplankton-zooplanktondetritus) models have been implemented in the GSL. Zakardjian et al. (2000) implemented a 1D and a 2D model in the LSLE to test the importance of diatom sinking speed and horizontal advection on the observed delay of the phytoplankton spring bloom in that region. Tian et al. (2000) implemented a 1D model in the center of the nwGSL cyclonic gyre. They highlighted the importance of organic matter remineralization and nitrification taking place below the euphotic layer, but above the winter mixed layer depth, as a nutrient source to the surface layer, as well as the importance of the sinking of large particles for the export of carbon below the winter mixed layer. A slightly modified version of the latter model was also implemented in Bonne Bay, Newfoundland, to look at the processes leading to the relative constancy of export fluxes from spring to fall (Tian et al., 2001). Le Fouest et al. (2005), hereafter identified as LF05, then implemented a moderately complex NPZD model into a 3D ocean circulation model of the whole GSL in order to identify the physical forcing governing primary production in different regions of the GSL. Subsequent improvements were brought to the original model: light attenuation by colored dissolved organic matter (CDOM) as a function of surface salinity and phytoplankton photoadapation (Le Fouest et al., 2006. Mei et al. (2010) implemented CDOM explicitly, as a passive tracer, rather than as a function of salinity as in Le Fouest et al. (2006). They also added a temperature sensitive maximum growth rate for phytoplankton and zooplankton. Their goal was to look at the importance of these parameters to the timing of the phytoplankton bloom.
However, although the LF05 model reproduces some characteristics of the system, a drift of the biogeochemical conditions is present, limiting its use for multi-year simulations aimed at analysing interannual variability or climate change. In this paper, we describe GSBM (the Gulf of St. Lawrence Biogeochemical Model), which is based on LF05, but to which numerous improvements were made to better represent the microbial loop and maintain realistic biogeochemical conditions over the long term. The importance of properly representing and/or parameterizing the processes (sinking speed, fragmentation and remineralization rates) affecting particulate organic matter in numerical models to reproduce the biogeochemical tracers distribution in the marine environment has been reported many times (Fasham and Evans, 1995;Kawamiya et al., 1995Kawamiya et al., , 1997Haigh et al., 2001;Kriest et al., 2012;Kriest and Oschlies, 2015;Baker et al., 2017;Lavoie et al., 2019). GSBM also includes a DO module, a carbonate chemistry module, as well as the supply of organic and inorganic nutrients, DIC and alkalinity by rivers. GSBM is coupled to the ocean circulation model CANOPA. In the following sections, we describe the different components of the model and show its skill at reproducing the important biogeochemical features of the ecosystem.

Ice-Ocean Model CANOPA
The ocean circulation model used in our study, described in detail in Brickman and Drozdowski (2012) and with additional validation in Lavoie et al. (2016), is based on the Nucleus for European Modelling of the Ocean (NEMO) system. The modelling system is based on the Océan PArallélisé code, version 9.0 (OPA; Madec, 2012). A thermodynamic-dynamic sea-ice model, Louvain-la-Neuve Sea Ice Model (LIM2; Madec et al., 1998), is coupled to the circulation model. The grid of the model covers the GSL, the Scotian Shelf, and the Gulf of Maine (GoM; Figure 1). The grid has a horizontal resolution of 1/12 • in latitude and longitude (about 6 × 8 km depending on the location) and 46 layers of variable thickness (from 6 m close to the surface to 250 m at depth in the ocean).
It is a prognostic model, meaning that the temperature and salinity fields are free to evolve with time and are only constrained through open boundary conditions, freshwater runoff, and surface forcing. We use the Brickman and Drozdowski (2012) monthly climatologies of temperature and salinity to initialize the model and set the open boundary conditions. These were built from extensive Fisheries and Oceans Canada databases covering the 1971-2000 period. An annual cycle of the barotropic transport is prescribed at the Strait of Belle Isle in addition to a smaller baroclinic transport calculated from the monthly temperature and salinity fields. Tidal components (M2, S2, N2, O1, and K1) are included in the model through surface elevation and barotropic transport at the open boundaries.
Freshwater enters the domain through precipitation and runoff from the 78 main rivers. The monthly runoff of the St. Lawrence River is estimated from sea level measurements at Lauzon using the approach of Bourgault and Koutitonsky (1999) and obtained from the St. Lawrence Global Observatory 3 . Runoff for all other rivers is obtained from a hydrological model described in Lambert et al. (2013). The hydrological model provides natural runoff curves (i.e., following precipitation and evaporation). Some of the simulated seasonal runoff curves were modified to take into account the presence of large hydro-electric dams (Saguenay, Manicouagan, Betsiamites and Outardes rivers, see Lavoie et al., 2017). We use three-hourly surface forcing data (air temperature, relative humidity, winds, cloud cover, and precipitation) obtained from the Canadian Meteorological Centre -Global Environmental Multiscale (CMC-GEM) atmospheric model  to force the ice-ocean model.

Biogeochemical Model
GSBM simulates the biogeochemical cycles of dissolved oxygen, inorganic carbon and nitrogen, and the biological components that determine the dynamic of the planktonic ecosystem (Figure 2). The model has 11 state variables. It includes both simplified herbivorous and microbial food chains typical of bloom and post-bloom conditions. The export of biogenic matter to depth is mediated by the herbivorous food web (nitrate, large phytoplankton (diatoms), mesozooplankton, detrital particulate organic matter, while the microbial food web (ammonium, small phytoplankton, microzooplankton, dissolved organic matter) is mainly responsible for nutrient recycling in the euphotic zone. Nitrate (NO 3 − ), ammonium (NH 4 ), dissolved organic nitrogen (DON), detrital particulate organic nitrogen (PON), DIC and total alkalinity (TA) are supplied by rivers. The tight coupling between small phytoplankton growth and microzooplankton grazing, microzooplankton egestion and detritus decay, DON remineralization to NH 4 , and nitrification represent the dynamic of the microbial food chain. Biological variables are calculated in nitrogen units and algal biomass and production converted to chlorophyll a (chl a) and carbon units using fixed ratios. Detrital particulate organic matter sinks toward the bottom with a speed that varies regionally as a function of surface salinity to account for the influence of riverine input on organic matter aggregation and sinking. Sinking particles are gradually degraded into DON via a temperature-dependant fragmentation rate. The phytoplankton growth rate is a function of light and nutrient availability. The available light for phytoplankton growth is a function of sea-ice cover, chl a and turbidity (Le Fouest et al., 2006).
Biological processes have an impact on the concentration of chemical constituents in the water column. DO and DIC production and consumption rates, as well as changes in TA are related to nitrogen transformations using constant ratios. Calcium carbonate saturations and pH are calculated using the carbonic acid dissociation constants of Lueker et al. (2000) and pH in the total scale . Phosphate and silicate concentrations are not included in the calculations. The impact of omitting these nutrients in the calculation was evaluated to be smaller than 1% for calcium carbonate saturations in Barrow and Davies Straits (Azetsu-Scott et al., 2010). The main changes made to the original version of LF05 are summarized in Appendix A. A sensitivity analysis is also presented in the Supplementary Material for some of these changes (see Supplementary Figures 1-6). The equations of the model are detailed in section "Biogeochemical Model Equations", while a list of parameters is given in Table 1.

Data Used for Initial and Boundary Conditions, and Validation
Monthly climatologies of nitrate, DO, DIC and TA are used to initialise the model and provide lateral boundary conditions. Monthly climatologies of nitrate, DIC, and TA are also used for river input, except for the St. Lawrence River near Quebec City were DIC and TA are annual means.  (Walsh and Vigneault, 1986). The data were provided in mg N L −1 and converted to mmol N m −3 using a factor of 71.3944 (1000/14.0067).
The frequency of sampling varied among rivers. Sampling was defined as excellent when the 12 months of the year were sampled during at least six years (Figure 3). The data were qualified as good when 8 to 11 months were sampled for at least six years. Below this sampling frequency, the data were qualified as poor. Those rivers represented in the model for which no data are available are represented as No data in Figure 3A. A climatology of monthly nitrate concentration was prepared for the rivers qualified as good or excellent ( Figure 3B). For the rivers where data are insufficient or lacking, we first defined regions FIGURE 3 | Contours of the six regions with rivers having similar characteristics for nitrate (A), dissolved inorganic carbon (C), and total alkalinity (E). The annual mean concentration are given for each rivers (colorbar) when data are available. The shape of the symbols indicates the quality of sampling in terms of frequency or when no data are available (purple circle). The panels to the right show the corresponding monthly climatology for each region and for Quebec City, or the constant annual mean when a climatology could not be calculated (B,D,F). The regions are different for nitrate than for dissolved inorganic carbon and total alkalinity. Frontiers in Marine Science | www.frontiersin.org with similar characteristics (in term of nitrate concentration, physiogeography and land use), and used the mean nitrate climatology of the region to fill in the gaps. The contours of the six regions are presented in Figure 3A.
The rivers flowing into the upper St. Lawrence estuary (region #1) display high nitrate concentrations due to the presence of numerous cities and farms in that area. The GSL's north shore rivers, east of the Saguenay Fjord, as well as the Newfoundland and southern Nova Scotia rivers display the lowest nitrate concentrations (region #2). The northern Nova Scotia (region #5) and Massachusetts rivers (region #6), on the other hand, display the highest nitrate concentrations due to farming activities and large human populations.

Dissolved Inorganic Carbon and Total Alkalinity
The same approach and sources of data described above were used for river DIC and TA, although much less data were available for DIC than for nitrate (Figures 3C-F). However, pH and TA are more often sampled and we used the carbon equilibrium equations mentioned in section "Biogeochemical Model" to derive DIC concentrations. We use an annual mean value of 1350 mmol C m −3 for DIC in the St. Lawrence River near Quebec City, as well as pH values provided by the MELCC to derive an annual mean TA of 1339 mmol C m −3 . These values are slightly higher than those reported in Dinauer and Mucci (2017) for DIC and TA at the same location (1242 mmol C m −3 and 1205 mmol C m −3 respectively).
The regions drawn for the estimation of the DIC and TA climatologies are slightly different than for nitrate (Figures 3C,E). Rivers in regions #1, #2, #4 and #5 present very low DIC and TA concentrations (as low as 50 mmol C m −3 for DIC) while concentrations in rivers from regions #3 and #6 can reach 1000 mmol C m −3 for DIC.
The DIC data were provided in mg C L −1 and converted to mmol C m −3 using a factor of 83.26 (1000/12.0107). The TA data were provided in mg CaCO 3 L −1 and converted to mmol C m −3 using a factor of 9.9913 (1000/100.0869).

Other Variables
Dissolved oxygen concentrations in rivers are calculated from the solubility equation (Benson and Krause, 1984), with river temperature set to the temperature of the nearest ocean grid cell of the model and salinity set to zero. Values for NH 4 , and PON were set to 1.0 and 0.05 mmol N m −3 respectively. DON was set to 3.0 mmol N m −3 in most rivers. These are conservative estimates as the concentrations of these constituents are higher in the rivers flowing in the upper St. Lawrence Estuary (e.g., Hudon et al., 2017) and in other areas with important anthropogenic activities. We kept values in the lower range found for the GSL's north shore river. DON measurements were made at a higher frequency in 8 rivers (St. Lawrence, Manicouagan, aux Outardes, Betsiamites, Mitis, Saguenay, Ouelle and Du Gouffre, Michel Lebeuf, Maurice Lamontagne Institute, personal communication, 2019) and were used to calculate annual mean values that are supplied to the model. These annual mean vary between about 6 mmol N m −3 for 6 of the rivers, 15.8 mmol N m −3 for the Ouelle river, and 19.1 mmol N m −3 for the St. Lawrence River.

Initial and Open Boundary Conditions
An essential element to properly simulate the biogeochemical conditions of the system is to obtain realistic initial and open boundaries conditions, especially for the water masses that flow into the domain (e.g., water properties along the shelf edge south of Newfoundland, and deep water mass of the Laurentian Channel). We gathered as much data as possible over the model domain and just outside for nitrate, DO, DIC and TA to build monthly climatologies. The data were then interpolated over the entire model grid with an optimal linear interpolation technique using the OAX5 program developed at the Bedford Institute of Oceanography, Fisheries and Oceans Canada 4 . The spatial and temporal parameters for the interpolation varied from one variable to the other as a function of data coverage (large for nitrate, average for DO and weak for DIC and TA).

Nitrate
We use nitrate data gathered for the atlas of Brickman and Petrie (2003), that covers the GSL and the Scotian Shelf, data collected during the JGOFS program in the GSL , and data collected during the COUPPB program (Devine et al., 1997) and by Therriault and Levasseur (1985) in the St. Lawrence Estuary. Additional data were obtained from the BioCHEM database , and the World Ocean Database (Garcia et al., 2010b). We used the period 1981 to 2010 to build the nitrate climatology.

Dissolved Oxygen
The DO measurements were obtained from the BioCHEM database , and from the World Ocean Database up to 2006 (Garcia et al., 2010a). The DO concentrations obtained from the World Ocean Database were converted from ml L −1 to mmol m −3 using a factor of 44.6596 (García and Gordon, 1992). As for nitrate, we used data covering the 1981 to 2010 period to build the climatology.

Dissolved Inorganic Carbon and Total Alkalinity
The DIC and TA concentrations used for the interpolation of the initial and boundary conditions were obtained from GLODAPv2 (Key et al., 2015;Olsen et al., 2016), BioCHEM , and from the St. Lawrence Global Observatory (see text footnote 3) for some older data in the LSLE. As for rivers, TA and pH are usually the variables that are measured and we used the carbon equilibrium equations mentioned in section "Biogeochemical Model" to derive DIC concentrations. Since there are fewer observations available compared to nitrate and DO, we kept all the available data to build the monthly climatology (1980 to 2016).

Biogeochemical Model Equations
The model equations for the NPZD model are based on LF05. The choice of the NPZD model parameters is discussed in detail in LF05 and will not be repeated here. Testing of the LF05 parameters and results were also performed in subsequent studies in the GSL (Le Fouest et al., 2006Lavoie et al., 2008;Mei et al., 2010) and in Hudson Bay (Sibert et al., 2011). The addition of a DO module was also tested by Nicot (2013) and Lavoie et al. (2015). However, since changes were made to these components (NPZD and DO) and that a carbonate system module was added, we detail all the equations below. The resulting model was tested over many years, one parameter at a time, to understand the functioning of this complex ecosystem, and whose large domain encompasses freshwater to oceanic conditions. The parameter values that were retained are found in Table 1, while the parameters tested and their ranges found in the literature are presented in Supplementary Table 1.

Large Phytoplankton
The rate of change of the large phytoplankton (LP) biomass is equal to (1) where µ LP is the large phytoplankton specific growth rate, m LP is the large phytoplankton mortality rate, gz MEZ is the mesozooplankton grazing rate (see section "Small phytoplankton"), frac LP is the fraction of the mesozooplankton diet composed of LP, MEZ is the mesozooplankton biomass, and sed LP is the large phytoplankton sinking rate. The phytoplankton growth rate is formulated as in LF05 and is a function of both light and nitrogen availability. However, we modified the formulation of the light-based (dtc) doubling time by using a variable light half-saturation (ke d ) for LP to account for spatial differences in photoacclimatation (e.g., Le Fouest et al., 2010). and where dtmin LP is the minimum doubling time of the biomass, and PAR is the photosynthetically available radiation experienced by phytoplankton. The light half-saturation, ke d , is equal to −26kcdom+16, where kcdom varies with surface salinity (between 0.2492 for salinities smaller than 26 and 0.0302 for salinities higher than 32) following Le Fouest et al. (2006). We use the formulation proposed by Vallina and Le Quéré (2008) for the uptake of dissolved inorganic nitrogen by phytoplankton to avoid the inhibition of ammonium uptake by nitrate, which is present with the O'Neill parameterization used by LF05. For large phytoplankton where NuNO3 LP and NuNH4 LP are the nitrate and ammonium uptake fractions respectively. Values for the half saturation constants for ammonium (k4 LP ) and nitrate (k3 LP ), as well as large phytoplankton mortality (m LP ) and sinking (sed LP ) rates, are found in Table 1. The fraction of the mesozooplankton food composed of large phytoplankton is equal to Mesozooplankton food (food MEZ ) is described in section "Small phytoplankton."

Small Phytoplankton
The rate of change of small phytoplankton (SP) biomass is equal to where µ SP is the small phytoplankton specific growth rate, m SP is the small phytoplankton mortality rate, gz MIZ is the microzooplankton grazing rate (equation 16), frac SP is the fraction of SP used as food by microzooplankton, and MIZ is the mesozooplankton biomass, The equations for small phytoplankton growth and nutrient uptake are similar to those described above for LP except for some of the parameters (see Table 1). The light half-saturation (ke f ) for SP is also constant (9 Ein m −2 d −1 ). SP is not sinking. The fraction of small phytoplankton that is grazed by microzooplankton is equal to Microzooplankton food (food MIZ ) is described in section "Microzooplankton."

Zooplankton Mesozooplankton
The rate of change of mesozooplankton biomass is equal to where ass MEZ is the mesozooplankton food assimilation efficiency (70%), gz MEZ is the mesozooplankton grazing rate, m MEZ the mesozooplankton mortality rate, and ex is the mesozooplankton ammonium excretion rate (0.05 d −1 ).
The mesozooplankton grazing rate is formulated using a sigmoidal "Holling's type III" functional response as in Lavoie et al. (2009).
where gmax MEZ is the maximum mesozooplankton grazing rate (0.5 d −1 ), k MEZ is the half-saturation constant for mesozooplankton grazing (0.9 mmol N m −3 ), and food MEZ is the mesozooplankton food. Mesozooplankton graze on both large phytoplankton and microzooplankton where frac MIZ is the fraction of microzooplankton on which mesozooplankton graze (0.5).
The mesozooplankton mortality, m MEZ , was modified to follow a temperature dependence as in Buitenhuis et al. (2006). This was necessary to improve the representation of the biogeochemical conditions in the western Scotian Shelf and GoM.
where mort MEZ 0 • C is the mortality rate at T = 0 • C (0.048 d −1 ) and b (1.056) is the temperature dependence of mortality (b 10 = 1.72).

Microzooplankton
The rate of change of microzooplankton (MIZ) biomass is equal to where ass MIZ is the fraction of food ingested by microzooplankton that becomes biomass (40%), gz MIZ is the microzooplankton grazing rate, m MIZ the microzooplankton mortality rate (0.02 d −1 ), and the last term is the loss of microzooplankton to mesozooplankton grazing. The microzooplankton grazing is formulated with a sigmoidal "Holling type III" functional response as in LF05, but we added grazing on PON, in addition to grazing on small phytoplankton, and we modified some of the parameters.
where frac PON is the fraction of PON on which microzooplankton graze (0.32) and k MIZ is the half-saturation constant for microzooplankton grazing. The maximum microzooplankton grazing rate (gmax MIZ ) and half-saturation constant (k MIZ ) were slightly reduced compared to LF05 (to 1.9 d −1 and 0.7 mmol N m −3 respectively).

Detrital Particulate Organic Nitrogen
Particular attention has been paid to the fate of PON in this study. To reproduce the main biogeochemical characteristics of the system in both short and long term simulations, the remineralization of sinking PON needs to occur at the right depth. This was achieved by improving the parameterization of the sinking speed and remineralization rate, and by taking into account the influence of freshwater outflow on the sinking rate.
The rate of change of PON is equal to where fg is the fragmentation rate of sinking PON to DON. It is a function of water temperature (T) as in Aumont et al. (2015), where fr PON is the fragmentation rate of PON at 0 • C (0.076 d −1 ) and b is the temperature sensitivity (1.056, b 10 = 1.72). The last term of equation 18 represents the PON changes resulting from sinking in the water column. The sinking rate of PON (sed PON ) varies spatially as a function of surface salinity to account for the influence of river input on the size of PON aggregates (see Lucotte et al., 1991;Al Azhar et al., 2017), These changes results in a greater amount of PON being remineralized in the water column compared to LF05.

Dissolved Organic Matter
Dissolved organic matter is supplied by rivers, and produced through microzooplankton excretion and mortality, small phytoplankton mortality, and through the fragmentation of PON. The rate of change of DON is equal to where eg. (30%) is the fraction of unassimilated food grazed by microzooplankton that is egested as DON, and rem (0.19 d −1 ) is the remineralization rate of DON towards ammonium by bacteria. A fraction of the PON and diatoms (kresed = 0.8) that is deposited to the sediments is returned instantaneously to the pelagic bottom layer (last grid cell in the water column) as ammonium (20%) and DON (80%). This remobilization was added to account for regeneration processes occurring in the sediment which are of importance in shallow areas such as the southern GSL and Scotian Shelf, and at the head of the deep channels.
where z is the thickness of the pelagic bottom layer.

Nutrients Nitrate
Besides the nitrate supply at the model boundaries (rivers and ocean), nitrate production occurs through the oxidation of ammonium, while losses occur through consumption by phytoplankton, where µ LP and NuNO3 LP are the diatom growth rate and NO 3 uptake fraction respectively, µ SP and NuNO3 SP are the flagellate growth rate and NO 3 uptake fraction, and Nitr is the ammonium oxidation rate (nitrification). We implemented nitrification as in Sibert et al. (2011) except that we added two additional terms to take into account the reduced nitrification in the euphotic zone (resulting either from light inhibition or from competition between phytoplankton and nitrifying bacteria) and under low DO conditions (Ward, 2000), which are found in bottom waters of the GSL.
where Nitr max is the maximum nitrification rate, PAR(z) is the light intensity in the middle of layer z, ke Nitr , kb Nitr , and ko Nitr are the half-saturation constants for the dependence of nitrification on light and NH 4 and DO concentrations respectively (see Table 1).

Ammonium
Changes in ammonium concentration result from phytoplankton uptake, excretion by meso-and microzooplankton, remineralization of DON and nitrification.
A small fraction (16%, or 20% of the 80% being returned to the pelagic bottom layer) of the PON and diatoms settling on the bottom is also remineralized to NH 4 ,

Dissolved Oxygen
The DO model follows the approach of Jin et al. (2007) such as described in Nicot (2013). The main change here, compared to the implementation of Nicot (2013), is in the addition of DO consumption by nitrification. The pathways of oxygen production and consumption are depicted in Figure 2.
The changes in DO result from the differences in biological production and consumption, as well as from air-sea fluxes in the surface layer. The changes resulting from biological processes are where the production rate of DO (P_O 2 ) is equal to the sum of new and regenerated production by large and small phytoplankton with different stochiometric O 2 :N molar ratios for new (r NP = 138:16) and regenerated production (r RP = 106:16). Dissolved oxygen is utilized during the conversion of organic nitrogen to ammonium, mainly via remineralization of DON but also through zooplankton metabolism (ammonium excretion by mesozooplankton and microzooplankton), and nitrification. The consumption rate is given by where r Nit , the O 2 :N molar ratio during nitrification, is equal to 2:1 as in Jin et al. (2007). According to Kortzinger et al. (2001) and references therein, the respirative ratio O 2 :C org for the remineralization of organic matter would have a value of about 1.34, which results in an O 2 :N (r RO ) of 8.88 using the redfield ratio for C:N. This is close to our O 2 :N molar ratios for new production defined above (8.625). We thus use the latter ratio (8.625) for the remineralization of organic matter. The air-sea oxygen flux (F O2 ) is parameterized using the standard transfer-velocity formulation of Wanninkhof (1992) for steady winds, based on the difference between the DO concentration at the surface and the oxygen saturation concentration.
where k w is the gas transfer velocity for oxygen where f is the fraction of the sea surface covered with sea ice, Sc O2 is the Schmidt number for oxygen in seawater computed using the formula of Wanninkhof (1992), Sc20 is the Schmidt number in seawater at 20 • C (589.39) and W is the wind speed at 10 m in m s −1 . The oxygen saturation, O 2 sat , is computed from the simulated surface temperature and salinity following the equations established by Benson and Krause (1984).

Dissolved Inorganic Carbon
The solution of the carbonate chemistry equations follows Signorini et al. (2001). The cycling of DIC is coupled to the ecosystem model through biological consumption (C_CO 2 ) and production (P_CO 2 ) rates in a similar fashion to DO, except that it is not affected by nitrification. As for DO, the different pathways are depicted in Figure 2, and the equations are given below, The C:N molar ratios are the same for new and recycled production (rC NP = rC RP = 106:16) and we use the same ratio (r RC = 6.625) for the remineralization of organic matter. The airsea CO 2 flux (FCO 2 ) is regulated by the gradient of the partial pressure ( pCO 2 ) across the air-sea interface, the piston velocity (k w ), and the solubility of CO 2 in seawater (α, in mmol m −3 µatm −1 ) which is a function of seawater temperature and salinity (Weiss, 1974).
The direction of the flux is determined by the gradient of the partial pressure, where a negative flux indicates CO 2 uptake by the ocean, and a positive flux outgassing. The atmospheric pCO 2 (pCO 2 atm ) time series used to force the model is based on the least square fit equation of Signorini et al. (2001) but with two different parameter values (A 0 = 334.05 µatm instead of 292.23, and A 1 = 0.0004 µatm month −1 instead of 0.134) to obtain an increasing linear trend from one year to the other. Their model is at a similar latitude (50 • N) as the northern portion of our model domain. Atmospheric pCO 2 is assumed to be the same over the whole domain. The partial pressure of dissolved CO 2 in seawater is defined by the relationship where [CO 2 ] is the concentration of carbon dioxide in solution. We calculate [CO 2 ] from DIC concentration, TA and the carbonic acid dissociation constants (K 1 and K 2 ) of Lueker et al. (2000). The detailed method for the calculation of [CO 2 ] is provided by Dickson et al. (2007). The piston velocity (k w ) is calculated as in equation 31 for DO but using the coefficients for CO 2 instead of O 2 in the calculation of the Schmidt numbers (Wanninkhof, 1992).

Total Alkalinity
Total alkalinity is unaffected by air-sea exchange of CO 2 , and changed only slightly by biological production and remineralization of organic carbon. The model includes biological sources and sinks of alkalinity associated with nitrogen cycle processes including phytoplankton uptake, remineralization of organic nitrogen, and nitrification. The stoichiometry of these reactions is given by Wolf-Gladrow et al. (2007). Production and dissolution of biogenic carbonate minerals are neglected as is sedimentary denitrification. In the GSL, water column sources and sinks of alkalinity are small relative to fluvial sources (Figure 3), so the main controls on its distribution are dispersion and dilution of river water.

GENERAL MODEL SOLUTION DESCRIPTION AND SKILL ASSESSMENT
The simulation is initiated May 1st 2005 and runs until the end of 2018. It takes a few years for the oceanic conditions near the open boundaries of the model to propagate up to the head of the deep Laurentian Channel (Bugden, 1991;Gilbert, 2004).
We thus compare the model results with observations starting in 2009. Lavoie et al. (2016) demonstrated the performance of the CANOPA model in reproducing the hydrographic conditions of the GSL. Here we make a similar analysis of the model performance for the biogeochemical conditions, i.e., nitrate, ammonium, DIC and DO concentrations, TA, pH, and calcium carbonate saturations (calcite and aragonite). We use data from the same databases that were used to build the model initial and boundary conditions but collected after 2010. We also use ammonium, pH, and calcium carbonate saturation states data obtained from the BioCHEM database and from Pierre Pepin (Northwest Atlantic Fisheries Center, Fisheries and Oceans Canada) to validate the model results. A pointto-point comparison (same time and approximate location) is done between the model output and the observations. We use a Taylor diagram (Jolliff et al., 2009) to visualise the model skill at reproducing the different variables (Figure 4).

Nitrate Concentrations
The spatial distribution of nitrate in the GSL has been mapped on the horizontal (Brickman and Petrie, 2003;Plourde et al., 2014) and on the vertical along the main channels (Coote and Yeats, 1979). The concentration of nitrate is sampled throughout the GSL twice a year, usually in June and late October, and in the surface layer in March, as part of the Atlantic Zone Monitoring Program (AZMP). The model reproduced well the observed horizontal and vertical distribution of nitrate reported in the above papers (correlation coefficient between simulated and observed values of about 0.85; Figure 4), with increasing concentrations towards the deeper layers and from the mouth of the Laurentian Channel towards its heads at depth (Figures 5A,  6A,E). The model also reproduced well the upwelling of nitraterich water at the head of the Laurentian Channel, the so called nutrient pump reported by Steven (1974). The average simulated nitrate flux for year 2010 across the 46-m horizon over the area depicted on Figure 1 at the head of the Laurentian Channel (magenta box) is equal to 18 mmol N m −2 d −1 . Cyr et al. (2015) calculated fluxes between 18 and 300 mmol N m −2 d −1 in the area of maximal upwelling, that correspond to the western end of the box depicted in Figure 1. These nutrient-rich waters are mixed with the nutrient-rich fluvial input that lead to high nitrate concentration in the surface layer. The build up of nitrate in the surface layer and its downstream transport is particularly apparent in winter when nutrient uptake in the mixed layer is reduced and nitrate concentrations in rivers are high (Figures 3A,  5A). A minimum concentration is found in the subsurface layer of the LSLE. This feature is present in both the observations (Greisman and Ingram, 1977) and the model (Figure 6A). Sinclair et al. (1976) and Cyr et al. (2015), contrary to Steven (1974), suggested that the LSLE nutrient pump was playing a minor role in primary production outside the LSLE. However, Sévigny et al. (1979) report that the estimate of Sinclair et al. for the seaward transport of nitrate in the upper layer of the LSLE needs to be corrected upward by a factor of 14 (this would result in a transport varying between 54.6 and 462 mol N s −1 depending on the month) and that the nutrient pump would contribute to primary production in the Gaspé Current. Savenkoff et al. (2001), using inverse modelling, estimated the mean summer seaward nitrate transport from the LSLE in the 0-30 m layer at 460 mol N s −1 , which is similar to the highest (corrected) value of Sinclair et al. (1976). This value corresponds approximately to the mean simulated nitrate transport calculated in the upper 46 m at Pointe-des-Monts between May and September of 2009 to 2016 (472 mol N s −1 ). The simulation shows that May to September is the period of the year when nitrate transport out of the LSLE in the upper layer is minimal (Pointe-des-Monts transect), while in fall and winter, transport over 1500 mol N s −1 can be reached during particularly strong wind forcing events.
Large nitrate transport out of the GSL was also estimated at Cabot Strait in winter. Petrie and Yeats (2000) report that the GSL is the primary source of nitrate for the Scotian Shelf in winter (9300 mol N s −1 ). The simulated nitrate fluxes towards the Scotian Shelf, at the same location and depth as in Petrie and Yeats (2000) are remarkably close to those estimated by the latter in both winter and summer months ( Figure 7E). Overall, the nitrate distribution and transport calculated by the model at various locations in the GSL are in good agreement with those estimated in the literature.

Ammonium Concentrations
Fewer data are available for ammonium concentrations compared to nitrate although they have been routinely collected as part of the AZMP in the last few years. Some measurements have also been reported in the literature for the LSLE (Therriault and Levasseur, 1985;Zakardjian et al., 2000) and GSL (Tremblay et al., 2000). The ammonium concentrations reported in these papers reached up to 4 mmol N m −3 in the LSLE, at a depth of about 50m, while in the GSL, measured values reached 2.5 mmol N m −3 at the base of the euphotic layer, with higher values during summer and on the Magdalen Shallows (up to 7 mmol N m −3 ). The model results are consistent with high ammonium concentration on the Magdalen Shallows and at the base of the euphotic layer following primary production (Figures 5B,D, 6B,F). The model results indicate that part of the organic matter produced in the LSLE and in the nwGSL is transported towards the southern GSL where it is recycled and converted to ammonium and then nitrate. The addition of rapid remineralization of sedimenting PON in the bottom layer is particularly important in the reproduction of high ammonium concentration on the Magdalen Shallows. The recycling of sedimenting material in this shallow region, coupled to the strong stratification, explains the high nutrient concentrations found in the subsurface layer compared to concentrations found at similar depths in the adjacent Laurentian Channel (Coote and Yeats, 1979). It also leads to high simulated ammonium outflow at Cabot Strait ( Figure 7C) in summer, even though the water transport is reduced at that time of the year (e.g., mean velocities in Figure 7B). The Taylor diagram (Figure 4) and the point to point comparison in Figure 6F indicate a low skill at reproducing the observed ammonium concentrations, even though we do reproduced its main features in the system (Figures 5B-D). The exact location of the base of the euphotic layer, were high ammonium values are found, might differ in the model. The timing for organic matter production and advection, and its recycling, can also be slightly different in the model. The analysis of ammonium concentrations was also done on frozen samples in many instances, which could lead to some contamination of the samples. There are thus many reasons that  could explain the reduced correlation between observed and simulated ammonium concentrations.

Timing of the Phytoplankton Spring Bloom
We evaluate the capacity of the model to reproduce the phytoplankton spring bloom by comparing the simulated surface chl a with satellite-derived data averaged over specific areas (Figure 8). We also calculate the timing of the beginning of the spring bloom over each grid cell of the domain (Figure 8H) as the first day when chl a exceeds 1.45 mg m −3 . The map obtained with the simulated data is similar to other maps found in the literature (Siegel et al., 2002;Henson et al., 2009) and reproduces the known progression of the bloom, starting on the Scotian Shelf and then proceeding more or less counterclockwise around the GSL (de Lafontaine et al., 1991) following sea-ice retreat (Figures 8A,C).
The main discrepancy occurs on the Magdalen Shallows, where the simulated bloom is delayed compared to the one estimated from the satellite data ( Figure 8C). For the simulated estimate, the whole area of the box is taken into account, even the areas that are ice-covered (and thus present low chl a values). In the satellite data estimate, measurements are taken only in ice-free water but considered to represent the whole box. These differences could lead to a high bias in chl a concentrations estimated from satellite in early spring. Another discrepancy is noted in the nwGSL, where the fall bloom is not simulated.
In the LSLE, the phytoplankton bloom phenology is controlled in part by the turbidity of freshwater flowing into the estuary, which limits light availability, as well as by the strong horizontal advection, which in conjunction with the diatom sinking rate (larger sinking rates lead to a delayed bloom) can have a great impact on the timing of the bloom (Zakardjian et al., 2000). Satellite chlorophyll algorithms are not reliable in that area due to the large colored dissolved organic matter (CDOM) content of surface waters (Laliberte et al., 2018), and we use in situ data from the AZMP Rimouski station to evaluate the characteristics of the simulated phytoplankton bloom in that region ( Figure 8G). The shape of the simulated chl a curve initially follows the in situ chl a curve, with a peak in June, as reported in the literature (end of June -early July, Levasseur et al., 1984;Therriault and Levasseur, 1985;Starr et al., 2004). However, the model does not capture the bloom characteristics from August to October with much smaller simulated values (falls outside of the observed standard deviation in October). The northeastern part of the LSLE, where inflow from the nwGSL occurs , presents characteristics that are more similar to those of the nwGSL with a bloom in April or May (Therriault and Levasseur, 1985;Roy et al., 2008;Xie et al., 2012).

Primary Production
Few measurements of primary production are available. Nevertheless, the model reproduced the known primary production spatial distribution (Figure 9), with higher primary FIGURE 9 | Averages, over the 2010-2016 period, of (A) primary production (background color, g C m −2 y −1 ), (B) percentage of primary production that reaches the sediments, (C) zooplankton production (g C m −2 y −1 ), and (D) particulate organic matter flux across the 150-m horizon. The average primary production calculated by the model as well as its standard deviation is displayed in italic in panel (A) for different regions (LSLE, NWGSL, AC, EC, CGSL, CB, MS, CS, SS, and GoM). The range of estimates from the literature is also displayed below the simulated average when available. The estimates are taken from Sinclair et al. (1976); Levasseur et al. (1984), and Therriault and Levasseur (1985) in LSLE, Sévigny et al. (1979) and Savenkoff et al. (1996) in the NWGSL, in Savenkoff et al. (1996) for AC+EC and CS, in Legendre (1971) for CB, in O' Reilly and Busch (1984) for the GoM, in Fournier et al. (1977); Mills and Fournier (1979), Mousseau et al. (1996), and Shadwick et al. (2011) for the SS. The 200 m isobaths is displayed in (B). Station S1 of Savenkoff et al. (1996) is displayed in panel (D). production in the LSLE compared to the rest of the GSL (Steven, 1974;de Lafontaine et al., 1991;Fuentes-Yaco et al., 1997), as well as larger primary production on the Scotian Shelf and GoM compared to the GSL (Table 7 of Therriault and Levasseur, 1985). The mean primary production simulated by the model for years 2010 to 2016 falls within the values reported in the literature (Figure 9A). The main discrepancies between the simulated primary production and those reported in the literature occur in the Gaspé Current region, where simulated production is too low, and in the northeast GSL where it seems higher than observed [see Figure 6 of Lavoie et al. (2008)]. The simulated secondary production is high in the Gaspé Current region (Figure 9C), which most likely leads to the low primary production simulated by the model in this area.

Particulate Organic Matter
To successfully reproduce nutrient and DO concentrations as well as pH and calcium carbonate saturation states in the intermediate and deep layers, sinking and fragmentation rates of particulate organic matter must be represented correctly. The model results are very sensitive to these parameters (see Supplementary Figures 1, 2 in the Supplementary Material).
Particulate and dissolved organic matter data are not available for a direct comparison but we can compare with concentrations and fluxes reported in the literature. Care must be taken when doing this comparison however, as only PON and DON of marine origin, as well as the labile fraction of riverine DON are considered in the model, while observations also include terrigeneous sources. Lucotte et al. (1991) estimated that less than 30% of the primary production in the LSLE reaches the bottom sediment (mostly as fecal pellets), most of it (70%) being recycled in the water column. In the GSL, a lower percentage of the primary production would reach the sediment (7-10%; Savenkoff et al., 1996;Silverberg et al., 2000). These results are in line with a faster PON sinking speed in the LSLE (see section "Detrital Particulate Organic Nitrogen"). The model results are in accordance with these estimates, with 5 to 20% of the primary production reaching the sediments in the deep channels of the GSL (Figure 9B), the highest values being found in the downstream part of the LSLE. In shallower regions, this percentage is higher and reaches up to 50% around the Gaspé Peninsula ( Figure 9B). This high percentage results in large part from the advection of phytoplankton and zooplankton in this region that leads to high PON sedimentation.
The particle flux calculated at 150 m in the model ( Figure 9D) does not correspond as well with those reported by Savenkoff et al. (1996) at their station S1 (see Figure 9D), with estimated and simulated values of 45 ± 12 mg C m −2 d −1 and 80 mg C m −2 d −1 respectively. However, it is to note that their estimate is based on multiple trap sampling that lasted a few days only over each of their nine cruises (between April and December), while the model value represents continuous sampling over a similar period. In Savenkoff et al. (2000), the mean particulate organic carbon flux (mean of the 5 JGOFS stations) at 50 m was 119 ± 28 mg C m −2 d −1 in winter-spring (November-December and April) and 161 ± 65 mg C m −2 d −1 in summer-fall (June-October).

Dissolved Oxygen Concentrations
The vertical and horizontal distributions of DO concentrations (Figures 6C,G, 10), as well as the bottom DO saturation (Figure 11) are successfully reproduced by the model. The Taylor diagram (Figure 4) shows a correlation coefficient of 0.85 for DO over the whole domain. This correlation climbs to 0.92 when comparing data from the GSL only. In the LSLE, PON sinking speed is enhanced by flocculation processes and thus more organic matter is remineralized at depth. The greater DO consumption combined with the reduced ventilation of the water below 175 m (Lavoie et al., 2000;Belzile et al., 2016) leads to low DO concentrations (50-100 mmol m −3 ) and saturation values (as low as 20%) in the deeper layers of the LSLE (Figures 6C, 11). The deeper parts of the GSL (below 325 m in the eastern part of the Laurentian Channel) contain water entering from the North Atlantic which are characterized by higher DO concentrations (160-180 mmol m −3 , Figure 6C). Figure 6C shows that upwelling of the deeper layers occurring at the head of the Laurentian Channel also brings water with lower DO concentration near the surface.
Besides the LSLE hypoxic zone, the model reproduced many other observed features of the DO distribution, such as the horse shoe shape of the 35% saturation contour in the nwGSL (although our annual mean does not extend as far along the southern edge of the Laurentian Channel, Figure 11), and the minimum DO along the edge of the Scotian Shelf (see Figure 2 of Bianucci et al., 2016) and in some basins on the Scotian Shelf (Figure 11). Observed and simulated DO profiles in the GoM are also in good agreement (Figures 10B,C). The simulated DO seems to be a little high in July ( Figure 10B) compared to observed profiles, although it remains within the values observed over a few years (shaded area).
Low DO saturation (below 35%) is also observed in the Esquiman and Anticosti Channels (white dots in Figure 11). Simulated DO saturation values are low in these channels but do not reach 35% (they rather oscillate around 40-45%). This is explained in part by a shallower depth in the model compared to reality. The comparison of bottom values is thus not made at the same depth. The amount of PON reaching the bottom might also be insufficient, due to a small amount of mesozooplankton, or to non-optimal PON sinking speed and remineralisation rates. The inclusion of a sediment box with slow release from the sediments rather than instant remineralisation could also improve DO concentrations and saturation levels in general, especially in regions where the productive period is short.
Our mean respiration rates in the water column (50-200 nmol O 2 L −1 per day, Supplementary Figure 7) are similar to those estimated by Bourgault et al. (2012) (30-120 nmol O 2 L −1 per day), although their vertical distribution differs. Bourgault et al. (2012) place minimum respiration rates in the cold intermediate layer, in accordance with the temperature dependency used for their pelagic respiration rate. Although we do include a temperature dependence of PON fragmentation rates, our minimum values are found in the deep layers in the eastern part of the GSL that receives a smaller amount of organic matter compared to other regions. However, high respiration rates are simulated in the LSLE and nwGSL near the bottom (Pointe-des-Monts area in Supplementary Figure 7) due to the greater accumulation of PON near the bottom in that region (Figures 9B,D). These respiration rates are, however, low compared to other estuarine systems (Fennel and Testa, 2019).

Carbonate System
In this section we evaluate the simulated air-sea CO 2 fluxes as well as pH and calcium carbonate saturation states at the bottom water that are particularly low in the LSLE (Mucci et al., 2011). Simulated air-sea CO 2 fluxes are compared with data reported in the literature. In the LSLE and GSL, air-sea CO 2 fluxes have been estimated by Dinauer and Mucci (2017). The LSLE is a source of atmospheric CO 2 in general due to the strong influence of freshwater inflow that carries organic matter into the LSLE where it is remineralized and leads to an increase in surface water pCO 2. . This process is counteracted during late spring and summer by the strong primary production taking place in the LSLE that consumes DIC and leads to low surface pCO 2 (Supplementary Figure 8) and a reversal of the flux. The range of simulated air-sea CO 2 flux calculated for the months of May  to July (−3.1 ± 5.2 mmol C m −2 d −1 , with values ranging between −43.8 to 26 mmol C m −2 d −1 over years 2010 to 2016) corresponds to the estimates of Dinauer and Mucci (2017) for the same months and years (0.8 ± 7.2 mmol C m −2 d −1 , range −21.9 to 15.1 mmol C m −2 d −1 ). The means have opposite sign but are small in both cases. A negative sign represents CO 2 being taken up by the ocean, while a positive sign represents outgassing of CO 2 to the atmosphere. The difference in the mean can result from our earlier bloom timing in the LSLE (Figure 8G), to the fact that the model covers the entire area of the LSLE, as well as the entire three-months period each year, which is not the case with the field sampling.
In contrast to the LSLE, the GSL is a sink of atmospheric CO 2 in winter and during the spring bloom (April-May). Outgassing occurs afterwards as the sea surface warms (see flux reversal at the end of May in Supplementary Figure 8). On average, Dinauer and Mucci (2017) estimated that the GSL would be a weak sink (−1.2 ± 3.0 mmol C m −2 d −1 , with values spatially ranging between −8.4 to 3.6 mmol C m −2 d −1 ). Here again, the model results show a similarly low mean (1.35 ± 6.5 mmol C m −2 d −1 ) but with a large range of fluxes (−43.8 to 17.8 mmol C m −2 d −1 ). The CO 2 fluxes calculated in the boxes NWGSL, CGSL, AS, EC and CA represented in Figure 9A were used to calculate the mean (rather than the GSL box of Figure 1) for a better comparison with the estimations of Dinauer and Mucci (2017) that do not cover the southern GSL.
The model also reproduces well the low pH in the Laurentian Channel (Figures 6D,H, 12A), and at the Shediac Station (7.7 to 7.9, Menu-Courey et al. (2019), see Figure 1 for the location of the Shediac Station). The model also reproduces the main features of bottom water aragonite saturation state. Most of the bottom waters of the GSL are undersaturated with respect to aragonite. This is also true for the shallower southern GSL where important remineralization of organic matter occurs, as mentioned previously. We do not have data for a direct comparison of aragonite saturation in the GoM but Wang et al. (2013) report aragonite saturation states below 1.3 in the deep part of the GoM, which corresponds to our simulated values ( Figure 12B). The simulated values on the Scotian Shelf are, however, slightly low ( Figure 12B).

MODEL APPLICATIONS
The performance of the biogeochemical model was demonstrated above. We now use the model results to look at the impact of the circulation modes described in Lavoie et al. (2016) on the concentration and distribution of the biogeochemical variables. We also use the model to estimate the contribution of rivers in the total supply of nitrate to the surface layer of the LSLE. Lavoie et al. (2016) described two mean seasonal circulation modes (winter-spring and summer-fall), different than the ones usually described as a function of surface stratification (springsummer and fall-winter), and driven mainly by local and large-scale wind forcing. In winter and spring, a "semi-closed" circulation mode is present, with many strong gyres and reduced exchanges between the western and eastern GSL. The inflow at the Strait of Belle Isle and at Cabot Strait is stronger in winter, which, in the simulation, ultimately leads to an increased mean water level in the LSLE (see Figure 11 of Lavoie et al., 2016). The increased water level creates a west to east gradient that combines with the strong northwesterly winds to oppose the stronger inflow at the GSL's bounding straits.

Impact of Circulation Modes on the Seasonality of Biogeochemical Variables
The nitrate flux from the St. Lawrence River is maximal in winter and spring, due to the combination of high runoff and high nitrate concentration ( Figure 3B). The depth from which water is upwelled at the head of the Laurentian Channel is also maximal in winter (Saucier et al., 2009), which combined with a shoaling of the upper limit of Atlantic water (Galbraith, 2006), leads to upwelling of water with higher salinity and temperature (Belzile et al., 2016), higher nitrate concentration, and lower DO and pH than in summer. These two processes lead to the delivery of a large amount of nutrients to the LSLE when primary production is at a minimum. Because of the "semi-closed" configuration of the system at this time of the year, with its many gyres, the nutrients accumulate mostly in the LSLE and nwGSL. Some of this nutrient-enriched surface water also makes its way around the Gaspé Peninsula ( Figure 5A). The downstream advection of nitrate directly from the LSLE or nwGSL towards the Scotian Shelf is comparatively low, but combined with the recycling of organic matter (see section "Ammonium Concentrations"), nitrate concentrations can reach about 6-8 mmol N m −3 in the southern GSL in winter. These concentrations combined with strong transport at this time of the year lead to important nitrate transport out of the GSL at Cabot Strait, both observed and simulated (Figures 7A,E). In spring, increased river runoff and surface heat flux increase stratification, light availability increases and surface nutrients are used up by primary producers in most regions, except in the LSLE due to greater light limitation caused by turbid waters.
During summer, the transport at the GSL's bounding straits is minimal, the northwest winds are weak and the westward coastal current along the northern coast of the GSL is stronger and often flows as far as the Saguenay fjord in the LSLE . The Anticosti gyre is not closed (more of a U shape), and the outflow of the St. Lawrence River remains close to the south shore of the LSLE. The circulation in the eastern GSL is weaker. Overall, there are fewer gyres present in the system in the second circulation mode. Primary production is still high in the LSLE and Gaspé Current during summer, due to the constant supply of nutrient by upwelling at the Laurentian Channel head. During this more "open" circulation mode, less regional retention is occurring and a larger part of the phytoplankton, zooplankton, and particulate and dissolved organic material produced in the LSLE and Gaspé Current are transported towards the Magdalen Shallow where they are recycled (Plourde and Runge, 1993;Zakardjian et al., 2000;Savenkoff et al., 2001; Figure 15 of Ouellet et al., 2013). This recycling leads to the high ammonium concentrations reported in section "Ammonium Concentrations" (up to 7 mmol N m −3 depending on the year, Figures 5B,D), and low DO, pH and aragonite saturation compared to other relatively shallow regions (Figures 11, 12). The high ammonium concentration observed in late summer and early fall are then oxidized to nitrate, which will then be available for transport towards the Scotian Shelf. Outflow of inorganic nutrient and organic material also occurs along the southern edge of the Laurentian Channel and eventually flows into the deep channels and basins of the Scotian Shelf. Remineralization of this organic matter along the way leads to a nitrate maximum and a DO minimum along the shelf edge, from within the GSL down to the Gulf of Maine (Figure 11).
In summary, nutrients are made available through upwelling and river flux to the surface of the LSLE, nwGSL and partly to the southern GSL. As these nutrients are used up, the organic matter produced is transported towards the southern GSL, where it is recycled to nitrate that is transported towards the Scotian Shelf where it will be used up in spring. The combination of the dynamics of the planktonic system with the circulation modes over the year leads to a cascade that eventually bring the nutrients from far upstream, down to the Scotian Shelf and most likely as far as the GoM, through one or two transformation cycles.

Nitrate Sources and Relative Contribution of Rivers
The demonstrated skill of the coupled biogeochemical model and the inclusion of riverine input allows us to revisit important questions about nitrate sources in the system and their relative contribution. The main sources of nutrients to the GSL are the St. Lawrence River (which combines high runoff with high nutrient concentrations), inflow at the Strait of Belle Isle and Cabot Strait, and atmospheric deposition. The importance of each of these contributions to the GSL nutrient budget is investigated in a separate paper in preparation.
Here we compare the riverine input with the upwelling flux at the Laurentian Channel head. The latter is estimated as the mean flux through the base of the 6th layer, located at 46 m, averaged over the purple box displayed in Figure 1.
The St. Lawrence River represents a sizable nutrient contribution to the surface layers of the LSLE, depending on the time of the year. Using monthly mean concentrations and runoff at Quebec City (corresponding to the model input at this boundary) we obtain a nitrate flux of 288 mol N s −1 (12.7 × 10 4 tons y −1 ). This is twice the flux of 5.7 × 10 4 tons y −1 estimated by Greisman and Ingram (1977) using a mean nitrate concentrations of 13 mmol N m −3 in Lac St. Pierre and a mean runoff of 10,000 m 3 s −1 . Their estimate is based on data sampled in June-July 1975. Our estimate is also higher than the estimate of Coote and Yeats (1979) of 8.4 × 10 4 tons y −1 who used annual means for nitrate concentration (14 mmol N m −3 ) and runoff (13,000 m 3 s −1 ) at Quebec City. Our estimate is higher because we use monthly values, where the maximum runoff corresponds to the maximum nitrate concentration (30.6 mmol N m −3 in April) and also because our annual mean nitrate concentration, based on more recent data, is higher (22.5 mmol N m −3 ). Our mean runoff is similar at 11,000 m 3 s −1 . The total nitrate contribution of all the rivers included in the model and that flow into the GSL is equal to 324 mol N s −1 (14.3 × 10 4 tons y −1 ). The St. Lawrence River is clearly the main riverine contributor, supplying nearly 90% of nitrate from rivers. Greisman and Ingram (1977) estimated that the St. Lawrence River contribution to the nitrate pool of the LSLE surface layer was less than 25% based on their June-July measurements. Using inverse modelling, Savenkoff et al. (2001) obtained a similar contribution, 18%, for the July to September period. However, Gilbert et al. (2007) estimated that in spring, when both nutrient concentration and runoff are at their maximum, this relative contribution could reach 50%. The other half would be supplied by the intense vertical mixing and upwelling at the Laurentian Channel head.
The model results are mostly in accordance with these estimates with a mean St. Lawrence River contribution of about 25% in summer (June to September) and a higher river contribution from late fall to early spring of about 36%. The minimum St. Lawrence River contributions were simulated in July 2012, August 2009, 2012, and September 2014. The largest St. Lawrence River contribution was simulated in December 2010 (51%) and in March and April 2017 (42 and 45% respectively). Our mean winter-early spring riverine contribution is smaller than the 50% estimated by Gilbert et al. (2007) because even though the riverine flux increases, water that upwells at the head of the Laurentian Channel during this period also contains a greater amount of nutrients (see section "Impact of Circulation Modes on the Seasonality of Biogeochemical Variables"). The very large contribution in December 2010 (51%) is anomalous (the December average is 36%) and results from strong easterly winds that led to high water levels in the LSLE and reduced upwelling.

CONCLUSION
A new biogeochemical model of the GSL, based on LF05, but to which numerous additions/changes were made, including the addition of a dissolved oxygen and dissolved inorganic carbon components, was presented. Changes to LF05 were made in order to maintain the biogeochemical characteristics over multiyear runs. This is essential if we want to use the model to study the inter-annual variability and make climate projections. These important changes included the addition of nitrification, modification to the particulate organic matter sinking and fragmentation rates, and instantaneous remineralisation of part of the organic matter that reaches the sediments. The impact of these changes is discussed in more detail in the Supplementary Material. The model skill was demonstrated by a comparison with numerous data as well as with fluxes and budgets reported in the literature. The stability of the model is demonstrated by the good comparison with observations (up to 2017), many years after the model initialisation (2005). Annual and inter-annual variability are also specifically represented in Figures 7, 10.
The results highlight the importance of the remineralization of organic matter, partly imported from upstream regions, taking place in the southern Gulf. This process leads to high ammonium concentrations and low DO, pH and calcium carbonate saturation states for such a shallow region.
The processes taking place in this region also contribute to a transformation loop that contributes to the transport of nutrients from the GSL to the Scotian Shelf and eventually to the GoM. Thus, the GSL does act as a nutrient trap as stated by Coote and Yeats (1979), with a nutrient increase as the water moves up towards the head of the Laurentian Channel, but it also represents a nutrient source to downstream regions. The transformation of ammonium into nitrate also maintains high primary productivity on the Magdalen Shallows, which support important fisheries. This process was not reproduced in LF05.
With the inclusion of many rivers, and their respective nitrate inputs, the model allowed us to calculate the relative contribution of rivers to the nitrate content of the mixed layer in the LSLE, compared to the upwelling at the head of the Laurentian Channel. Our estimates were similar to those estimated in the literature for the summer period. However, we showed that the river contribution was larger in winter and early spring, but not as much as had been previously estimated.
The model represents a tool that can be used to elucidate the complex functioning of the GSL ecosystem and answer different questions to support fisheries and ecosystem management. Two additional studies using this numerical model are under way to analyse the seasonal variability of the carbonate system, and the processes governing the interannual variability of primary production in the GSL. Next steps for model development include the addition of the effect of temperature on all biological processes to study the impact of climate change on the biogeochemical state of the system, the inclusion of a simple sediment model to represent the temporal delay in remineralization of sedimenting PON, and forcing with annually variable boundary conditions.

DATA AVAILABILITY STATEMENT
The observations used in the paper are available through public repositories. The source code, boundary conditions and other databases as well as simulation raw data supporting the conclusions of this paper are very large and thus difficult to drop in a public repository. However, all are archived on the Fisheries and Oceans Canada servers and available on request without undue reservation.

AUTHOR CONTRIBUTIONS
DL, MS, and JC conceived the project. JC helped with the coupling between NEMO and GSBM. YL implemented the oxygen module. JRC provided the carbonate system code that BB implemented into GSBM. MS and KA-S provided data and advice with the carbonate system. DG provided dissolved oxygen observational data. NL implement total alkalinity as a tracer. DL, NL, JC, and OR were responsible for the coupled model improvements. DL, NL and OR were responsible for the simulations and analysis. DL wrote the first draft of the manuscript. All authors provided feedback and contributed to the final manuscript.