Biogeosciences Present nitrogen and carbon dynamics in the Scheldt estuary using a novel 1-D model

A 1-D, pelagic, reactive-transport model of a completely mixed, turbid, heterotrophic estuary – the Scheldt estuary – is presented. The model resolves major carbon and nitrogen species and oxygen, as well as pH. The model features two organic matter degradation pathways, oxic mineralisation and denitrification, and includes primary production as well as nitrification. Apart from advectivedispersive transport along the length axis, the model also describes O2, CO2, and N2 air-water exchange. The aim of this study is to present a model which is as simple as possible but still fits the data well enough to determine the fate and turnover of nutrients entering the estuary and their spatial patterns in the years 2000 to 2004. Nitrification is identified as one of the most important processes in the estuary, consuming a comparable amount of oxygen as oxic mineralisation (1.7 Gmol O2 y−1 vs. 2.7 Gmol O2 y−1). About 10% of the 2.5 Gmol of nitrogen entering the estuary per year is lost within the estuary due to denitrification. Nitrogen and carbon budgets are compared to budgets from the seventies and eighties, showing that nitrification activity has peaked in the eighties, while denitrification steadily declined. Our model estimates an average CO 2 emission of 3.3 Gmol y −1 in the years 2001 to 2004, which is a comparatively low estimate in the context of previous estimates of CO 2 export from the Scheldt estuary. Correspondence to: A. F. Hofmann (a.hofmann@nioo.knaw.nl)


Introduction
Estuaries play an important role in the transfer of land derived nutrients and organic matter to the coastal ocean as they act as bio-reactors where both the quantity and the quality of the constituents are altered (Soetaert et al., 2006).Additionally, Borges et al. (2006) and Frankignoulle et al. (1998) report that estuaries are globally important sources of CO 2 , by ventilating riverine dissolved inorganic carbon (DIC) as well as DIC originating from the degradation of riverine organic matter.This estuarine filter function is not a static property and evolves with changing forcings (Cloern, 2001).Thorough assessment of evolving estuarine filter functioning requires not only access to long-term datasets documenting changes in concentrations and loadings, but also biogeochemical models that allow reproducing these trends and deriving transformation intensities.
Long term changes in estuarine biogeochemical cycling have been well documented for the Scheldt estuary (Soetaert et al., 2006).Billen et al. (1985) predicted that increased oxygen levels due to lowered organic loadings would increase delivery of nitrogen to the North Sea.Soetaert and Herman (1995a) indeed found that an increased percentage of riverine nitrogen was exported to the North Sea in the eighties as compared to the seventies.Recently, Soetaert et al. (2006) report major changes for the Scheldt estuary, not only in nutrient loadings, but also changes in nitrogen and phosphorus retention and regeneration, from the mid sixties to the beginning of the 21st century.Against this background it is of interest to investigate the interplay of different biogeochemical processes as well as to quantify import from and export to both the North Sea and the atmosphere in the Scheldt estuary during the years 2000 to 2004.
For macrotidal estuaries such as the Scheldt estuary, Vanderborght et al. (2007) identified reactive-transport modelling as a powerful approach to investigate nutrients and Published by Copernicus Publications on behalf of the European Geosciences Union.carbon budgets and cycling at the estuarine scale.To answer questions about seasonal dynamics, a fully coupled, two-dimensional hydrodynamic and reactive-transport model (e.g.Vanderborght et al., 2007) is the method of choice.However, for the assessment of an average situation for a certain period of time that can be compared to other decades, a simple 1-D model reproducing yearly averaged values is sufficient (cf.also 1-D approaches by Soetaert andHerman, 1994, 1995a,b,c;Regnier et al., 1997).
The objective of this study was to construct a 1-D, biogeochemical, pelagic, reactive-transport model of the mixed, turbid, heterotrophic Scheldt estuary.The model is deliberately kept as simple as possible (such that the performance is still sufficient) to optimize the ratio between the number of parameters and data available for calibration and validation.Furthermore, the level of complexity of the physical transport processes is kept comparable to the level of complexity of the biogeochemical processes.We use this model to obtain whole estuarine budgets for dissolved inorganic carbon, oxygen and nitrate, and ammonium, determining the fate and turnover of nutrients entering the estuary and describing the spatial patterns of nutrient concentrations and fluxes in the Scheldt estuary in the years 2000 to 2004.We compare our nitrogen budget to budgets from the seventies (Billen et al., 1985) and eighties (Soetaert and Herman, 1995a).The obtained value for an average CO 2 export to the atmosphere in the years 2001 to 2004 will be compared to reported CO 2 water-air flux estimates for the Scheldt estuary (Frankignoulle et al., 1998;Hellings et al., 2001;Vanderborght et al., 2002;Gazeau et al., 2005).
Recently, Vanderborght et al. (2002) and Vanderborght et al. (2007) also published modelling studies for the Scheldt estuary.Both of those models include the physical transport of solutes in the estuary in great detail (Vanderborght et al. (2002) describes a tidally resolved 1-D model and Vanderborght et al. (2007) a tidally resolved 2-D model).This rather complex description of physical processes leads to severe limitations of the models: 1. there is a high demand for high resolution data which makes it difficult to run the model for past decades with scarce data coverage or for predictive future scenarios, 2. the computational demand of the models is rather high, rendering it difficult to run the model for longer model times, 3. and to port these models to other systems, detailed bathymetrical maps are needed which might not always be available.
Biogeosciences, 5, 981-1006, 2008 www.biogeosciences.net/5/981/2008/4. Furthermore, the complex representation of physics together with a rather crude representation of biogeochemistry in the models of Vanderborght et al. (2002) and Vanderborght et al. (2007) implicitly puts emphasis on the role of physics for the estuarine ecosystem functioning, which might not correspond to reality.
Since our model suffers from none of these drawbacks, our model code is public domain (the model codes of Vanderborght et al. (2002) and Vanderborght et al. (2007) cannot be publicly verified), and our model uses more recent data, we feel that this paper provides a complementary contribution to the scientific literature.

The Scheldt estuary
The Scheldt estuary is situated in the southwest Netherlands and northern Belgium (Fig. 1).
The roughly 350 km (Soetaert et al., 2006;Van Damme et al., 2005) long Scheldt river drains a basin of around 21 500 km 2 (average from numbers given in Soetaert et al., 2006;Vanderborght et al., 2007;Van Damme et al., 2005 andMeire et al., 2005) located in the northwest of France, the west of Belgium and the southwest of the Netherlands (Soetaert et al., 2006).The hydrographical basin of the Scheldt contains one of the most densely populated areas in Europe (Vanderborght et al., 2007) with about 10 million inhabitants (Meire et al., 2005;Soetaert et al., 2006), resulting in an average of 465 inhabitants per km 2 .Anthropogenic eutrophication and pollution of the Scheldt estuary therefore are of considerable magnitude, especially due to the poor waste water treatment in upstream areas, e.g. in Brussels (Meire et al., 2005;Van Damme et al., 2005;Soetaert et al., 2006;Vanderborght et al., 2007).The water movement in the Scheldt estuary is dominated by huge tidal displacements with around 200 times more water entering the estuary during a flood than the freshwater discharge during one tidal cycle ( Vanderborght et al., 2007).The average freshwater flow is around 100 m 3 s −1 (Heip, 1988).The cross sectional area of the estuarine channel shows a quite regular trumpet-like shape opening up from around 4000 m 2 upstream to around 75 000 m 2 downstream (Fig. 2; Soetaert et al., 2006) whilst the mean water depth varies quite irregularly between values of 6 m and 14 m with the deepest areas towards the downstream boundary (Soetaert and Herman, 1995b).The estuary has a total tidally averaged volume of about 3.619 10 9 m 3 and a total tidally averaged surface area of 338 km 2 (Soetaert et al., 2006;Soetaert and Herman, 1995b), the major parts of which are situated in the downstream area.Peters and Sterling (1976), as cited by Vanderborght et al. (2007), divide the Scheldt river into three zones: The first zone, between the estuarine mouth at Vlissingen and Walsoorden, consists of a complex system of flood and ebb channels and a moderate longitudinal salinity gradient within the polyhaline range.The second zone from Walsoorden to Rupelmonde is characterised by a well defined river channel and a steep salinity gradient with salinities at Rupelmonde between 0 and 5 (Meire et al., 2005).The third zone, upstream from Rupelmonde, consists of the Scheldt freshwater river system and various tributaries, with tidal influence in the Scheldt up to the sluices of Gent (Van Damme et al., 2005).The model presented here comprises the first and the second zone with an upstream boundary at Rupelmonde (river km 0) and a downstream boundary at Vlissingen (river km 104).

