Factors governing the pH in a heterotrophic, turbid, tidal estuary

A method to quantify the inﬂuence of kinetically modelled biogeochemical processes on the pH of an ecosystem with time variable acid-base dissociation constants is presented and applied to the heterotrophic, turbid Scheldt estuary (SW Netherlands, N Belgium). Nitriﬁcation is identiﬁed as the main process governing the pH proﬁle of 5 this estuary, while CO 2 degassing and advective-dispersive transport “bu ﬀ er” the e ﬀ ect of nitriﬁcation. CO 2 degassing accounts for the largest proton turnover per year in the whole estuary. There is a clear inverse correlation between oxygen turnover and proton turnover. The main driver of long-term changes in the mean estuarine pH from 2001 to 2004 is a changing freshwater ﬂow which inﬂuences the pH “directly” via [ P CO 2 ] and 10 [TA] and to a signiﬁcant amount also “indirectly” via [ P NH + 4 ] and the nitriﬁcation rates in the estuary.


Introduction
The pH is often considered a master variable to monitor the chemical state of a natural body of water, since almost any process affects the pH either directly or indirectly 15 (e.g. Stumm and Morgan, 1996;Morel and Hering, 1993). This textbook knowledge has rarely been applied in studies of natural ecosystems due to limited understanding of the complex interplay of factors controlling the pH of natural waters. While current approaches do allow for modelling the pH of complex ecosystems (e.g. Boudreau and Canfield, 1988;Regnier et al., 1997;Vanderborght et al., 2002;Hofmann et al., 20 2008b), the influences and relative importances of the different physical and biological processes on the pH in those systems remain unquantified.
Especially considering the acidification of the ocean (e.g. Orr et al., 2005) and coastal seas (e.g. Blackford and Gilbert, 2007) and potential impacts of pH changes on biogeochemical processes and organisms (e.g. Gazeau et al., 2007;Guinotte and Fabry, the pH in natural aquatic systems. Estuarine ecosystems are suitable testbeds for methods quantifying the influences of especially biological processes on the pH due to their role as bio-reactors (Soetaert et al., 2006) and associated large biological influences on the pH. Mook and Koene (1975) suggested that the characteristic pH profile observed in estuaries simply re-5 sults from chemical equilibration following the mixing of freshwater and seawater. They stated that, due to the rapid increase of the dissociation constants of the carbonate system with salinity, estuaries like the Scheldt estuary (SW Netherlands and N Belgium), with high river inorganic carbon loadings and associated low riverine pH, exhibit a distinct pH minimum at low salinities. However, Mook and Koene (1975) assumed a 10 closed system and conservative mixing of total dissolved inorganic carbon ([ CO 2 ]) and total alkalinity ([TA]), i.e. they were neither considering carbon dioxide exchange with the atmosphere nor processes changing total alkalinity. Although Wong (1979) obtained reasonable agreement applying this approach to measurements in the Chesapeake Bay and it is still used to predict estuarine pH profiles (Spiteri et al., 2008), it is a rather crude approximation of reality. Whitfield and Turner (1986) showed that assuming an open system, i.e. allowing for CO 2 exchange with the atmosphere, results in significantly different pH profiles with differences up to 0.7 pH units at low salinities for systems that are fully equilibrated with the atmosphere. Furthermore, biogeochemical processes can play a significant role in influencing the pH of aquatic ecosystems 20 (e.g. Ben-Yaakov, 1973;Regnier et al., 1997;Soetaert et al., 2007). As can be seen in Fig. 1, the importance of both considering an open system with air-water exchange and including biogeochemical processes is obvious: the distinct pH minimum at low salinities found for a closed system by Mook and Koene (1975), but doubted for an open system by Whitfield and Turner (1986), can be clearly confirmed with a full biogeochem- 25 ical model. However, the relative importance of single biogeochemical processes, CO 2 air-water exchange, and transport for the pH of the system still remains unknown. Regnier et al. (1997), Vanderborght et al. (2002), and Hofmann et al. (2008b) show that reaction transport models with gas exchange, including the effects of biogeochem-199 ical processes consuming or producing protons, can reproduce the longitudinal pH profile of the Scheldt estuary fairly well. However, in none of these studies the influences of transport, CO 2 air-water exchange, and biogeochemical processes on the pH are quantified independently, because in those studies the pH was calculated using an implicit numerical approach 1 which did not allow for such a quantification. 5 While Jourabchi et al. (2005) and Soetaert et al. (2007) took steps in that direction, Hofmann et al. (2008a) present a comprehensive step by step method to set up a biogeochemical model that allows for the quantification of the influences of kinetically modelled processes (e.g. transport, CO 2 air-water exchange, biogeochemical processes) on the pH. Their direct substitution approach describes the pH evolution explicitly using an expression for the rate of change of the proton concentration over time. Hofmann et al. (2008a) introduce their explicit approach to pH modelling for systems where the dissociation constants 2 of the involved acid-base reactions are considered constant over time. However, in studies of the Scheldt estuary (e.g. Vanderborght et al., 2002;Hofmann et al., 2008b) the dissociation constants are calculated dynam-15 ically as functions of salinity and temperature to obtain reasonable pH values. While the explicit approach presented in Hofmann et al. (2008a) can be applied to a system with a spatial gradient in the dissociation constants which remains constant over time, the approach needs to be extended for application to a system where the dissociation constants vary over time, e.g. due to changes in temperature and salinity.

20
Hence, this study has four objectives: 1) the extension of the explicit pH modelling approach presented by Hofmann et al. (2008a) such that it can be applied to systems where the dissociation constants are variable over time, 2) the validation of this explicit approach by comparing predicted pH values to those obtained with an implicit approach (Hofmann et al., 2008b), 3) the quantification of proton production and consumption 25 along the Scheldt estuary by transport, CO 2 air-water exchange, and biogeochemical 1 Hofmann et al. (2008a) call this approach the operator splitting approach 2 Throughout the paper "dissociation constant" means the stoichiometric equilibrium con- 200 processes independently, given a certain freshwater flow and boundary conditions, and 4) an exploration of factors governing the mean estuarine pH in changing estuarine systems such as the Scheldt over the years 2001 to 2004.