Biogeochemical processes
The model describes four main biogeochemical processes, oxic mineralisation (R Ox ), denitrification (R Den ), nitrification (R Nit ), and primary production (R PP ) as given in Table 1.Although primary production is limited in the Scheldt estuary (much lower than respiration in Gazeau et al. (2005) and one order of magnitude lower than respiration in Vanderborght et al., 2002), it is included in the model since it can be of significance in the downstream stretches of our model domain (Soetaert and Herman, 1995c).The primary production description considered in this model lumps the carbon fixing activity of several groups of benthic and pelagic organisms.Oxic mineralisation and denitrification are included as the two pathways of organic matter degradation that are energetically most favourable (Canfield et al., 2005).Nitrification is included in the model since it is one of the most important O 2 consuming processes in the Scheldt estuary (Soetaert and Herman, 1995a;Andersson et al., 2006).
Denitrification is a heterotrophic process and consumes organic matter.The nitrogen of the organic matter mineralised by denitrification can be assumed to be directly oxidised with nitrate to form N 2 (Froelich et al., 1979), which represents an implicit coupling of denitrification and the anammox (anaerobic ammonium oxidation) process.However, denitrification can also be modelled as a source of ammonium if all the organic nitrogen is assumed to be released as ammonium (Froelich et al., 1979;Joergensen, 1983).To keep the model as simple as possible and because of the absence of data on anammox, we model denitrification as a source of ammonium.
The incorporation of sulfate reduction, sulfide oxidation, biogenic calcification, as well as calcite and aragonite precipitation has been tested, but model performances did not increase significantly.In order to meet the objective of constructing a model as simple as possible we opted against inclusion of these processes in the model.Reduction of manganese and iron oxyhydroxides is neglected because of their low concentrations.Similarly, due to methane concentrations on the nano molar scale (Middelburg et al., 2002), methanogenesis is also neglected.
Modelling net respiration (respiration -primary production) yields approximately the same model performance as modelling respiration and primary production individually.Following reviewer feedback, we opted for the latter due to the increased amount of information that can be extracted from the model in that way.
Table 2 shows kinetic formulations for all processes considered in the model.As in Soetaert and Herman (1995a), organic matter has been split up into two different fractions, one fast decaying fraction FastOM with a low C/Nratio of γ FastOM and one slow decaying fraction SlowOM with a high C/N-ratio of γ SlowOM , respectively.R Ox and R Den are modelled as first order processes with respect to organic matter concentration in terms of organic nitrogen as done by Soetaert and Herman (1995a).As in Regnier et al. (1997) and Soetaert and Herman (1995a), oxygen concentration inhibits denitrification.Nitrate concentration influences denitrification according to a Monod relationship as given by Regnier et al. (1997).The molecular nitrogen produced by denitrification is assumed to be immediately lost to the atmosphere.Mineralisation processes are thus modelled in the same order of sequence as in diagenetic models, depending on the availability of different oxidants with different free energy yields per mole organic matter oxidised (e.g.Froelich et al., 1979;Soetaert et al., 1996;Canfield et al., 2005).
Total mineralisation rates for fast and slow degrading organic matter are assumed to be only dependent on temperature and not on the mineralisation pathway.To ensure that the total temperature dependent mineralisation rates rt Min X , for FastOM and SlowOM respectively, are always realised, the independently calculated limitation factors for the two potentially simultaneously occurring mineralisation processes are rescaled as in Soetaert et al. (1996) to partition the total mineralisation rates.
Note that denitrification is an anaerobic process (Canfield et al., 2005) and does not occur in an oxic water column but in the anoxic layers of the sediment.To keep the model as simple as possible, we decided to not include sedimentwater exchange of solutes.As a result of that, we model denitrification with parameters allowing a certain extent of denitrification in an oxic watercolumn as a proxy for benthic denitrification (following Regnier et al., 1997).Therefore, the parameters k Inh O 2 and k NO − 3 should not be considered true microbiological rate parameters.Since we did not include sediment-water exchange in our model, also other pelagically modelled processes incorporate in part sedimentary processes.
Nitrification is expressed as a first order process with respect to ammonium concentration and with a Monod dependency on oxygen as in Soetaert et al. (1996).The same oxygen half saturation parameter k O 2 as for oxic mineralisation is used for nitrification.The nitrification may depend on the microbial community of nitrifiers.de Bie et al. (2001) documented dramatic shifts in nitrifier populations along the estuary from the freshwater to the marine part.Therefore, we assume that the nitrification activity in the freshwater (upstream) part of the modelled system is performed by freshwater nitrifying organisms that experience increasing stress as salinity increases.As a consequence, their activity collapses in the downstream region due to an insufficient adaptation to marine conditions (see also Helder and de Vries, 1983) which is incompletely compensated by activity of marine nitrifying species, so that the overall nitrification activity in the estuary decreases with increasing salinity.This is expressed with an inverse dependency of nitrification on salinity, according to a sigmoid function based on a Holling Type III functional response (Gurney and Nisbet, 1998), modified to have a value of 1 for zero salinity and to asymptotically reach an offset value of o Nit for high salinities.Biogeosciences, 5, 981-1006Biogeosciences, 5, 981- , 2008 www.biogeosciences.net/5/981/2008/Table 2. Kinetic process formulation with two organic matter fractions (X ∈ {FastOM, SlowOM }).D signifies the mean water depth of the estuary, and T urb its relative turbidity (T urb=0 downstream; T urb=1 at the turbidity maximum for medium dry periods at river km 20 (Meire et al., 2005); T urb=0.6 upstream (calibrated)).Note that the notation [Z] signifies the concentration of species Z.  : Millero (1995)), converted to free proton scale, pressure corrected according to Millero (1995) using the mean estuarine depth for each model box, and converted to volumetric units using the temperature and salinity dependent formulation for seawater density according to Millero and Poisson (1981).
To simulate light attenuation for primary production, a maximal reaction rate r PP is scaled by a 1 -Holling Type III (Gurney and Nisbet, 1998) dependency term for mean estuarine depth D and for relative turbidity T urb (T urb=0 downstream; T urb=1 at the turbidity maximum for medium dry periods at river km 20 (Meire et al., 2005); T urb=0.6 upstream (calibrated)).Furthermore primary production depends on a Monod term for DIN (Dissolved Inorganic Nitrogen= ).The fraction of ammonium that is used by primary production p PP NH + 4 is calculated with a Monod term as well.Primary production is assumed to produce reactive organic matter (FastOM) only.
The temperature dependency for all processes is modelled with a Q 10 formulation using a standard Q 10 value of 2, all rates are expressed at a standard base temperature of 15 , and the upstream relative turbidity have been calibrated such that the primary production in our model fits the combined carbon fixation activity of the groups of organisms given in Soetaert and Herman (1995c) In any natural aqueous system, and in saline systems in particular, a certain set of chemical acid-base reactions has to be taken into account if pH is to be modelled.Due to their fast reaction rates compared to all other modelled processes, these reversible acid-base reactions are considered to be in local equilibrium at any time and at any point in the estuary (Stumm and Morgan, 1996;Hofmann et al., 2008).For the model presented here, the set of acid-base equilibria given in Table 3 has been chosen.