Materials and methods
2.1 The Scheldt estuary 5 The turbid tidal Scheldt estuary is situated in the southwest Netherlands and northern Belgium (Fig. 2). The roughly 350 km (Soetaert et al., 2006) long Scheldt river drains a basin of around 21 500 km 2 (Soetaert et al., 2006) located in the northwest of France, the west of Belgium and the southwest of the Netherlands. The water movement in the Scheldt estuary is dominated by huge tidal displacements with around 10 200 times more water entering the estuary during a flood than 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 (Soetaert et al., 2006) whilst the mean water depth 15 varies quite irregularly between values of 6 m and 14 m with the deepest areas towards the downstream boundary (Soetaert and Herman, 1995). 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, 1995), the major parts of which are situated in the downstream area. The model presented here comprises the stretch 20 of river between the upstream boundary at Rupelmonde (river km 0) and the downstream boundary at Vlissingen (river km 104). elled processes oxic mineralisation, denitrification, nitrification, and primary production (for details see Hofmann et al., 2008b). Furthermore air-water exchange of carbon dioxide and oxygen as well as advective-dispersive transport of all chemical species are included. Acid-base equilibria as given in Table 1 have been considered for the pH calculation. Note that the dissociation constants (K * ) of the acid-base reactions are 5 calculated dynamically as functions of salinity, temperature and hydrostatic pressure, where salinity and temperature vary over time while the mean estuarine depth and thus the hydrostatic pressure remains constant over time. Furthermore all dissociation constants are converted to the free pH scale (Dickson, 1984). Organic matter has been split into a reactive (FastOM) and a refractory (SlowOM) fraction, entailing two different rates for the two fractions for oxic mineralisation and denitrification. The resulting mass balances for the state variables of the model are given in Table 2. Note that X signifies the set of all total quantities except for total alkalinity in the model (total quantities are called equilibrium invariants in Hofmann et al., 2008a,b).