Air-water exchange
While N 2 is assumed to be instantaneously lost to the atmosphere, O 2 and CO 2 are exchanged with the atmosphere according to a formulation given in Thomann and Mueller (1987): Saturation concentrations are calculated using Henry's law (e.g.Atkins, 1996): with f C (atm) being the fugacity of C, K 0 C (mmol (kgsoln atm) −1 ) being the Henry's constant for C, and ρ SeaWater (kg-soln m −3 ) being the temperature and salinity dependent density of seawater.Henry's constants are calculated according to Weiss (1974) for CO 2 and based on Weiss (1970) for O 2 .Both formulations can be found in appendix A.
The atmospheric fugacities of CO 2 and O 2 are assumed to be the same as their atmospheric partial pressures.f CO 2 over the Scheldt waters is assumed to be 383 µatm (averaged from Scheldt area specific values for pCO 2 (partial pressure) given in Borges et al., 2004b), f O 2 is assumed 1 to be 0.20946 atm as given by Williams (2004).The temperature and salinity dependent density of seawater is calculated according to Millero and Poisson (1981).
A variety of different empirical relationships between K L C and wind speed have been proposed in the literature (e.g.Wanninkhof, 1992;Borges et al., 2004b;Raymond and Cole, 2001;McGillis et al., 2001;Clark et al., 1995;Liss and Merlivat, 1986;Kuss et al., 2004;Nightingale et al., 2000;Banks and Herrera, 1977;Kremer et al., 2003).All these relationships have been implemented and tested.However, none of these formulations yielded a model fit significantly better than with a constant piston velocity.To keep the model as simple as possible and considering the high degree of uncertainty associated with K L C -wind speed relations, we used a constant K L C =2.7 cm h −1 (calibrated) in our model.
Inclusion of air sea exchange of NH 3 has been tested but discarded since it did not change the model performance.

Advective-dispersive transport
Since this model focuses on a period of several years, tidally averaged one-dimensional advective-dispersive transport of substances (all modelled substances/chemical species are assumed to be dissolved) is assumed.This transport approach incorporates longitudinal dispersion coefficients which parametrise a variety of physical mechanisms, including effects of either vertical or horizontal shear in tidal currents (Monismith et al., 2002).Advective-dispersive transport Tr C of chemical species C is expressed as given in Thomann and Mueller (1987) and used amongst others by Soetaert and Herman (1995b) and Ouboter et al. (1998) for the Scheldt estuary: The cross sectional area of the channel A (m 2 ), the tidal dispersion coefficient E (m 2 s −1 ) and the advective flow Q (m 3 s −1 ) are functions of the position x along the estuary.Due to the lack of experimental data for E in the Scheldt estuary, a relationship between local water depth D and E has been developed as part of this work.
The model has been spatially discretised according to a finite differences approach given in Thomann and Mueller (1987).Equation ( 3) can be spatially discretised for concentrations [C] i in mmol m −3 of any modelled species C in spatial model box i, describing the flow of matter across a set of model boxes with homogeneous contents and constant volume over time.This is done by applying a first order centred differencing scheme with a step-width of half a model box to the dispersive term, first to the "outer" derivative, then to the obtained "inner" derivatives.In the discretisation of the advective term, the same centred differencing scheme is used for the flux Q, while a backwards differencing scheme 2 2 Using the same centred differencing scheme for concentrations as well would mean calculating concentrations at the boundary of model boxes.This can lead to a non mass conservative behaviour as the example of a zero concentration in the previous model box and a non-zero one in the current model box shows.
is used for the concentration [C].This approach leads to a transport formulation of: and Q i−1,i (m 3 s −1 ) being the flow over the interface between box i−1 and i; E i−1,i (m 3 s −1 ) the bulk dispersion coefficient at the interface between box i−1 and i; E i−1,i (m 2 s −1 ) the tidal dispersion coefficient at the interface between box i−1 and i; A i−1,i (m 2 ) the cross sectional area of the interface between box i−1 and i; x i (m) the length of model box i; and x i−1,i (m) the length from the centre of box i−1 to the centre of box i.

Model state variables and their rates of change
The chemical equilibrium reactions happen on a much faster timescale than the biogeochemical and physical processes (Zeebe and Wolf-Gladrow, 2001).To avoid numerical instabilities while keeping the solution of the model computationally feasible, a set of model state variables that are invariant to the acid-base equilibria has been devised.This was done using the transformation into the canonical form and operator splitting approach (Hofmann et al., 2008).The resulting state variables of the model (kinetic species and equilibrium invariants) are given in Table 4 (definitions) and Table 5 (rates of change).Note that N 2 is neglected, salinity (S) and dissolved organic carbon (DOC) are added to the list of state variables, and that organic matter ((CH 2 O) γ (NH 3 )) is split into a reactive (FastOM) and a refractory part (SlowOM).
Note further that "organic matter" ((CH 2 O) γ (NH 3 ), Fas-tOM, SlowOM) refers to particulate organic matter.Dissolved organic matter exhibits quasi conservative mixing in the Scheldt estuary (Soetaert et al., 2006) with a C/N ratio γ DOM of 13.5 (average value for Scheldt estuary dissolved organic matter, T. van Engeland personal communication) and is conservatively modelled as [DOC].Note that particulate matter is transported in the same manner as other dissolved state variables.

pH pH, or the free proton concentration [H +
], is modelled since it can be used as a master variable to monitor the chemical state of a natural body of water, as almost any biogeochemical process occurring in such an environment affects ] either directly or indirectly (Stumm and Morgan, 1996;Soetaert et al., 2007).The free pH scale is used here, since all components of the systems in question are considered explicitly, including [HF] and [HSO − 4 ] (cf.Dickson, 1984).
www.biogeosciences.net/5/981/2008/Biogeosciences, 5, 981-1006, 2008 [FastOM] mmol N m −3 fast decaying particulate organic matter: salinity equilibrium invariants: Our model contains a pH calculation routine as described by the "solution method 3b" given in Hofmann et al. (2008) which has been inspired by Luff et al. (2001) and Follows et al. (2006).Since pH data for the Scheldt estuary have been measured on the NBS pH scale (Durst, 1975), the modelled free scale pH was converted to the NBS scale using the activity coefficient for H + , which is calculated by means of the Davies equation (Zeebe and Wolf-Gladrow, 2001).The use of the Davies equation is assumed to be a sufficient approximation, since according to Zeebe and Wolf-Gladrow (2001) it is valid up to ionic strengths of approximately 0.5 and yearly averaged salinity values at the mouth of the Scheldt estuary (the downstream boundary of our model) are about 28 resulting in a ionic strength of approximately 0.57, while in the remaining estuary salinity values, hence ionic strengths, are lower.

Data
All data that are referred to as monitoring data were obtained by the Netherlands Institute of Ecology (NIOO) for 16 stations in the Scheldt between Breskens/Vlissingen (The Netherlands) and Rupelmonde (Belgium) by monthly cruises of the NIOO RV "Luctor".River kilometres (distance from the start of the model domain at Rupelmonde) and locations of the stations can be found in Table 6, a map of the Scheldt estuary indicating the positions of the sampling stations is given in Fig. 1.

Boundary condition forcings
Values measured at monitoring station WS1 are considered to be downstream boundary conditions and values of station WS16 to be upstream boundary conditions.For [DOC] and salinity, monitoring data were used for upstream and downstream boundaries.Organic matter (OM) is partitioned in a reactive (FastOM) and a refractory (SlowOM) fraction.At the downstream boundary a fraction of 0.4 of the organic matter (Soetaert and Herman, 1995a) upstream: 4700 mmol m −3 ; downstream: 2600 mmol m −3 ; constant for all modelled years) have been calibrated using data from the year 2003.For [ HSO − 4 ], [ HF] and [ B(OH) 3 ] they were calculated from salinity monitoring data using formulae given in appendix B (DOE, 1994).To ensure consistency, total alkalinity boundary concentrations were calculated from pH, salinity and temperature monitoring data and [ CO 2 ] boundary conditions according to the total alkalinity definition used in our model (Dickson, 1981).