The implicit pH modelling approach
In Hofmann et al. (2008b) the pH is modelled implicitly by numerically solving a system of equations constructed from the equilibrium mass action laws of the acid-base reactions given in Table 1 and the concentrations of the total quantities in X at every time step of the numerical integration of the equations given in Table 2. This implicit 20 pH modelling approach (operator splitting approach, Hofmann et al., 2008a), is equivalent to the approach presented by Follows et al. (2006) and similar to the approaches presented by Luff et al. (2001). Furthermore, it is inspired by classical pH calculation methods as given by Ben-Yaakov (1970) and Culberson (1980) and variations of it are used by Regnier et al. (1997) and Vanderborght et al. (2002). Due to its implicit nature, 25 this approach does not allow for quantifying the individual influences of the kinetically modelled processes.
202 2.4 The explicit pH modelling approach With their direct substitution approach, Hofmann et al. (2008a) present a new methodology for pH modelling that describes the pH evolution over time with an explicit expression for the rate of change of the proton concentration. Since all the kinetically modelled processes are independent from one another, they separately contribute to 5 the rate of change of the proton concentration d t i expresses the contribution of process i to the rate of change of the proton concentration d t into terms due to the kinetically modelled processes provides a quantitification of their influences on the pH.  (Dickson, 1981) is used 3 Assuming constant acid-base dissociation constants entails 15 which allows formulating a total derivative of total alkalinity  Table (2) into Eq. (5) and rearranging the terms, we arrive at an equivalent to Eq. (1) for the given model. This allows us to individually quantify the influence of oxic mineralisation, denitrification, 5 nitrification, primary production, air-water exchange and advective-dispersive transport on the pH if the acid-base dissociation constants are assumed to be constant over time. (Note again that a spatial gradient in the dissociation constants which is constant over time does not pose a problem.) In the following we describe how to apply the explicit pH modelling approach to a 10 system with time variable acid-base dissociation constants.

The explicit pH modelling approach with time variable dissociation constants
Letting the dissociation constants vary over time entails in X and the dissociation constants in K * . Obviously, the dissociation constants are functions of temperature T , salinity S and pressure P K * i = f i (T, S, P ) Since the mean depth in the model does not vary over time, we consider constant pressure P . However, the functions for temperature and salinity dependence of some 20 dissociation constants are expressed on the seawater pH scale (K * ,SWS ) or the total pH 204 scale (K * ,tot ) (Dickson, 1984) and not on the free pH scale (K * ,free ) which is consistently used in the model presented here. These dissociation constants were converted to the free pH scale, without loss of generality from the seawater scale by (Dickson, 1984;Zeebe and Wolf-Gladrow, 2001) This means the total derivative of [TA] considering dissociation constants variable over time can be written as Appendix A details how the partial derivatives of [TA] and of the dissociation constants can be calculated 4 . 4 Note that Eq. (10) contains partial derivatives of [TA] and K * i with respect to one of their variables. This entails that all other variables of these quantities, as defined by Eqs. (6) and (9)