Physical condition forcings
River flows at the upstream boundary are obtained from the Ministry of the Flemish Community (MVG).Since the advective flow increases from the upstream boundary towards the north sea due to inputs by amongst others the Antwerp harbour and the channel Gent-Terneuzen (van Eck, 1999;Soetaert et al., 2006) and no downstream flow values for the modelled years were available, data from the years 1980 to 1988 obtained from the SAWES model (van Gils et al., 1993;Holland, 1991) have been used to calculate a flow profile along the estuary.This is done by calculating percentages of mean flow increase between SAWES model boxes from the 1980 to 1988 data and scaling the upstream-border data for the years 2001 to 2004 accordingly, implementing a total flow increase of around 45% from Rupelmonde to Vlissingen.This practice implies that the lateral input into a particular part of the estuary has the same concentration as the input from upstream of this area.
The average water depth D (Fig. 2) was obtained for 13 MOSES model boxes from Soetaert and Herman (1995b) and interpolated to centres and boundaries of the 100 model boxes used in this work.The cross sectional area A (Fig. 2), has been obtained as a continuous function of river kilometres from Soetaert et al. (2006).Temperatures of the water column are monthly monitoring data and range from around freezing in winter to roughly 25 • C in summer.Along the estuary there is a small spatial gradient in yearly temperature averages from around 13 • C close to Rupelmonde to 12 • C at Vlissingen.At 40 km from Rupelmonde, roughly at the position of the entrance to the harbour of Antwerp (Zandvliet lock) the yearly temperature averages show a pronounced maximum around 1 • C higher than values at the upstream boundary.This is probably due to a combined warming effect of the Antwerp harbour and the nuclear power plant at Doel roughly 4 km upstream of the temperature maximum.

Implementation and calibration
The spatial dimension of the model area along the estuary from Rupelmonde to Vlissingen is discretised by means of 100 supposedly homogeneous model boxes of 1.04 km length.These model boxes are assumed to have a tidally averaged volume (constant over time) and are numbered from i=1 (upstream) to i=100 (downstream).
The model has been implemented in FORTRAN using the modelling environment FEMME 3 (Soetaert et al., 2002) and numerically integrated over time with an Euler scheme using a time-step of 0.00781 days.Seasonal dynamics for the four model years were resolved but only yearly averages will be presented.Most post-processing of model output and creation of figures has been done with the R statistical computing environment (R Development Core Team, 2005).Boundary and physical conditions (weather) are forced onto the model with monthly data from the respective years as described in Sect.2.2.6.
An artificial spin-up year has been created, starting with arbitrary initial conditions for all state variables and containing only the initial values for the year 2001 for each forcing function.After running this spin-up year to create suitable initial conditions, the years 2001 to 2004 were run consecutively using forcing data from the respective year.Since FEMME requires a value for each forcing function at the beginning and at the end of each year, data for the first and the last day of each year were calculated by interpolation.For these calculations also data from 2000 and 2005 were used.
After the fitting of primary production in the year 2003 to values from Soetaert and Herman (1995c)

A relation between mean water depth D and tidal dispersion coefficients E
Following the idea put forward by Monismith et al. (2002) that there is proportionality between tidal dispersion and water depth and by comparing model output to field data for the conservative tracer salinity (see Fig. 3), we devised a linear relationship between tidal dispersion coefficients E i−1,i and mean water depth D i for model box i: with Values for E Max and E Min were calibrated (within the range of E values given in Soetaert and Herman, 1995b), while values for D Max and D Min were obtained for the 13 MOSES model boxes from Soetaert and Herman (1995b).

PSfrag replacements
River km (distance from Rupelmonde) E values obtained with our formula are comparable to the ones used by Vanderborght et al. (2007) at the respective grid size.
While transport coefficients were calibrated for the year 2003, yearly averaged longitudinal salinity profiles could be reproduced for the years 2001, 2002, and 2004 (see Fig. 3) by imposing respective advective flows Q and boundary conditions for salinity.
In all years, salinity more than linearly increases from around 1 upstream to values between 14 and 22 at km 60 and subsequently increases in linear fashion to values between 26 and 30 at the downstream boundary.Figure 5 compares the modelled nitrification rate (R Nit ) for the year 2003 with field data obtained in the same year by Andersson et al. (2006) with the 15 N method.This is an independent model validation as these data have not been used in any way for calibration.It shows excellent agreement between measured and modelled nitrification rates, with an initially more than but overall approximately linear decline in nitrification from 13.1 mmol m −3 d −1 to 0.4 mmol m −3 d −1 from km 0 to km 60, followed by a gradual decline to 0.12 mmol m −3 d −1 in the most downstream model box.4 consumption upstream and 41% downstream.The budget for [ CO 2 ] (Fig. 6b; Table 8) is characterised by CO 2 loss to the atmosphere via air-water exchange (E CO 2 ), advective-dispersive CO 2 input (Tr CO 2 ), as well as CO 2 production by oxic mineralisation (R OxCarb ) and denitrification (R DenCarb ).Primary production accounts for only 10% of the ≈6800 mmol C m −3 y −1 CO 2 consumption in the upstream region, while its relative importance increases in the downstream area where it accounts for almost 50% of the approximately 700 mmol C m −3 y −1 of CO 2 consumption.CO 2 production by oxic mineralisation (R OxCarb ) and denitrification (R DenCarb ) steadily decrease from values around 3700 and 2300 mmol C m −3 y −1 to values around 700 and 30 mmol C m −3 y −1 , where the relative importance of oxic mineralisation increases from 55% to 96% and that of denitrification decreases from 34% to 4%.The net advective-dispersive CO 2 input (Tr CO 2 ) accounts with about 800 mmol C m −3 y −1 for 12% of the CO 2 input at the upstream boundary, has a local maximum at river km 22 with ≈3300 mmol C m −3 y −1 , and has another one at river km 48 with ≈2700 mmol C m −3 y −1 , where it accounts for 70% of the total CO 2 input.At river km 60, Tr CO 2 exhibits a local minimum with ≈200 mmol C m −3 y −1 , followed by another local maximum at km 67 with ≈800 mmol C m −3 y −1 and, finally, a steady decrease reaching negative values at the downstream border.

Volume integrated budgets
As the estuarine cross section area increases from 4000 m 2 upstream to around 76 000 m 2 downstream, there is a much larger estuarine volume in downstream model boxes than in upstream model boxes.Thus, volume integrated production or consumption rates (rates "per river km") are similar in the upstream and the downstream part of the estuary (in accordance with findings of Vanderborght et al., 2002), unlike in the volumetric plots, where the upstream region was clearly dominant.Figures 8 and 9 show volume integrated budgets along the estuary averaged for the years 2001 to 2004.While for NH + 4 (Fig. 8a) the upstream region can be identified as most important in terms of total turnover, for CO 2 (Fig. 8b) and O 2 (Fig. 9a) a pronounced maximum of total turnover can be distinguished at around km 50 (intertidal flat system of Saeftinge) next to high total turnover values in the downstream region from around river km 70 on.For NO − 3 the upstream area is most important, similar to its volumetric budget.
Figure 8a and Table 11 in combination with the percentages given in Table 7 show that, in contrast to the volumetric budget, volume integrated NH + 4 production due to oxic mineralisation (R Ox ) is maximal at the downstream border with ≈12.2 mmol N km −1 y −1 .All other processes still contribute most to [ NH + 4 ] at the upstream boundary, ex-cept for ammonium consumption due to primary production which also has its volumetric maximum in the downstream region.In the volume integrated plot, it can be clearly seen that the net influence of advective dispersive transport on ) is negative at the downstream boundary (≈−4.6 mmol N km −1 y −1 ).The volume integrated budget for CO 2 (Fig. 8b; Table 11 in combination with Table 8) shows a distinct maximum of total CO 2 turnover at km 48 with a total CO 2 production of 88.5 mmol C km −1 y −1 and a total CO 2 consumption of 90.2 mmol C km −1 y −1 .This maximum is followed by a slight dip in total turnover at km 60.At river km 86, total CO 2 turnover rates reach values slightly lower than at the maximum of km 48 and gradually decline towards the downstream boundary where 60% of the maximal values are reached.Around the downstream boundary, CO 2 production by oxic mineralisation has its maximum at ≈50.6 mmol C km −1 y −1 while volume integrated CO 2 consuption due to primary production has its maximum at km 67 with ≈43.4 mmol C km −1 y −1 .Especially at kilometres 48 to 67 it can be seen that the volume integrated total loss of CO 2 is higher than the total CO 2 production.River km (distance from Rupelmonde) River km (distance from Rupelmonde) The volume integrated budget for O (Fig. 9a; Table 11 in combination with Table 9) shows similar process patterns and basic shape as the volume integrated budget for CO 2 , with a maximum of absolute values of total O 2 turnover at km 48 followed by a dip in absolute values of total O 2 turnover at km 60, with the exception of Tr O 2 changing sign twice.
Unlike for the other chemical species, the volume integrated budget for NO − 3 (Fig. 9b; Table 11 in combination with Table 10) shows distinct maxima in total turnover at the upstream boundary and decreases towards the downstream boundary.However, the resulting trumpet-like shape is not as pronounced as for the volumetric budget for [NO − 3 ].