205
In the same way as done in Hofmann et al. (2008a), we can derive a rate of change of the proton concentration from Eq. (10) which can be partitioned into contributions by the different kinetically modelled pro-5 cesses and by the influences of changes in the dissociation constants due to changes in their four variables T , S, , and [ HF]. This can be done by plugging in expressions for d [TA] d t , and d t as given in Table 2 and rearranging the terms. The result is an equivalent to Eq. (1) for the given model that takes into account the respective contributions of transport, air-water exchange of CO 2 , oxic mineralisation, denitrifica-10 tion, nitrification, primary production, the temperature and the salinity effect on the dissociation constants, as well as two terms for pH scale conversions 2.6 Implementation The model including the implicit and explicit pH modelling methods (Sect. 2.5) has been coded in FORTRAN within the ecological modelling framework FEMME (Soetaert et al., over the years 2001 to 2004. As mentioned above, the changes in the boundary conditions, the temperature forcing and the freshwater discharge (Table 3) are responsible for trends in the model results. Due to the minimal change in the mean estuarine temperature, the effect of changes in the temperature forcing has been neglected. In our simulation runs (as described above), boundary conditions and freshwater discharge 20 vary simultaneously, obscuring the effect of boundary conditions and the effect of freshwater flow change for single chemical compounds. Therefore, we executed a number of explorative runs in which freshwater discharge or boundary values (upstream and downstream) for individual state variables or groups of them were allowed to vary while freshwater discharge and boundary conditions for all other state variables remained 25 at 2001 values (Table 4). This has been done to investigate their individual effect on 208 the annual whole estuarine mean pH. End of the year 2001 conditions were used as initial conditions for these explorative runs and as a result the mean pH value for 2001 slightly differed from the one obtained from the simulation runs described above. describing the variations in the dissociation constants due to changes in S and T but without the pH scale conversion related terms; in the lower row considering all terms as described in Sect.  Influences of changes in the dissociation constants are about three orders of magnitude smaller than the influences of kinetically modelled processes. Furthermore their patterns along the estuary (not shown) depend on the respective implementation of the 10 model (e.g. on which pH scale the dissociation constants are calculated and to which pH scale they are converted) and are rather erratic and of limited scientific value: they are therefore not presented in Figs. 4 and 5. Yet, incorporation of these influences is necessary to obtain the excellent agreement between the explicitly calculated pH and the implicitly calculated pH as shown in Fig. 3. 15 Figure 4a exhibits a trumpet-like shape due to pronounced activity in the upper estuary, i.e. between river km 0 and 60. In this stretch of the estuary, the absolute influences of most kinetically modelled processes decline to stay at low levels until the mouth of the estuary. The most important proton producer at the upstream boundary is nitrification and its relative importance drops from 77% upstream to 11% downstream.

20
The proton production of oxic mineralisation also decreases from upstream to downstream. However, its relative importance as a proton producer increases from 23% at the upstream boundary to 64% at the downstream boundary. The most important proton consuming process is CO 2 degassing and its relative importance first increases from 50% at the upstream boundary to 92% at km 32 and then decreases again to 25 65% at the downstream boundary. Compared to CO 2 degassing, the proton consumption by primary production is rather small. It shows a steady downstream decrease with local maxima in the zone of maximal volumetric primary production in the estu-ary around km 48 and around km 67. The relative importance of primary production as a proton consumer increases from 4% at the upstream boundary to 38% at km 67 and decreases again to 33% at the downstream boundary. Denitrification is a proton consuming process with relatively low importance in the Scheldt estuary. Its relative importance is 2% at the upstream boundary, 1% at river km 32 and 0% along the 5 rest of the estuary. Advective-dispersive transport counteracts the dominant proton consuming or producing processes, exporting protons from the model boxes between the upstream boundary and around km 32 and importing protons from km 32 on until the downstream boundary. It shows a maximum of proton import around river km 48 and a secondary maximum around river km 67. At the upstream boundary advective-10 dispersive transport accounts for 44% of proton consumption while at river kilometres 48 and 67 it delivers about 50% of the protons. Figure 4b shows longitudinal profiles of volume integrated ("per river kilometre") influences on the proton concentration as calculated with Eqs. (13) to (22), averaged over the four modelled years. Table 6 shows selected values of those volume integrated 15 influences on the proton concentration.
As the estuarine cross section area increases from around 4000 m 2 upstream to around 76 000 m 2 downstream, while the mean estuarine depth remains at around 10 m, there is a much larger estuarine volume in downstream model boxes than there is in upstream model boxes. As a consequence, volume integrated proton production 20 or consumption rates in Fig. 4b are similar in the upstream and downstream region of the estuary. This is in contrast to volumetric rates which are much larger upstream than they are downstream for all processes. The mid-region of the estuary (between kms 30 and 60) can be identified as the most important region for volume integrated proton turnover. The volume integrated proton turnover of oxic mineralisation, primary 25 production and CO 2 degassing, is clearly larger downstream than it is upstream, while the volume integrated proton turnover of nitrification is still larger upstream than downstream. 211 Figure 5 shows a budget of proton production and consumption over the whole model area and one year, averaged over the four modelled years. It can be seen that CO 2 degassing and primary production are (except for the minor contribution of denitrification) the only processes that net consume protons in the estuary. Advective-dispersive transport, oxic mineralisation, and nitrification all net produce protons. CO 2 degassing 5 has the largest influence on the pH by causing the largest proton consumption, while nitrification is the main proton producer, closely followed by oxic mineralisation.
3.3 Factors responsible for the change in the mean estuarine pH from 2001 to 2004 Figure 6 shows the trend in the overal volume averaged pH in the estuary and the associated influences of the major kinetically modelled processes on the proton con- We use model scenarios to investigate the sensitivity of the estuarine pH to changes in freshwater flow and boundary conditions. Fig. 7 shows the results of the different 20 model scenarios summarised in Table 4.
Around 59% of the pH change in the system from 2001 to 2004 can be attributed to the change in the freshwater discharge (Fig. 7a) which also reproduces the general trend of decreasing absolute influences on the proton concentration. Especially the decline in the influence of nitrification can be clearly seen. However, the steep decrease in the influence of transport from 2001 to 2002 and its slight increase from 2003 to 2004 is not reproduced. As shown in Fig. 7h, about 44% of the pH change in the system from 2001 to 2004 can be attributed to the change in boundary conditions 212 (Note that in this complex non-linear model the pH changes due to separate freshwater discharge and boundary condition changes are not necessarily additive). In the more erratic pattern of influences displayed in Fig. 7h, one can identify the steep decrease in the influence of transport from 2001 to 2002, as well as its increase between the two following years. 5 Influences via [ CO 2 ] and [TA] are most important for changes in the pH, as changes in freshwater flow for these quantities account for 49% (Fig. 7b) and changes in boundary conditions account for 28% (Fig. 7i) of the total pH change to the system. However, the pattern of influences in Fig. 6 cannot fully be explained by just influences via [ CO 2 ] and [TA]. Especially the decrease in the influence of nitrification and the 10 distinct pattern of the influence of advective-dispersive transport is missing.
Influences via [ NH + 4 ] are also substantial: the change of freshwater flow accounts for 22% (Fig. 7d) and the change of the boundary conditions accounts for 19% (Fig. 7k) of the total pH change to the system. It can be seen that the influence is indirect as the influence of nitrification is decreased. Influences via [ NH Freshwater flow changes for S do not lead to a pH increase but a decrease of 22% primarily via changes in the influence of transport (Fig. 7c).

Quantification of proton production and consumption along the estuary
To our knowledge, we are the first to quantify the influences of kinetic processes on the pH for an entire estuarine ecosystem like the Scheldt estuary. Although pH profiles have been simulated quite accurately (e.g. Regnier et al., 1997;Vanderborght 5 et al., 2002Vanderborght 5 et al., , 2007Blackford and Gilbert, 2007;Hofmann et al., 2008b), the attribution of pH changes to specific biogeochemical processes has only been done qualitatively. Neglecting minor contributions like primary production and denitrification and the attenuation of gradients due to advective-dispersive transport, the pH profile in the Scheldt estuary is mainly the result of a balance between two biogeochemical reactions, ni-10 trification and oxic mineralisation, which produce protons and CO 2 degassing which consumes protons. This is fully consistent with and supports the findings by Regnier et al. (1997), but our treatment yields more and quantitative information.
The budget of volumetric influences on the proton concentration ( Fig. 4a; Table 5) exhibits the same trumpet like shape (higher values in the upstream region than in 15 the downstream region) as budgets for total ammonium, dissolved inorganic carbon, oxygen, and nitrate (Hofmann et al., 2008b). This confirms that high proton turnover is associated to high activity in (kinetically modelled) biogeochemical processes.
Moreover, the budgets for proton production and consumption presented here (Fig. 4) are approximately mirror images of the budgets for oxygen sources and sinks (Hofmann 20 et al., 2008b). That suggests that proton production and consumption are inversely correlated with oxygen production or consumption. The underlying cause is that oxygen consumption reactions transfer electrons to the oxygen atoms producing reduced oxygen (for example in nitrate or in water). The chemical species that is oxidised (for example the nitrogen in ammonia upon nitrification) is electron-rich before the reaction and electron-depleted and bound to oxygen after the reaction. An electron-rich species, however is more prone to bind electron-depleted protons than an electron poor species. Thus there is a general trend that upon oxidation of a chemical species 214 the electron-depleted protons are produced, and this provides a direct link between the oxygen and the proton budgets in our model. While CO 2 degassing accounts for the largest total proton turnover per year in the whole estuary, this process acts as a "buffer" for the effects of other processes on the proton concentration, since its magnitude is very sensitive to the current pH of 5 the system. The same holds for the influence of advective-dispersive transport with its buffering character being so pronounced that it changes sign along the estuary. This entails that, given a certain freshwater flow and certain boundary conditions, nitrification and to a lesser extent oxic mineralisation and primary production are the prime factors influencing the pH profile along the estuary, while CO 2 degassing and 10 advective-dispersive transport counteract their effects. This is consistent with findings of (Vanderborght et al., 2002) who also identify nitrification in the Scheldt as a process influencing CO 2 degassing via the pH.
The influences of changes in the dissociation constants are several orders of magnitude smaller than the influences of kinetic processes on the proton concentration 15 (Tables 5 and 6). Therefore, when describing the factors that govern the order of magnitude of the proton concentration of a system, they can be neglected. However, to describe the pH accurately, i.e. more accurate than 0.1 pH units, they should be included. This is especially important for modelling the proton concentration explicitly over a longer period of time, since deviations in

Synopsis
The main factors governing the pH in the heterotrophic, turbid, tidal Scheldt estuary 15 can be summarised as given in Fig. 8. Within the estuary, i.e. with given boundary conditions and freshwater flow, the dependencies depicted with red arrows govern the pH: mainly nitrification and oxic mineralisation (both producing protons) and primary production (consuming protons) influence the proton concentration, an effect which is "buffered" by the effect of CO 2 degassing and advective-dispersive transport. Consid-20 ering changes of the mean pH in the estuary over the years 2001 to 2004 the dependencies depicted with blue arrows are the governing factors. Changes in boundary conditions and freshwater flow mainly influence [ CO 2 ] and [TA] which can be considered a "direct" effect on the proton concentration. However, changes in boundary conditions and freshwater flow also change [ NH + 4 ] which in turn influences the ef-25 fect of nitrification, an "indirect" effect on the proton concentration. This "indirect" effect of changes in boundary conditions and freshwater flow is significant as it amounts to about 50% of the "direct" effect. 216 1. A method to quantify the influences of kinetically modelled processes on the pH of a system with time variable acid-base dissociation constants was presented and verified against an existing pH modelling approach.
2. By applying this method to a model of the Scheldt estuary we have identified nitri-5 fication as the main process governing the pH profile along the estuary while CO 2 degassing and advective-dispersive transport (given a certain freshwater flow and certain boundary conditions) "buffer" its effect. However, CO 2 degassing accounts for the largest total proton turnover per year in the whole estuary.
3. A clear inverse correlation between oxygen and proton turnover was found, con-10 sistent with theoretical considerations of redox chemistry.
4. While the influences of changes in the dissociation constants might be neglected in approximate whole estuarine budgets, they are important to correctly model the proton concentration explicitly in systems where the acid-base dissociation constants are assumed to be variable over time. and also -to a significant amount -"indirectly" via [ NH + 4 ] and the nitrification rates in the estuary. ∂[X j ] ) have been calculated analytically: Note that this list is system specific, e.g. since H 2 SO 4 dissociation is not considered an acid-base reaction in our system, HSO − 4 is considered a monoprotic acid.
A2 Partial derivatives with respect to and of the dissociation constants 5 In general, the terms K * 1 and K * 2 can be calculated from temperature (T , in Kelvin) and salinity (S), following e.g. Roy et al. (1993), in the form which allows to write Without loss of generality, we assume K * 1 and K * 2 to be calculated on the seawater pH scale and then converted to the free proton scale. According to Dickson (1984) and Zeebe and Wolf-Gladrow (2001)    Jourabchi, P., Van Cappellen, P., and Regnier, P.: Quantitative interpretation of pH distributions in aquatic sediments: A reaction-transport modeling approach, Am. J. Sci., 305, 919-956, available at: GotoISI ://000236181100002, 2005. 200 Luff, R., Haeckel, M., and Wallmann, K.: Robust and fast FORTRAN and MATLAB (R) libraries to calculate pH distributions in marine systems, Comput. Geosci., 27, 157-169, available at:  OxFastOM and R OxSlowOM are the reaction rates of oxic mineralisation for the reactive and refractory organic matter fraction respectively. Similarly, R DenFastOM , R DenSlowtOM , R Nit , and R PP are the rates of denitrification, nitrification and primary production. E C and Tr C express the air-water exchange and advective-dispersive transport rates of the respective chemical species. p P P NH + 4 is the fraction of NH + 4 usage of primary production as explained in Hofmann et al. (2008b).      −7.65 10 −4 (0%) −6.94 10 −3 (0%) −1.60 10 −2 (0%) −1.61 10 −2 (0%) −1.60 10 −2 (0%) 6.83 10 −3 (0%) pH profiles along the Scheldt estuary salinity gradient. The blue line represents the pH calculated with a closed system model (comparable to Mook and Koene, 1975); the red line represents the pH calculated with an open system model (comparable to Whitfield and Turner (1986) but with realistic kinetic CO 2 air-water exchange instead of a fully equilibrated system); the black line represents the pH calculated with the full biogeochemical model as presented in Hofmann et al. (2008b).