Whole estuarine budgets
Figures 10 and 11 show budgets of NH + 4 , CO 2 , O 2 , and NO − 3 production and consumption, integrated over the whole model area and one year (averaged over the four modelled years).

PSfrag replacements
Mmol N river-km −1 River km (distance from Rupelmonde)  and advective-dispersive transport (≈0.5 Gmol).Denitrification plays a minor role producing less than 0.1 Gmol of NH + 4 per year.Total estuarine CO 2 (Fig. 10b) is prominently influenced by out-gassing of CO 2 to the atmosphere, which consumes ≈3.3 Gmol, while primary production consumes only ≈2.0 Gmol.Advective-dispersive input and oxic mineralisa-tion supply CO 2 at values of ≈2.2 Gmol and ≈2.7 Gmol respectively.Again, denitrification has a minor influence, producing less than 0.3 Gmol of inorganic carbon per year.

Interannual differences
Although the freshwater discharge Q of the Scheldt shows no consistent trend during the years 1990 to 2004, it peaked in 2001 and fell rapidly until 2004 (Meire et al., 2005;Van Damme et al., 2005;van Eck, 1999) River km (distance from Rupelmonde) ) R PP This decrease in Q is most likely the reason for several observed trends in our model.E.g. salinity S increases from 2001 to 2004 which is clearly a result of lowered freshwater discharge, as more saline seawater can enter the estuary (Meire et al., 2005;Van Damme et al., 2005).Similarly, the observed decrease in [ CO 2 ] can be explained by less CO 2 being imported into the estuary from the river at lower freshwater discharge ([ CO 2 ] not shown).A decrease in [ CO 2 ] also means a decrease in [CO 2 ].However, via its influence on the dissociation constants K * , salinity also influences [CO 2 ], higher salinity meaning lower [CO 2 ], reinforcing its downward trend from 2001 to 2004.Decreasing levels of [CO 2 ] lower the CO 2 saturation state of the water and eventually lead to less CO 2 export to the atmosphere (The total amount of CO 2 export to the atmosphere and the volume averaged saturation states for the four

PSfrag replacements
Gmol O 2 estuary −1 Furthermore, pH influences [CO 2 ], higher pH implying lower [CO 2 ].There was an upward trend in pH during our model years, which reinforces the downward trend in CO 2 export to the atmosphere.We believe that there is a relation between the decrease in freshwater flow and the upward trend in pH from 2001 to 2004, however, the exact mechanism for this relationship is not straight-forward.A mechanistic model able to quantify the influence of different modelled kinetic processes on the pH (Hofmann et al., 2008) will shed some further light on this issue.

Model performance: data-model validation
Our objective was to construct a simple model reproducing observed data and interannual differences, allowing the establishment of annual budgets.Advective-dispersive transport is accurately reproduced (confer the fit of yearly averaged logitudinal profiles of the conservative tracer salinity, Fig. 3).
The model also reproduces  This level of performance has been achieved by using 7 biochemical parameters from literature and calibrating the other 5 using data for the year 2003 (plus the calibration of 5 biochemical parameters to fit primary production to values from Soetaert and Herman, 1995c).To run the model for the years 2001, 2002, and 2004, the upstream advective forcing, the temperature forcing, and the boundary concentrations have been adapted, but no further calibration was involved, making the fits for those years a model validation.Furthermore the model has been independently validated by comparing nitrification rates to field data from the year 2003 (Andersson et al., 2006) realising the objective of creating a tool to examine interannual differences and annual budgets of key chemical species in the Scheldt estuary.
The very good overall performance of our rather simple model confirms the notion of Arndt et al. (2007) that for estuaries biogeochemical model complexity can be kept low as long as physical processes (i.e.transport) perform sufficiently well.(Note that a good model performance in terms of physics does not imply a complex representation of physical processes.)This argument is strengthened by the fact that the fits of a model including NH 3 air-sea exchange, sulfate reduction, sulfide oxidation, biogenic calcification as well as calcite and aragonite precipitation did not significantly differ from the fits of the model presented here.However, some features of the fits of state variables given in Fig. 3 suggest that the inclusion of further processes might make the model even more accurate.Modelled [ NH + 4 ] for example is slightly too high between river km 30 and km 50, which could be explained by microbial ammonium uptake (cf.Middelburg and Nieuwenhuize, 2000).Furthermore, a reason for an underestimated organic matter concentration around the intertidal flat area of Saeftinge (≈river km 35 to 50) might be the fact that the model does not include any explicit organic matter input from this ecosystem consisting mainly of vascular plants.Although these and other arguments can be made about details, we consider the fit of our model adequate for our aims and for the simplified conceptual model.Denitrification is not strongly constrained in our model because of the lack of direct measurements.
Figure 13 shows the model fit for [NO − 3 ] with three different parametrisations for denitrification: our parametrisation based on literature values (solid black line), no denitrification (dashed red line) and denitrification maximised (dotted blue line) by using a very small k NO − 3 (10 −8 mmol N m −3 ) and a very large k Inh O 2 (10 8 mmol N m −3 ), resulting in ≈1.3 Gmol NO − 3 consumption due to denitrification per year (compared to 0.2 Gmol y −1 with our parametrisation).Although the effect of denitrification on [NO − 3 ] along the estuary is small, this plot shows that our parametrisation based on literature values gives the best model performance.Andersson (2007) reports 16 mmol m −2 d −1 of NO − 3 consumption due to denitrification for sediment from one lower estuary location in the Scheldt.Considering 338 km 2 of estuarine surface area, our model result of 0.2 Gmol NO − of this number.Even with maximal denitrification and the associated worse model performance (Fig. 13), only a denitrification of 10.5 mmol NO − 3 m −2 d −1 can be achieved.This difference between our top-down whole estuarine estimate and the estimate of Andersson (2007) clearly shows the difficulties one encounters when upscaling sedimentary biogeochemical process rates over large areas and long timescales.
Furthermore we model denitrification pelagically as a proxy for a sedimentary process.As a result of this practice the effect of denitrification on the concentrations in the water column is not geographically fixed but moves up and down the river.This would have a large effect on a tidally resolved model as parts of the sediment may fall dry and loose contact with the water column.However, for a tidally averaged model as ours and for yearly averaged values in particular, this effect is rather small.Nevertheless, we are aware that modelling pelagic denitrification as a proxy for benthic processes introduces an error in our volumetric and volume integrated logitudinal estuarine budget.But since the region of main denitrification activity always stays within our model domain, our total stock budgets are not significantly affected by this error.

Volumetric budgets for [ NH
According to the volumetric budgets for

and [NO −
3 ], the estuary can be divided into two parts, an upstream region (river km 0 until roughly km 55 at Walsoorden) and a downstream area, with process rates being up to an order of magnitude higher in the upstream part.This division, which is in accordance to findings of Vanderborght et al. (2007), Regnier et al. (1997), andSoetaert andHerman (1995a,c), is more pronounced for [NO − 3 ] and [ NH + 4 ] than for [O 2 ] and [ CO 2 ], because the latter two quantities also depend on gas exchange with the atmosphere which in turn only depends on the concentration of the species to be exchanged and not on the concentrations of other species.Biogeochemical rates, in contrast, usually depend on the concentrations of several substances.Therefore the rates of gas exchange processes decrease less than biogeochemical rates as nutrient levels decrease when going from upstream to downstream in the Scheldt estuary.However, E CO 2 is also strongly influenced by the pH in the estuary, which rises from ≈7.6 in the upstream regions to ≈8.1 downstream.This is because E CO 2 depends linearly on [CO 2 ] which in turn depends on pH.A rising pH induces a declining [CO 2 ].Considering the estuarine averages of [ CO 2 ]≈3400 mmol m −3 , T≈13 • C, and S≈10, [CO 2 ] declines by a factor of about 3.25 from around 130 mmol m −3 to around 40 mmol m −3 for pH increasing from around 7.6 to 8.1.This means, that the increase in pH along the estuary contributes to the decline in the absolute values of E CO 2 from upstream to downstream.
While the volumetric budgets for [ NH + 4 ] (Fig. 6a) and [NO − 3 ] (Fig. 7b) show relatively smooth features, the volumetric budget for [ CO 2 ] (Fig. 6b) exhibits a more spiked pattern.Its shape is a result of the dependency of [ CO 2 ] on CO 2 air-water exchange (E CO 2 ) which in turn is influenced by the estuarine depth D. This fact makes estuarine depth patterns visible in the budget for [ CO 2 ].To a lesser extent the same can be seen in the volumetric budget for [O 2 ] (Fig. 7a), which in the upper reaches of the estuary, between river km 0 and about 22, smoothly follows the [ NH + 4 ] budget (dominated by nitrification), but in the downstream part of the estuary also shows estuarine depth patterns, attributable to its dependency on O 2 air-water exchange (E O 2 ).

Oxygen budget
It is clear that oxygen consumption by nitrification and oxic mineralisation, and export by advective-dispersive transport processes are balanced not only by oxygen production due to primary production but mainly by oxygen import from the atmosphere.Similarly to Ouboter et al. (1998), in our model about one third of the total oxygen consumption in the estuary is due to nitrification; oxygen consumption by oxic mineralisation is comparable to oxygen consumption by nitrification, similar to conditions reported by Regnier et al. (1997).To be noted, however, is that integrated over the whole estuary, oxic mineralisation consumes more oxygen per year than nitrification which is in accordance with experimental findings of Gazeau et al. (2005).In accordance with its heterotrophic nature (Gazeau et al., 2005), the estuary is a net consumer of oxygen, only about 43% of the oxygen entering the estuary by advection-dispersion and reaeration leaves the estuary at the mouth.86% of the total O 2 input into the estuary is due to re-aeration.1986.This is most likely due to reduced riverine nutrient nd organic loadings and resulting lower [ NH +  4 ] in the ears 2001 to 2004 as compared to 1980 to 1986 (Soetaert t al., 2006).Due to lower [ NH +  4 ] in 2001 to 2004, voluetric nitrification rates in the upstream region were 77%, at m 60 roughly 11%, and in the downstream region roughly % of the values in the early eighties as reported by Soetaert nd Herman (1995a).Due to the large estuarine volume in e downstream region, the drop in nitrification in this area mostly responsible for the drop in total ammonium conumption by nitrification in the whole estuary from the early ighties to our model time period (2001 to 2004).
This downward trend in total nitrification shows that the itial intensification of nitrification in the Schelde due to inreasing oxygen levels since the second half of the seventies Van Damme et al., 2005;Soetaert and Herman, 1995c) has eaked and subsequently decreased again, most likely due to educed ammonium concentrations in the estuary (Soetaert t al., 2006), suggesting a shift from initial oxygen limitation f nitrification towards ammonium limitation now.A simir shift as has happened in time can be observed longitudially during our model time period, as nitrification is oxyen limited upstream and ammonium limited downstream.owever, as in the seventies (Billen et al., 1985) (Vanderborght et al., 2007), ary, as is advectively imported (confer Fig. 11b).Yet, the nitrate producing character of the Schelde estuary diminished due to reduced nitrification rates relative to the early eighties as Soetaert and Herman (1995a) reported that three times as much nitrate was exported to the sea as entered the estuary in the eighties.
Only 10% of the total N input in the system (nitrate, ammonium and particulate and dissolved organic nitrogen together) is lost to the atmosphere as N 2 due to denitrification (Fig. 14), while the rest is exported to the Southern Bight of the North Sea.This shows a clear downwards trend in the percentage of N 2 production, as in the eighties 21-25% of the total nitrogen imported into the Schelde (Soetaert and Herman, 1995a;Ouboter et al., 1998), and in the seventies around 40% (Soetaert and Herman, 1995a;Billen et al. 1985) of the total nitrogen import into the estuary was removed within the estuary, mainly due to denitrification.This phenomenon is likely due to improved oxygen conditions in the Schelde from the seventies until now (Soetaert and Herman, 1995a;Soetaert et al., 2006), moving the zone of denitrification more into the narrow upstream regions and generally allowing for less denitrification because of the limited area of sediments, the prime location of denitrification.A downward trend in the percentage of N 2 generated and escaping means an upward trend in the percentage of N expor to the sea.In absolute values, however, this resulted initially in an increasing N export to the North Sea from the seventies to the eighties from 1.9 Gmol y −1 to 3.6 Gmol y −1 (Soetaer    (Soetaert et al., 2006).Due to lower [ NH + 4 ] in 2001 to 2004, volumetric nitrification rates in the upstream region were 77%, at km 60 roughly 11%, and in the downstream region roughly 5% of the values in the early eighties as reported by Soetaert and Herman (1995a).Due to the large estuarine volume in the downstream region, the drop in nitrification in this area is mostly responsible for the drop in total ammonium consumption by nitrification in the whole estuary from the early eighties to our model time period (2001 to 2004).This downward trend in total nitrification shows that the initial intensification of nitrification in the Scheldt due to increasing oxygen levels since the second half of the seventies (Van Damme et al., 2005;Soetaert and Herman, 1995c) has peaked and subsequently decreased again, most likely due to reduced ammonium concentrations in the estuary (Soetaert et al., 2006), suggesting a shift from initial oxygen limitation of nitrification towards ammonium limitation now.A similar shift as has happened in time can be observed longitudinally during our model time period, as nitrification is oxygen limited upstream and ammonium limited downstream.However, as in the seventies (Billen et al., 1985), eighties (Van Damme et al., 2005;Soetaert and Herman, 1995c;Regnier et al., 1997), and nineties (Vanderborght et al., 2007) Furthermore, the estuary is a net producer of nitrate in the year 2001 to 2004.Roughly 1.4 times as much nitrate leaves the estuary by advection-dispersion at the mouth of the estuary, as is advectively imported (confer Fig. 11b).Yet, the nitrate producing character of the Scheldt estuary diminished due to reduced nitrification rates relative to the early eighties, as Soetaert and Herman (1995a) reported that three times as much nitrate was exported to the sea as entered the estuary in the eighties.
Only 10% of the total N input in the system (nitrate, ammonium and particulate and dissolved organic nitrogen together) is lost to the atmosphere as N 2 due to denitrification (Fig. 14), while the rest is exported to the Southern Bight of the North Sea.This shows a clear downwards trend in the percentage of N 2 production, as in the eighties 21-25% of the total nitrogen imported into the Scheldt (Soetaert and Herman, 1995a;Ouboter et al., 1998), and in the seventies around 40% (Soetaert and Herman, 1995a;Billen et al., 1985) of the total nitrogen import into the estuary was removed within the estuary, mainly due to denitrification.This phenomenon is likely due to improved oxygen conditions in the Scheldt from the seventies until now (Soetaert and Herman, 1995a;Soetaert et al., 2006), moving the zone of denitrification more into the narrow upstream regions and generally allowing for less denitrification because of the limited area of sediments, the prime location of denitrification, in the upstream regions.A downward trend in the percentage of N 2 generated and escaping means an upward trend in the percentage of N export to the sea.In absolute values, however, this resulted initially in an increasing N export to the North Sea from the seventies to the eighties from 1.9 Gmol y −1 to 3.6 Gmol y −1 (Soetaert and Herman, 1995a), before it decreased to 2.5 Gmol y −1 in our model period (value incl.DON export).This decrease in N export is due to approximately halved input loads (4.7 Gmol y −1 in the eighties (Soetaert and Herman, 1995a), 2.3 Gmol y −1 in our model).Yet the N export to the sea in 2001 to 2004 is still higher than in the seventies, in spite of the fact that the input into the system was much higher then (3.7 Gmol y −1 , Billen et al., 1985).Table 13 summarises these numbers.Note that these results should be considered tentative, since denitrification is poorly constrained by data and in our model (Sect.4.1.1).

Carbon
Estuaries are a significant source of CO 2 to the atmosphere and are even important on a global scale (Borges et al., 2006).Our model suggests an averaged CO 2 export to the atmosphere of 3.3 Gmol y −1 , which is about 13% of the total carbon input (including DOC) into the system (Fig. 15).This value is much lower than values reported by others (Table 14).These differences represent both a true trend in time as well as uncertainties in the estimations.Frankignoulle et al. (1998), Hellings et al. (2001), and Vanderborght et al. (2002) use data from the early to mid nineties of the 20th century.pCO 2 levels in the Scheldt estuary considerably dropped from the nineties to our model time period (Borges and Middelburg, unpublished data), entailing a true downwards trend in CO 2 degassing.One reason for this might be decreasing nitrification rates resulting in a higher pH amd thus lower pCO 2 .
However, there are several sources of uncertainties associated with estimating CO 2 degassing which could contribute in part to the differences shown in Table 14: 1. differences in the riverine discharge between the estimation periods and resulting carbon import into the system (confer a factor 2 difference in E CO 2 from 2001 and 2004 with an approximately halving freshwater discharge), 2. differences in estimates of the estuarine surface area which contributes to CO 2 air-sea exchange (see also Borges et al., 2006), 3. former overestimations of the piston velocity parametrisation mostly due to overestimated gas exchange flux measurements with the floating chamber method (Raymond and Cole, 2001), and 4. the use of budgeting approaches without the employment of mechanistic models with rigorous mass conservation (e.g.upscaling of discrete flux measurements).Frankignoulle et al. (1998) reported that in the Scheldt estuary two thirds of the CO 2 flux to the atmosphere results from heterotrophy and only one third from ventilation of riverine DIC and Abril et al. (2000) estimate that riverine DIC contributes only 10% to the total CO 2 degassing of the Scheldt estuary.In contrast, our model suggests that, averaged over the four modelled years, 67% of the CO 2 export to the atmosphere can be attributed to ventilation of riverine DIC (based on results of model runs with and without biogeochemical processes).

Summary
The Scheldt estuary is an active biogeochemical reactor where, all along its path, the transformations occur at a similar magnitude.This is due to the combined effects of very high volumetric process rates and a small estuarine volume upstream, giving way to gradually decreasing volumetric rates and an increasing estuarine volume downstream.
With respect to the nitrogen cycle, the estuary has evolved towards a more and more neutral passage way of total nitrogen, with only 10% of the total imported nitrogen being lost from the estuary.This is in sharp contrast to the situation in the eighties and seventies where this loss amounted to more than 20% and 40%, respectively.Coinciding with a reduced total nitrogen import in the estuary, this has led to a current Note that the budget is not fully closed.There is an overall loss term of 0.12 Gmol-C y −1 .This is consistent with poral downwards trend in [ P CO 2 ].Note also that organic carbon (Org-C) refers to particulate organic matter (OM).Dissolved orga on (DOC) is not shown in this budget.Around 2.2 Gmol of conservatively modelled DOC enters and leaves the estuary on average r.
rly linked to the decreased freshwater discharge, that was ved in that period.

pendix A Gas exchange constants (f (T ,S))
pirical formulations for the temperature and salinity dedency of the gas exchange constants used here, can be ught into the generic form: 8.314 Nm (K•mol) −1 for the gas constant R and an ambie pressure of 101 300 N m −2 .The expression for the Henr constant has then been created by dividing the expressi for the saturation concentration by an atmospheric oxyg fugacity of f O 2 =0.20946 atm (Williams, 2004).

Appendix B Total concentrations of seawater comp nents (f (S))
In DOE (1994) a table of concentrations of seawa components relative to chlorinity is given that allows infer formulations for total concentrations of all seawa components.
The formulae for the concentrations needed here, rewritt as functions of salinity with S=1.80655Note that the budget is not fully closed.There is an overall loss term of 0.12 Gmol C y −1 .This is consistent with the temporal downwards trend in [ CO 2 ].Note also that organic carbon (Org-C) refers to particulate organic matter (OM).Dissolved organic carbon (DOC) is not shown in this budget.Around 2.2 Gmol of conservatively modelled DOC enters and leaves the estuary on average per year.
total nitrogen export to the North Sea which is smaller than the export in the eighties but still larger than what was observed in the seventies.However, the estuary has a large effect on the nitrogen speciation, especially by transforming ammonium and, indirectly, the imported organic nitrogen to nitrate.Similarly to previous budget estimates, nitrification remains one of the most important oxygen consuming processes in the estuary.
The loss of imported carbon in the estuary amounts to about 13%, and occurs through physical ventilation of CO 2 to the atmosphere.Two thirds of this lost C is riverine-borne DIC, one third of the ventilated CO 2 originates from heterotrophic production in the estuary itself.Whilst the estuary remains a significant source of CO 2 to the atmosphere, our results suggest that there is a downward trend in CO 2 degassing from the nineties to our model time period.This result should be considered tentative due to the high degree of uncertainty associated with CO 2 degassing estimations.
Finally, in the four-year period (2001)(2002)(2003)(2004)) during which our model was applied, clear trends in the chemical concentrations and budgets were observed.These trends were clearly linked to the decreased freshwater discharge, that was halved in that period.

Fig. 2 .
Fig. 2. Mean depth D (black broken line) and cross sectional area A (red solid line) of the Scheldt estuary.
with [C] (mmol m −3 ) signifying the actual concentration of O 2 or CO 2 respectively, [C] sat (mmol m −3 ) being the saturation concentration of chemical species C in the water.K L C (m d −1 ) is the piston velocity, A (m 2 ) the horizontal surface area of the model box in question, V (m 3 ) the volume, and D (m) the mean water depth.
organic matter, pH (on the NBS scale) and salinity are available from the monitoring program.Data for 2003 have been used to calibrate the model and data for 2001, 2002, and 2004 have been used to independently validate it.The model does not distinguish between NO − 3 and NO − 2 and all data referring to the state variable [NO − 3 ] are the summed values of [NO − 3 ] plus [NO − 2 ].Data for organic matter (OM) have been calculated by assuming the measured percentage of nitrogen in the suspended particulate matter concentration to be available as particulate organic nitrogen.Nitrification rates for the year 2003 were obtained from Andersson et al. (2006).They have been measured with the 15 N method at salinities 0, 2, 8, 18 and 28 in January, April, July and October 2003.Yearly averaged values of these rate measurements were used to independently validate the model, as these data have not been used in model development and calibration.Total alkalinity [TA] data for the year 2003 have been obtained from Frederic Gazeau (personal communication and Gazeau et al., 2005).These yearly averages are based on four to five measurements along the estuary from January to May 2003 and monthly measurements for all 16 Scheldt monitoring stations from June to December 2003.[ CO 2 ] data for 2003 have been calculated from [TA] data, monitoring pH and salinity, temperature forcing data, and [ B(OH) 3 ], and [ NH + 4 ] from the model.
is considered FastOM, and for the upstream boundary a FastOM fraction of 0.6 has been derived Biogeosciences, 5, 981-1006, 2008 www.biogeosciences.net/5/981/2008/ (calibration of 5 biochemical parameters: r PP , k Inh DIN , k Inh D , k Inh T urb , and k NH + 4 ; and 1 boundary condition parameter: the relative turbidity upstream; fit not shown), the model was calibrated against data from the year 2003.During this model-data fitting procedure, 11 parameters were calibrated: 5 biochemical parameters, 3 boundary condition parameters, and 3 transport parameters.The transport formulation (parameters E max and E Min , see below) has been calibrated by comparing the thor or from the FEMME website: http://www.nioo.knaw.nl/ceme/femme/ model output to field data for S. [ CO 2 ] boundary conditions were calibrated by comparing the model output to [TA] data for 2003, and the parameters r Min FastOM , r Nit , p, k Inh S , o Nit as well as the fraction of OM that is considered FastOM at the upstream boundary were calibrated by comparing the model output to field data for [OM], [ NH + 4 ], [NO − 3 ], [O 2 ], and pH from the year 2003.

Fig. 3 .Fig. 4 .
Fig. 3. Fit of the biogeochemical model for 4 consecutive years.Calibration was done on data for 2003 only.Data and model output are yearly averages.

3. 2
Comparison of yearly averaged longitudinal concentration and rate profiles to measured data Several parameters have been manually calibrated to improve the fit of the biogeochemical model against field data from the year 2003 (Fig. 3).Yearly averaged values for data are time weighted averages.While biogeochemical model parameters (and transport coefficients) were fitted for 2003, the model predictions for 2001, 2002, and 2004 did not involve further calibration.Model fits to data from those years can thus be considered as model validation.It can be seen that the model reproduces the data reasonably well; not only for the calibration year 2003 but also for the validation years, with all years showing similar patterns.[ NH + 4 ] values are between 70 and 115 mmol m −3 upstream, falling almost linearly to values around 10 mmol m −3 at km 50 and staying below this value in the downstream stretches.Nitrate concentrations initially rise from concentrations approximately between 322 and 343 mmol m −3 to concentrations approximately between 325 and 348 mmol m −3 at km 16, subsequently fall more than linearly to values of 133 to 202 mmol m −3 at km 60, and then further decline quasilinearly to values between around 60 and 76 mmol m −3 in the most downstream model box.Oxygen concentrations start off at values between 64 and 91 mmol m −3 upstream, stay constant or even decline between river kilometres 0 and 20, followed by a quasi linear increase to values of ≈265 mmol m −3 at river km 60, and staying in this realm until the downstream border.The sum of the concentrations of both fractions of particulate organic matter shows larger discrepancy between model and data.It declines over the stretch of the estuary from values of 41 to 53 mmol m −3 upstream down to values of approximately 9 mmol m −3 downstream.The NBS pH shows a sigmoidal increase from values around 7.57 to 7.63 upstream to values between 8.07 and 8.12 in the downstream part of the estuary.In all years a slight dip in the order of 0.03 to 0.05 pH values can be observed between km 0 to 30 preceding the sigmoidal increase.

Figure 4
Figure 4 shows the fit of predicted [TA] and [ CO 2 ] against observed data for the year 2003 ([ CO 2 ] data calculated from [TA] data).Note that these data (from the year 2003) have been used to calibrate only the CO 2 boundary conditions for all four years.Both [TA] and [ CO 2 ] decrease in a sigmoidal fashion from upstream to downstream, [TA] from between 4430 and 4480 mmol m −3 to between 2700 and 2740 mmol m −3 , and [ CO 2 ] from around 4690 mmol m −3 to around 2600 mmol m −3 .
of the four years, model rates can be used to compile budgets for these quantities.Figures6 and 7show cumulative plots of volumetric budgets along the estuary averaged over the years 2001 to 2004.A common feature is the pronounced activity in the upper estuary, i.e. between river km 0 and 60.In this stretch of estuary, the absolute values of almost all rates decline in a quasi linear fashion to stay at low levels until the mouth of the estuary.It can be seen that [ NH + 4 ] (Fig.6a;

Fig. 10 .
Fig. 10.Total annual budgets for NH + 4 and CO 2 (b) for the whole estuary, averaged over 2001-2004.(Tr x ) up and (Tr x ) down signify advective-dispersive import or export at the upstream and downstream boundary.The error bars give the standard deviation σ , obtained from averaging over the four model years.

Fig. 11 .
Fig. 11.Total annual budgets for O 2 (a) and NH − 3 (b) for the whole estuary, averaged over 2001-2004.(Tr x ) up and (Tr x ) down signify advective-dispersive import or export at the upstream and downstream boundary.The error bars give the standard deviation, σ , obtained from averaging over the four model years.
and pH versus river kilometres very well.

Fig. 12 .
Fig. 12. Trends in key quantities from 2001 to 2004 (Q signifies the freshwater flow at the upstream boundary; S, [CO 2 ], E CO 2 : volume averaged, solid red line: total estuary, dashed blue line: upper part (to box 50), dotted green line: lower part).

4. 4
Synopsis of single species budgets: elemental budgets and comparisons with earlier decades Nitrogen and carbon budgets have been constructed for the entire estuary (Figs.14 and 15).4.4.1 NitrogenThe Scheldt estuary is a net consumer of ammonium.The total advective-dispersive input at the upstream boundary per year ((Tr NH + .F. Hofmann et al.: Nitrogen and carbon dynamics in the Schelde estuary 21 ig.14.Tentative nitrogen budget per year over the whole model area, averaged over 2001 to 2004.Values are given in Gmol y −1 and in t-N y −1 (in brackets).Note that the budget is not fully closed, there is an overall loss term of 0.13 Gmol-N y −1 .This is consistent with e temporal downward trend in [ P NH + 4 ] and [NO − 3 ].Note also that organic nitrogen (Org-N) refers to particulate organic matter (OM) issolved organic nitrogen (DON) is not shown in this budget.Using γ DOM and conservatively modelled [DOC], one can estimate tha round 0.16 Gmol of DON enters and leaves the estuary on average per year.es of process rates affecting [ NH + 4 ] in the years 2001 to 004 are at about 25% of the values during the years 1980

Fig. 14 .
Fig. 14.Tentative nitrogen budget per year over the whole model area, averaged over 2001 to 2004.Values are given in Gmol y −1 .Note that the budget is not fully closed, there is an overall loss term of 0.13 Gmol N y −1 .This is consistent with the temporal downward trend in [ NH + 4 ] and [NO − 3 ].Note also that organic nitrogen (Org-N) refers to particulate organic matter (OM).Dissolved organic nitrogen (DON) is not shown in this budget.Using γ DOM and conservatively modelled [DOC], one can estimate that around 0.16 Gmol of DON enters and leaves the estuary on average per year.
−1 exported to the North Sea 1.9 3.6 2.4 solute values of process rates affecting [ NH + 4 ] in the years 2001 to 2004 are at about 25% of the values during the years 1980 to 1986.This is most likely due to reduced riverine nutrient and organic loadings and resulting lower [ NH + 4 ] in the years 2001 to 2004 as compared to the years 1980 to 1986

Biogeosciences. 15 .
et al.: Nitrogen and carbon dynamics in the Schelde estuary 2 Tentative carbon budget per year over the whole model area, averaged over 2001 to 2004.Values are given in Gmol y −1 and C y −1 (in brackets).

Fig. 15 .
Fig. 15.Tentative carbon budget per year over the whole model area, averaged over 2001 to 2004.Values are given in Gmol y −1 .Note that the budget is not fully closed.There is an overall loss term of 0.12 Gmol C y −1 .This is consistent with the temporal downwards trend in

Table 4 .
State variables of the biogeochemical model.

Table 5 .
Rates of change of model state variables.

Table 8 .
Budget for [ CO 2 ]; values in mmol C m −3 y −1 ; percentages are of total production (positive quantities) or consumption (negative quantities), respectively.

Table 12 .
CO 2 export to the atmosphere (-E CO 2 ) in Gmol y −1 per modelled area, [CO 2 ] and [CO 2 ] sat are volume and yearly averaged values for the whole estuary in mmol m −3 .

Table 13 .
Trends in denitrification from the seventies to our model time period.Note that DON is included in the values.

Table 14 .
Export of CO 2 to the atmosphere (-E CO 2 ) in Gmol y −1 , integrated over our model area.