Biogeosciences pH modelling in aquatic systems with time-variable acid-base dissociation constants applied to the turbid , tidal Scheldt estuary

A new pH modelling approach is presented that explicitly quantifies the influence of biogeochemical processes on proton cycling and pH in an aquatic ecosystem, and which accounts for time variable acid-base dissociation constants. As a case study, the method is applied to investigate proton cycling and long-term pH trends in the Scheldt estuary (SW Netherlands, N Belgium). This analysis identifies the dominant biogeochemical processes involved in proton cycling in this heterotrophic, turbid estuary. Furthermore, information on the factors controlling the longitudinal pH profile along the estuary as well as long-term pH changes are obtained. Proton production by nitrification is identified as the principal biological process governing the pH. Its acidifying effect is mainly counteracted by proton consumption due to CO2 degassing. Overall, CO 2 degassing generates the largest proton turnover in the whole estuary on a yearly basis. The main driver of long-term changes in the mean estuarine pH over the period 2001 to 2004 is the decreasing freshwater flow, which influences the pH directly via a decreasing supply of dissolved inorganic carbon and alkalinity, and also indirectly, via decreasing ammonia loadings and lower nitrification rates.


Introduction
pH is often considered a master variable for the chemical state of aquatic ecosystems, since almost any process affects pH either directly or indirectly (e.g.Stumm and Mor-Correspondence to: A. F. Hofmann (a.hofmann@nioo.knaw.nl)gan, 1996; Morel and Hering, 1993).But as a result of this complex interplay of different processes, we only have limited understanding of the factors controlling the pH of natural waters (the water column of rivers, lakes, and oceans, as well as the pore water of sediments).Given the ongoing acidification of the surface ocean (e.g.Orr et al., 2005) and coastal seas (e.g.Blackford and Gilbert, 2007), and the large impact that future pH changes may have on biogeochemical processes and organisms (e.g.Gazeau et al., 2007;Guinotte and Fabry, 2008), it is desirable to obtain a better quantitative understanding of factors that regulate pH in natural aquatic systems.This requires modelling tools, which basically satisfy two criteria: (1) models should be able to accurately predict the time evolution of pH under certain biogeochemical scenarios, and (2) models should be able to assess the contribution of individual biogeochemical processes in observed pH changes.
The first aspect has received quite some attention over the last two decades.A number of modelling approaches have been presented that allow to predict the time evolution of pH in aquatic ecosystems (e.g.Boudreau and Canfield, 1988;Regnier et al., 1997;Vanderborght et al., 2002;Hofmann et al., 2008b).The second aspect, however, has received far less attention.Improving and extending previous work by Jourabchi et al. (2005) and Soetaert et al. (2007), we have recently presented an extension of the conventional pH modelling approach, developing a procedure that allows to quantify how different biogeochemical processes contribute to the overall proton cycling in an aquatic ecosystem (Hofmann et al., 2008a).This method basically splits the total rate of change of protons into individual contributions of the processes that are driving the pH change (details are given below).However, our original presentation contained an important restriction: acid-base dissociation constants were assumed to remain constant in time.Although this assumption is reasonable in some aquatic systems, it is not valid in environments that experience substantial temporal changes in temperature and salinity.
One such environment where temporal variability in salinity and temperature may strongly affect pH are estuaries.Estuarine systems are characterized by pronounced salinity gradients, marked seasonality, a large exchange of CO 2 with the atmosphere, and intensive biogeochemical processing (Soetaert et al., 2006;Regnier et al., 1997;Vanderborght et al., 2002).This makes them suitable (though demanding) testbeds for pH modelling methods.Hence, this study has two main objectives: (1) the extension of the 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, and (2) the validation of this approach for an estuarine system (the Scheldt estuary).We will compare calculated pH values with available data, quantify the proton production or consumption due to different processes (advective-dispersive transport, CO 2 air-water exchange, primary production, nitrification etc.), and investigate the factors that are responsible for the observed trend in the mean estuarine pH over the years 2001 to 2004.

The Scheldt estuary
The turbid, tidal Scheldt estuary forms the end part of the Scheldt river (Fig. 1), which is about 350 km long in total and drains a basin of 21 500 km 2 in the northwest of France, the northern part of Belgium and the southwest of the Netherlands (Soetaert et al., 2006).The model domain here includes the final stretch of river, located between the upstream boundary at Rupelmonde (river km 0) and the downstream boundary at Vlissingen (river km 104).Water movement in the Scheldt estuary is dominated by large tidal displacements.The volume of water that enters during flood tide at the downstream boundary of the estuary is 200 times larger than the freshwater discharge integrated over a whole tidal cycle at the upstream boundary (Vanderborght et al., 2007).Yearly averaged, the freshwater discharge is around 100 m 3 s −1 (Heip, 1988).The cross sectional area of the estuarine channel shows a regular trumpet-like shape, gradually widening from around 4000 m 2 upstream to around 75 000 m 2 downstream (Soetaert et al., 2006).The mean water depth varies irregularly between values of 6 m and 14 m with the deepest areas located towards the downstream boundary (Soetaert and Herman, 1995).The estuary has a total tidally averaged volume of 3.62 km 3 and a total tidally averaged surface area of 338 km 2 (Soetaert et al., 2006;Soetaert and Herman, 1995).
Table 1.Biogeochemical processes included in the model formulation.Process rates: R Ox =oxic mineralisation, R Den =denitrification, R Nit =nitrification, R PP =primary production.γ denotes the C/N ratio of organic matter: γ =4 for the reactive organic matter fraction, γ =12 for the refractory fraction.p N denotes the fraction of NH + 4 in the total DIN uptake of primary production.
2.2 Biogeochemical model formulation Hofmann et al. (2008b) provided a detailed model assessment of the carbon and nitrogen cycling in the Scheldt estuary.Their model and results provide a solid basis and are used as a starting point.The one-dimensional reactive transport model here incorporates the same biogeochemical reaction set, and uses the same forcings and kinetic descriptions of process rates as in Hofmann et al. (2008b).
The biogeochemical reaction set accounts for oxic mineralisation, denitrification, nitrification, and primary production with stoichiometries as given in Table 1 (for details see Hofmann et al., 2008b).Furthermore, the model includes two physical transport processes: the exchange of carbon dioxide and oxygen across the air-water interface and advectivedispersive transport of dissolved constituents.The model tracks the concentration of twelve chemical species in total.Organic matter is split into a reactive (FastOM) and a refractory (SlowOM) fraction, thus leading to two separate rates for oxic mineralisation and denitrification.The resulting mass balance equations for all state variables are given in Table 2.More detailed information on the model formulation can be found in Hofmann et al. (2008b).We deliberately did not change the design and parameterization of the model presented in Hofmann et al. (2008b) to show that the pH modelling procedure presented here can be implemented as an extension to existing biogeochemical models.The only, but important, difference between the two model implementations is the pH calculation procedure (as explained below).
Seven acid-base equilibria are accounted for in the pH calculations (Table 3).The respective dissociation constants (K * 's) are calculated as functions of salinity (S), temperature (T ) and hydrostatic pressure (P ).Salinity and temperature may vary over time, while the mean estuarine depth, and so the hydrostatic pressure, remains constant over time.The (S,T ,P ) relations that have been reported for the dissociation constants in the literature differ in their associated pH scale (Dickson, 1984).Before being used in the model calculations, all dissociation constants are systematically converted to the free pH scale.

Data set and model parameterization
Data on state variables, model parameters and process rates were collected from various sources and spanned the four year period [2001][2002][2003][2004].Monthly values of the discharge at the upstream boundary were obtained from the Ministry of the Flemish Community (MVG).Percental flow increases along the estuary were obtained from the SAWES model (van Gils et al., 1993;Holland, 1991) and are used to calculate the lateral input along the estuary as percentages of the upstream flow.These percentages have been assumed constant in time.The average water depths D of the 13 boxes in the Scheldt model MOSES from Soetaert and Herman (1995) were interpolated to the 100 model boxes used here.The cross sectional area along the estuary was specified as a continuous function of the distance along the estuary (river kilometres) as given in Soetaert et al. (2006).
Data for , organic matter, pH (on the NBS scale), temperature and salinity were obtained for 16 stations between the upstream boundary at Rupelmonde (Belgium) and the downstream boundary at Breskens/Vlissingen (The Netherlands) by monthly cruises of the RV "Luctor" from the Netherlands Institute of Ecology (NIOO).
The model requires suitable boundary values for each state variable at the upstream and downstream boundary.
, temperature and salinity, data collected at the monitoring station close to Breskens were taken as downstream boundary conditions.Similar data collected at the monitoring station close to Rupelmonde were used as upstream boundary conditions.The boundary values for [ HSO − 4 ], [ HF] and [ B(OH) 3 ] were estimated from salinity at the boundaries assuming standard seawater proportions.
One problem was that suitable boundary value data were available for pH, but not for dissolved inorganic carbon ([ CO 2 ]) or total alkalinity ([TA]).However, [TA] data were obtained for the year 2003 for various stations within the estuary from Frederic Gazeau (personal communication and Gazeau et al., 2005).We used these [TA]   are the reaction rates of oxic mineralisation for the reactive and refractory organic matter fraction respectively.Similarly, R Den FastOM , R Den SlowOM , 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.[TA] denotes total alkalinity.Roy et al. (1993), with extension for whole salinity range as given in Millero (1995); K * W : Millero (1995); Millero (1995)).All stoichiometric equilibrium constants were converted to the free proton scale, and pressure corrected according to Millero (1995) using the mean estuarine depth for each model box.Subsequently, they are converted to volumetric units using the temperature and salinity dependent formulation for seawater density according to Millero and Poisson (1981).
available [TA] data, and subsequently, the [ CO 2 ] values were adjusted to provide the best fit.These [ CO 2 ] values calibrated for the year 2003 were subsequently used for all other years.To be consistent with pH boundary conditions available for all modelled years, [TA] boundary conditions were not directly calibrated from available 2003 [TA] measurements within the estuary, but calculated from calibrated [ CO 2 ] values and available pH boundary conditions according to the total alkalinity definition used in our model (see below).

pH modelling approaches
In a recent review paper on pH modelling, we showed that there are two main approaches to calculate pH in a reactive transport model of an aquatic ecosystem (Hofmann et al., 2008a).These two methods can be referred to as implicit and explicit.The implicit method is the conventional pH modelling approach.This paper advances a new (extended) version of the explicit approach.To verify their consistency, we implemented both approaches side-by-side.

The implicit pH modelling approach
The implicit approach is the traditional method, which goes back to the work of Ben-Yaakov (1970) and Culberson (1980), and also forms the standard approach included in textbooks on aquatic geochemistry (Millero, 2006;Sarmiento and Gruber, 2006).The more recent treatments of the pH modelling problem by Luff et al. (2001) and Follows et al. (2006) are also based on this implicit method.Furthermore, within the context of estuarine modelling, the implicit approach has been used by Regnier et al. (1997) and Vanderborght et al. (2002), and we also implemented it in our previous paper on the biogeochemistry of the Scheldt estuary (Hofmann et al., 2008b).The implicit approach is termed "implicit" because it essentially takes a detour via alkalinity to arrive at pH.In other words, one first expresses how biogeochemical processes affect the alkalinity mass balance (last equation in Table 2), and subsequently, at every time step, one calculates the pH from alkalinity and other total quantities (i.e. total carbonate, total borate, etc.).The latter step comes down to the solution of a non-linear equation system governed by the mass action laws of the acid-base reactions (as given in Table 3).Because of this two-step procedure, there is no direct link between biogeochemical process rates and the actual proton production or consumption these processes generate.As discussed extensively in Hofmann et al. (2008a), this is a clear drawback of the implicit approach, since it does not allow the modeller to quantify how important an individual process really is in the total proton balance of the system.

The explicit pH modelling approach
Extending and improving the work of Jourabchi et al. (2005) and Soetaert et al. (2007), we recently presented a new approach to pH calculations in reactive transport models (Hofmann et al., 2008a).This method, called explicit here, allows us to describe the pH evolution over time (just like the implicit method), but in addition, it derives the contribution of separate biogeochemical processes to the overall rate of change of the proton concentration.First, we briefly review the main features of the explicit pH method as it was presented in (Hofmann et al., 2008a), where it was assumed that dissociation constants remain constant in time.Subsequently, we discuss how the method can be generalized to time-variable dissociation constants, which is the novel aspect here.
The central idea of the explicit pH modelling method is that the total rate of change of the proton concentration can be written as where dt i expresses the contribution of i-th process (e.g.nitrification, CO 2 degassing, etc.).This partitioning of the total rate of change of the proton concentration into separate contributions by different processes thus quantifies how important a given process is in the overall production and consumption of protons.
The definition of total alkalinity [TA] (Dickson, 1981) that is associated with the acid-base reaction set (Table 3) in our biogeochemical model is If one assumes that the acid-base reactions are in equilibrium, then each of the dissociated species on the right hand side can be expressed as a function of the proton concentration [H + ], the concentrations of the total quantities and some associated dissociation constants K * i .
As a consequence, total alkalinity can be written as If one assumes that the dissociation constants K * i do not depend on time, the time derivative of the total alkalinity can be specified as This expression can be directly rearranged so that it provides an expression for the total rate of change of protons By plugging the expressions for d [TA]  dt and dt given in Table (2) into Eq.( 7) and regrouping terms, we arrive at where the terms on the right hand side are given in Table 4 (Eqs.15 to 20).This expression matches our sought-after expression Eq. (1).If the acid-base dissociation constants can be assumed constant over time, this expression allows us to individually quantify the contributions of oxic mineralisation, denitrification, nitrification, primary production, airwater exchange and advective-dispersive transport to the rate of change of the proton concentration.Note that the acidbase dissociation constants should only be independent of time (accordingly, a stationary spatial gradient in these constants does not pose a problem).In the following section, we introduce the novel aspect of this publication: we describe how to apply the explicit pH modelling approach to a system with time variable acid-base dissociation constants.

The explicit pH modelling approach with time variable dissociation constants
For a number of problems in aquatic biogeochemistry, the assumption of time-invariant acid-base dissociation constants is not tenable.Estuaries form a prime example of this.Because salinity and temperature show a marked seasonal dependence, the acid-base dissociation constants will also change over time.When accounting for this timedependency, one should account for the direct effect of temperature T , salinity S and pressure P on the dissociation constants By assuming a constant pressure, one can remove P from the list of variables, as we will do here.However, some temperature and salinity dependencies of dissociation constants are only available on other pH scales (without loss of generality the seawater pH scale (K * ,SWS )) and not on the free pH scale (K * ,free ) as used here.These dissociation constants can be converted to the free pH scale using the relation (Dickson, 1984;Zeebe and Wolf-Gladrow, 2001) This implies that, in general, the dissociation constants are also functions of [ HSO − 4 ] and [ HF] Note that if one assumes that [ HSO − 4 ] and [ HF] are simply proportional to S, [ HSO − 4 ] and [ HF] need not to be treated as independent variables.However, to preserve the generic nature of the presented method, we treat S, [ HSO − 4 ], and [ HF] as independent variables in the derivations below.Accordingly, if one assumes that T , S, , and [ HF] are dependent on time, and so the dissociation constants K * i become dependent on time, the time derivative of total alkalinity (Eq.5) now provides the more elaborate expression Appendix A details how the partial derivatives in this expression can be explicitly calculated as a function of the proton concentration, the total quantities X j and the dissociation constants K * i .Again, Eq. ( 12) can be directly transformed into following expression for the total rate of change of protons In a similar way as above, one can plug the expressions for dt , and dt as given in Table 2 into Eq.( 13) and rearrange the terms so that terms associated with a given process are grouped together.This leads to with process specific terms as given in Table 4.In addition to the processes already featuring in Eq. ( 8), there are now additional contributions related to the effect of changing temperature and salinity (Eqs.21 and 22), as well as two terms linked to pH scale conversion (Eqs.23 and 24).

Model implementation
The biogeochemical reactive transport model and the pH calculation procedures were implemented in FORTRAN within the ecological modelling framework FEMME (Soetaert et al., 2002).The model code is an extension of the code discussed in (Hofmann et al., 2008b).As noted above, the implicit pH calculation routine was already present in this code and has been retained.Here, the code was extended with the explicit pH modelling method as described in Sect.2.4.3.Accordingly, the model code implements two different pH calculation procedures, which independently predict the pH evolution of the estuary, thus allowing for a consistency check between these two methods.The model code can be obtained from the corresponding author or from the FEMME website: http://www.nioo.knaw.nl/projects/femme/.Post processing of the model results and visualization of the model output was done using the programming language R (R Development Core Team, 2005).

Baseline simulation, calibration, generation of output
The baseline simulation is the model solution that best fits the data.It is used as the starting point for the subsequent sensitivity analysis.To arrive at the baseline simulation, a time dependent simulation over a four year period was performed spanning the years from 2001 to 2004.Boundary conditions (values for [TA], S, , and [ HF], temperature, and the freshwater flow) were imposed onto the model with monthly resolution, using the data as described in Sect.2.2.The initial conditions for this time dependent simulation were obtained by a steady state model run, where all boundary conditions were fixed at January 2001 values.The state variable [H + ] (for the explicit pH modelling approach) has been initialized with an implicitly calculated initial pH.The model was numerically integrated over time with an Euler scheme using a time-step of 0.00781 days.To be comparable to NIOO monitoring pH data (measured on the NBS pH scale), the simulated pH values (calculated on the free pH scale) were converted to the NBS scale using the Davies equation (cf.discussion of the use of the Davies equation in Hofmann et al., 2008b).
A sequential calibration-validation procedure was employed using field data for total ammonium, nitrate + nitrite, oxygen, organic matter, salinity, pH, total alkalinity and nitrification rates.First, the data for the year 2003 were used to calibrate the process rate parameters of the model in a run that covered the year 2003 only.Subsequently, data for the remaining years (2001, 2002, and 2004) were used to independently validate the model in the above mentioned fouryear run (keeping rate parameters at values calibrated for 2003, but changing boundary conditions and physical forcings to values of the respective years).In addition, nitrification rates for the year 2003 were taken from Andersson et al. (2006) and were used to independently validate the model.Additional details on data sources, analytical methods and model calibration, as well as fits of model output to data other than pH can be found in Hofmann et al. (2008b).
To arrive at yearly-averaged, longitudinal profiles, the model output (concentrations, process rates, pH) of the timedependent simulation for every model box was suitably averaged over the desired period.A similar procedure was applied to the monthly monitoring data.We also created volume integrated, "per model box" longitudinal profiles.To this end, the volumetric rates as calculated with Eqs. ( 13) to ( 22) were multiplied with the volume of the respective model box along the river.Subsequently, we then integrated longitudinally to arrive at "whole estuarine" process rates, i.e. we multiplied the rates for each model box with the volume of the respective box and summed up.We also averaged longitudinally to arrive at "mean eastuarine" values of concentrations, pH, and process rates.To this end, the value in each model box was weighted with the volume of that box, summed up and divided by the total estuarine volume.Additionally, all obtained profiles, mean and total values have been averaged over the four modelled years.
Figure 2 shows the resulting yearly-averaged pH profiles of the baseline simulation.The final pH fit shows a slight underestimation of pH in the upstream region and a slight overestimation in the downstream region.This could most probably be remedied with a more elaborate model including additional biogeochemical processes and a more elaborate hydrology (migrating from a one-dimensional model to two or even three dimensions).However, to get an idea about the dynamic equilibrium that maintains the pH along the estuary, we consider the present model elaborate enough and the resulting fit good enough.

pH modelling verification runs
The implicit and explicit pH modelling approaches are two different methods for calculating the same quantity.This means that if these two methods are implemented correctly, they must produce exactly the same pH response (provided the underlying biogeochemical model is exactly the same, as is the case here).This was verified by executing the baseline simulation with both the implicit method and the explicit pH modelling method with time-variable dissociation constants as in Sect.2.4.3.These two methods provide indeed the same response within the numerical precision of the code, as shown in Fig. 3c for the year 2003 (the other 3 years show a similar response).Some details warant further discussion.The extended explicit pH modelling method presented here explicitly accounts for time-variable acid-base dissociation constants.This adds a number of d[H + ] dt i terms to the proton balance equation (compare Eq. ( 8) to Eq. ( 14)).It is useful to rewrite Eq. ( 14) in following form The P P and P C terms gather all proton producing and consuming processes given in Eq. ( 8).The stands for all terms that are present in Eq. ( 14) but not in Eq. ( 8).The term includes terms that account for two effects: timedependent changes in salinity and temperature, which influence the dissociation constants directly, and time-dependent changes in [ HSO − 4 ] and [ HF], which influence the dissociation constants via pH-scale conversion.As shown below, the proton production P P and the proton consumption P C are very large compared to the net rate of change of protons d[H + ] dt , and as a result, P P and P C almost balance each other.Additionally, the term is also small compared to P P and P C. Accordingly, one is tempted to omit from Eq. ( 14).Yet, we found that such an omission introduces substantial deviations as shown in Figs.3a,  b.Omitting all terms in from Eq. ( 14), yields pH profiles that are substantially different from the correctly calculated ones (Fig. 3a).Including the terms that account for the direct effect of temperature and salinity (but not including the pH scale effect) yields a better result (Fig. 3b).However, a small deviation is still present.Additionally including the pH scale conversion related terms makes the response of the explicit method identical to the the implicit one, as required (Fig. 3c).
To explain why the neglect of very small terms in (so small that they irrelevant when comparing P P and P C and the contributions of individual processes therein) can have relatively large consequences, one must realize that P P and P C are very large compared to the net rate of change of protons d[H + ] dt .Accordingly, a small disturbance of this balance (like ignoring ) can have a huge impact on d[H + ] dt , and hence, on the model predicted pH value.This explains the deviations in Fig. 3.In conclusion, the explicit pH modelling method is powerful, but one needs to ensure that it is consistently implemented and no terms are omitted.

Sensitivity analysis
To obtain an idea about what factors (estuarine mixing, salinity gradient, air-water gas exchange, and biological processes) control the longitudinal pH profile in the estuary, we performed a number of exploratory model simulations.sequentially activated to arrive at more biogeochemically complex scenarios.The basic scenario simply served to investigate the effect of estuarine mixing on the pH profile.To this end, we purposely neglected gas exchange (CO 2 and O 2 ) and biological processes, and assumed no dependence of the stoichiometric constants on salinity (values were calculated at T =15 • C and S=15).Accordingly, in this basic scenario, the pH profile along the estuary is solely determined by conservative mixing of [TA] and [ CO 2 ], driven by the difference in upstream and downstream boundary values (Fig. 4).
Starting from this basic scenario, we conducted four additional simulations where groups of processes where sequentially activated.In a first scenario, gas exchange and biological processes where still neglected ("closed system, no biology"), but the salinity dependence of the stoichiometric constants was now explicitly accounted for.In two subsequent simulations, either gas exchange was additionally activated ("open system, no biology") or biological processes were included ("closed system, biology").In a fourth and final simulation, all processes were activated thus leading to the baseline simulation as discussed above ("full biogeochemical model").To be comparable to the work of Mook and Koene (1975), the resulting pH profiles are plotted against the salinity gradient in the estuary (Fig. 5).  5 shows the inter-annual variations of these factors as well as variations in values of important state variables volume averaged over the whole estuary.We carried out a sensitivity analysis to find out which particular factors could explain the observed four-year pH trend.

Factors governing changes in the mean
Because the temperature follows a predictable seasonal cycle, inter-annual changes in the temperature forcing are negligible.Model runs with all parameters fixed at the 2001 values, but with the actual time-variable temperature forcing from 2001 to 2004, show no discernably different pH profile.This leaves to investigate the effect of changing freshwater input and changing boundary concentrations.We suspected that the biogeochemistry of a particular compound was predominantly influenced by changes in its total "load", i.e., the total input at the upstream boundary.For example, the ammonia load is simply defined as the freshwater discharge Q times the upper boundary concentration [ NH + 4 ] up .In our sensitivity analysis, we were particularly interested which loading changes (i.e. of which chemical variable) were responsible for the observed inter-annual pH changes.
In the baseline simulation described above, boundary conditions and freshwater discharge vary simultaneously over time, forced by the monthly monitoring data.To disentangle the effects of freshwater discharge and boundary changes on the total loading, we executed 14 simulations in which the freshwater discharge or the boundary concentration values (upstream and downstream) were independently varied.In these simulations, all other state variables remained at 2001 values.Table 6 lists the groups of state variables for which the loading has changed, either by changing the freshwater flow (left column) or the boundary conditions (right column).Each time the resulting pH change of the four year period was expressed as a fraction of the pH change arrived at in the baseline scenario to quantify the importance of the given parameter change in explaining the four-year pH trend.

Factors controlling the pH profile along the estuary
We performed a number of exploratory simulations to investigate the major controls on the longitudinal pH profile in the Scheldt estuary, which is also characteristic for other estuaries.A first striking aspect is that the pH increases from the upstream freshwater boundary to the downstream marine boundary.A "skeleton" simulation (with no gas exchange, no biogeochemistry and no dependence of acid-base constants on salinity) shows that this pH increase is simply the result of conservative mixing of [TA] and [ CO 2 ], as the ratio [ CO 2 ]/[TA] decreases (Fig. 4).
On top of the overall increase in pH, the observed pH profile shows a changing curvature with a distinct pH minimum at low salinities.Via a sensitivity analysis, we examined how strongly different groups of processes (salinity dependence of acid-base constants, gas exchange, biology) influence the shape of the longitudinal pH profile.Figure 5 shows that the simulation, which considers a closed system (no CO 2 and O 2 exchange with the atmosphere) and no biological processes, exhibits a distinct pH minimum at low salinities (blue line).However this simulated pH minimum is located at higher salinities (more downstream) than the observed pH minimum.Enabling gas exchange, while keeping biological processes inactivated, results in a concave pH profile with high pH values and no minimum (red line).Conversely, activating biological processes in a closed system results in a convex pH profile with a minimum but with pH values far below observations (magenta line).Finally, a full model simulation with all biological processes and gas exchange toggled on (black line) shows a profile with changing curvature and the distinct pH minimum at low salinities.This full model fits the observed data best, although some discrepancies remain (a more sharply recovering pH minimum in the data

Proton production and consumption along the estuary
The prime advantage of the explicit pH modelling method is that it calculates the proton production/consumption rates associated with each individual reactive transport process.Figure 6a shows the contributions of individual processes to the net rate of change of protons as calculated with Eqs. ( 13) to ( 22).Longitudinal profiles of proton production or consumption rates (per unit of solution volume) were extracted from the baseline simulation, averaged over the four year period (for every model box), and plotted cumulatively.Table 7 lists resulting values at selected positions (model boxes) in the river where the pH profiles in Fig. 6 show interesting features (these locations are also indicated in Fig. 1).
In a first step, we can look at the overall proton cycling, i.e., the total proton production (P P ) and total proton consumption (P C) along the estuary.A first observation is that the P P and P C terms are always four to five orders of magnitude larger than the actual rate of change of protons (which is in the 10 −5 mmol m −3 range).This implies that proton production and consumption are nearly balanced, and that the internal cycling of protons far outweighs the net proton change over time.A second aspect is that the proton production rates, which are linked to changes in the dissociation constants, are about three orders of magnitude smaller than those of other processes.Therefore, they are not presented in Figs. 6 and 7. Yet, as noted above, incorporation of these influences is necessary for a consistent implementation of the explicit pH approach (i.e. to model absolute pH values correctly).
As shown in Fig. 6a, the proton production and consumption per unit volume shows a marked decrease in the upper half of the estuary (between river km 0 and 60), after which the decrease proceeds more gradually.Over the whole estuary, P P and P C decrease by a factor of 20, from around 1.3 mmol m −3 yr −1 to 0.06 mmol m −3 yr −1 .This decrease in proton turn-over generates the trumpet-like shape in Fig. 6a, and can be attributed to a similar decreasing trend in the biogeochemical activity per volume of solution (as discussed below).
Along the estuary, the magnitude of the individual contributions d[H + ] dt i generally follows the decreasing trend in the total proton production/consumption.The exceptions are primary production, whose proton consumption remains relatively constant, and advective-dispersive transport, which shows a noticeable profile (further discussed below).However, there are some marked changes in the relative importance of processes in the overall proton cycling.Nitrificaton and oxic mineralisation produce protons; CO 2 degassing, primary production, and denitrification consume protons.The dominant proton producer at the upstream boundary is nitrification (77%), yet its relative importance drops to 11% downstream.In parallel, the relative importance of oxic mineralisation as a proton producer increases from 23% upstream to 64% at the downstream boundary.In terms of proton consumption, the most important process is CO 2 degassing.Its relative importance increases from 50% at the upstream boundary to 92% at km 32 and then decreases again to 65% at the downstream boundary.Compared to CO 2 degassing, the proton consumption by primary production is small.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.The proton consumption due to denitrification is not important in the Scheldt estuary (around 1% along the estuary).
The role of advective-dispersive transport in proton transport is markedly different from that of the biogeochemical reactions.Advective-dispersive transport counteracts the dominant proton consuming or producing processes.As a result of that, it switches sign.Around river km 32, it switches from proton consumption (importing protons into a model box) to proton production (exporting protons from a model box).Moreover, its rate does not change monotonically.Proton production due to advective-dispersive transport shows a maximum around river km 48 and a secondary maximum around river km 67.At the upstream boundary advectivedispersive transport accounts for 44% of the proton consumption, while downstream it accounts for about 25% of the proton production.
Figure 6b shows the longitudinal profile of the volumeintegrated proton production and consumption rates (expressed "per model-box").Table 8 lists selected values along the estuary.The cross section area increases from around 4000 m 2 upstream to around 76 000 m 2 downstream, while the mean estuarine depth remains approximately constant around 10 m.As a result, there is a much larger estuarine volume in the downstream model boxes than in the upstream model boxes.This volume increase per box compensates for the decrease in the rates per unit volume.As a consequence, the volume-integrated proton production or consumption rates in Fig. 6b remain similar along the estuary.The mid-region of the estuary (between kms 30 and 60) emerges as the most important region for volume integrated proton cycling.In this area the proton budget is dominated by the physical transport processes: CO 2 air-water exchange and advective-dispersive transport.The volume-integrated proton production/consumption of oxic mineralisation, primary production and CO 2 degassing is clearly larger downstream than upstream.In contrast, the volume-integrated proton production of nitrification is larger upstream than downstream.

Whole estuarine proton budget
Figure 7 shows a proton budget integrated over the whole model area and averaged over the four modelled years.It can be seen that CO 2 degassing and primary production are the processes that net consume protons in the estuary (disregarding the minor contribution of denitrification).Nitrification is the main proton producer, followed by oxic mineralisation and advective-dispersive transport.Note that the proton budget does not add up to zero, but a small negative value remains (−20 kmol [H + ] estuary −1 y −1 ).This indicates that the estuary as a whole, averaged over a year, is not com-pletely in steady state.This is consistent with the upwards trend in the mean estuarine pH, which is observed in the data and in the model simulations.
As shown in Fig. 6a, the relative importance of processes in the overall proton cycling changes along the estuary.To capture this in a proton balance, the estuary was divided into three zones; the upstream region between river km 0 and 30 which exhibits a pH below ca.7.6 (cf.Fig. 2), the midstream region between river km 30 and 60 where the pH exhibits a marked gradient, and the downstream region from river km 60 onwards where the pH remains approximately constant around 8. for each zone (Fig. 8).Because the decrease in the process rates from upstream to downstream is compensated by an increase in estuarine volume, the total proton production is of the same order of magnitude in the three zones.However, there are marked changes in the relative importance of processes.
Upstream, the proton production due to nitrification, and to a lesser extent oxic mineralisation, are counteracted by CO 2 degassing and advective-dispersive transport.In the midstream section, the contribution of aerobic respiration and advective-dispersive transport (which has now become a proton producer) are of similar magnitude as that of nitrification.In this midstream part, CO 2 degassing is also the major proton consuming process.In the downstream part of the estuary, nitrification is even less important, and the net proton production by oxic mineralisation as well as the net consumption of protons by primary production gain more importance.

Factors responsible for the trend in the mean estuarine pH from 2001 to 2004
Figure 9 shows the trend in the mean estuarine pH over the four year period.Both the trend as derived from the data as well as the trend simulated by the model are shown.The pH changed by about 0.085 from 8.010 to 8.095 in the model, which is backed by a similar pH increase of 0.095 in the data.There is a small offset between model and data.This is because the model slightly underestimates the pH in the upstream region with low estuarine volume and overestimates the pH in the downstream region with high estuarine volume.
Because of volume averaging, the latter dominates the mean estuarine pH.As a first step to investigate the underlying causes of this trend, we have plotted the volume integrated proton production/consumption for the individual years 2001 to 2004 alongside each other (the values displayed in Figs. 4, 5 and  6 were time averages over the whole four-year period).This shows that the overall proton turn-over decreased over the four-year period.The contributions of CO 2 degassing and nitrification steadily declined from 2001 to 2004 (with the decline being more pronounced for CO 2 degassing).The contribution from oxic mineralisation showed no clear trend, while the contribution of transport declined from 2001 to 2003 and then slightly increased again from 2003 to 2004.
As a second step, we investigated the factors responsible for the observed pH trend by simulating the various model scenarios in Table 6. Figure 10 shows the results of these different model scenarios.As discussed above, the loading changes of compounds can be caused by either changes in the freshwater discharge or changes in boundary concentrations.We examined these two effects separately.A large    .Blue line: pH data (NBS scale).NIOO monitoring data has been time-averaged using a center weighting scheme and results have been interpolated to every model box and then volume averaged.Red bars: indicator for the temporal and spatial variability of the modelled pH in the estuary.In every model box, the standard deviation of the temporal variability has been calculated, and these results are subsequently volume averaged over all model boxes.)portion of the pH change from 2001 to 2004 can be attributed to the decrease in the freshwater discharge alone (Fig. 10a).This decrease in the freshwater discharge generates a similar pattern of decrease in the overall proton cycling as found in the baseline simulation, albeit not on the same magnitude.Particularly, the decline in the proton production from nitrification is noticeable.The mean estuarine pH increases with 59 % of the increase in the baseline simulation.On the other hand, about 44% of the pH change can be attributed to the changes in boundary concentrations alone (Fig. 10h).However, the pH does not increase monotonously in this scenario.Note that pH changes due to freshwater discharge and boundary condition changes should not be necessarily additive.Furthermore, we can investigate which chemical species are dominating the long-term pH trend.A change in the loading of [ CO 2 ] and [TA] are most important.For these species, changes in the loading due to freshwater flow decrease account for 49% of the pH change (Fig. 10b), while changes in boundary conditions account for 28% of the pH change (Fig. 10i).Changes in the loading of [ NH +  4 ] are also influential and particularly modulate the proton production by nitrification (Fig. 10d and k).Due to a reduced freshwater flow, less ammonium is imported into the estuary, resulting in lower nitrification rates, especially in 2004 as compared to 2003.Lower nitrification rates, in turn, lead to a higher pH.Influences via [ NH + 4 ] are thus mainly indirect effects via changed nitrification rates.
The effect of decreased freshwater input on the salinity S does not lead to a pH increase but a decrease of 22% (Fig. 10c Fig. 10.Results for the model sensitivity scenarios given in Table 6.These simulations are used to investigate the factors governing the change in the mean estuarine pH from 2001 to 2004.See Fig. 9 left side for legend.salinity, which then stimulates outgassing, resulting in a decrease of the pH over time.As shown in Fig. 10e, f, g, j, l,  m, and n

pH modelling in aquatic systems
Estuarine systems, like the Scheldt estuary, are characterized by strong geochemical cycling, and hence, intense consumption and production of protons.Accordingly, a true challenge is to develop accurate pH models for such dynamic systems.Biogeosciences, 6, 1539Biogeosciences, 6, -1561Biogeosciences, 6, , 2009 www.biogeosciences.net/6/1539/2009/ In the past, a number of reactive transport models have been developed that accurately reproduce the longitudinal pH profile of the Scheldt estuary (Regnier et al., 1997;Vanderborght et al., 2002;Hofmann et al., 2008b).However, these models implemented the implicit pH modelling approach, and so, they were not able to quantify the contribution of specific biogeochemical processes to the overall proton cycling.One (crude) way to investigate the importance of a single process in the overall proton cycling, is to switch this process off, and examine the effect on pH.However, because processes interact, switching off a given process will influence the rates of the other processes, and hence, this could complicate the interpretation.Preferably, one wants a method that quantifies the contribution of individual processes concurrently and without changing the rates of biogeochemical processes.
The explicit approach to pH modelling allows just that.This explicit approach was originally pioneered by Jourabchi et al. (2005) and Soetaert et al. (2007), though these treatments were partially incomplete.The approach proposed in Jourabchi et al. (2005) is only applicable to steady state conditions, and it also treats the effect of advective-diffusive transport on proton cycling incompletely (the proton transport term is omitted).Similarly, Soetaert et al. (2007) introduced a method to quantify the influence of processes on pH, one process at a time, but this approach did also not account for advective-dispersive transport (which is important for proton cycling as we show here).
Recently, Hofmann et al. (2008a) have reviewed the construction of pH models, and advanced an explicit pH modelling approach that allows to quantify the contribution of all processes (that is physical transport and biogeochemical reactions) on the proton cycling within an aquatic system.This approach describes the pH evolution explicitly using an expression for the rate of change of the proton concentration over time.However, the approach advanced in Hofmann et al. (2008a) assumed that dissociation constants remained constant over time.In previous modelling studies of the Scheldt estuary (e.g.Vanderborght et al., 2002;Hofmann et al., 2008b), it was shown that the dissociation constants must be calculated as functions of the time-variable salinity and temperature in order to obtain reasonable fit to the pH data.In response to this, we have here extended the explicit method so that the dissociation constants can change over time.Note that both the method given in Hofmann et al. (2008a) and the generic treatment presented here are applicable to both water column and sediment pore water ecosystems.

pH profile along the estuary
In the Scheldt estuary, pH increases from upstream to downstream.In Fig. 4 we show that this is an effect of a decreasing [ CO 2 ] / [TA] ratio from upstream to downstream, which itself is an effect of conservative mixing of [ CO 2 ] and [TA] with riverine and marine end-members.However, the pH increase induced by this conservative mixing exhibits no distinct minimum at low salinities and no change in curvature.The distinct shape of the pH increase observed in the Scheldt estuary and other estuaries, exhibits a clear minimum at low salinities and has a distinct quasi sigmoidal shape.Accordingly, the pH profile along the estuary must be determined by other factors.
Over the past few decades a discussion has taken place about what these controlling factors could be.Mook and Koene (1975) suggested that the characteristic pH profile in estuaries with high inorganic carbon loadings from upstream (like the Scheldt estuary) only results from acid-base equilibration following the mixing of fresh water and seawater, due to the salinity dependency of the dissociation constants of the carbonate system.The conceptual model of Mook and Koene (1975) is very simple, it only accounts for conservative mixing of [ CO 2 ] and [TA] and salinity dependent acid-base equilibration.Accordingly, it ignores CO 2 exchange with the atmosphere as well as biological processes (e.g.aerobic respiration, primary production, nitrification) that may affect [TA] and [ CO 2 ].The blue line in Fig. 5 represents the model of Mook and Koene (1975).It confirms that in a closed system with no biology a pH minimum at low salinities is created by the salinity dependency of acid-base dissociation constants (as the pH profile in the case with salinity independent dissociation constants does not show a minimum -Fig.4).However, this does not exclude that other processes could be important in shaping the pH profile along the estuary.
Although the model of Mook and Koene (1975) provided reasonable agreement with measurements in the Chesapeake Bay (Wong, 1979), and is still used to predict estuarine pH profiles (Spiteri et al., 2008), it is thought to be a rather crude approximation of reality.Whitfield and Turner (1986) showed that CO 2 exchange with the atmosphere results in significantly different pH profiles.Differences up to 0.7 pH units were found at low salinities for systems that are fully equilibrated with the atmosphere.Enabling gas exchange while keeping biological processes inactivated in our model results in a concave pH profile with high pH values and no minimum (red line in Fig. 5).As also argued by Whitfield and Turner (1986), this pH profile strongly deviates from the one that considers no gas exchange, but it severely overestimates observed pH values.This shows that neglecting gas exchange introduces an error into the pH calculations in the Scheldt estuary.Moreover, the effect of gas exchange on pH profile is stronger than the salinity dependency effect of the acid-base dissociation constants.
As mentioned above, estuaries are zones of intense biological activity, so neglecting biological processes in pH calculations seems unjustified.The strong effect of biology (or equally biogeochemistry) is shown by the magenta line in Fig. 5).This simulation produces a convex pH profile with a clear pH minimum at low salinities, which strongly underestimates observed pH values.Finally, a full model simulation with all biological processes and gas exchange toggled on (black line in Fig. 5) again shows a distinct pH minimum at low salinities, but the resulting pH profile follows much more closely the observed data.This shows that the conceptual models which ignore gas exchange or biology, oversimplify the pH dynamics in an estuarine system.Effectively, the conceptual model of Mook and Koene (1975) predicts a pH minimum at low salinities similar to the one observed in the data, but for incorrect reasons.
In summary, one can conclude that the overall pH increase in estuaries similar to the Scheldt estuary (high CO 2 loadings) results from a decreasing [ CO 2 ]/[TA] ratio between the riverine and marine end-members.The characteristic sigmoidal shape of the pH increase results from an interplay between CO 2 outgassing and biological activity.Overall, biological activity has a clearly acidifying effect.Nitrification seems to be the main cause for the distinct dip in pH values at low salinities.The salinity dependency of acid-base dissociation systems only has a minor effect on the characteristic pH profile in estuaries.
Note that a quantitative mathematical link between spatial pH gradients along the estuary and process rates at different points of the estuary will be the subject of another publication.

Proton production and consumption along the estuary
The pH profile in the Scheldt estuary (Fig. 2) is the result of a balance between proton production and proton consumption processes.Proton production P P and the proton consumption P C (expressed per volume unit of solution -Eq.25) are very large compared to the net proton production rate.Accordingly, the magnitude of either P P or P C can be used as a measure for the intensity of proton cycling.This proton cycling intensity decreases by a factor of 20 along the estuary, going from around 1.3 mmol of protons m −3 yr −1 at the upstream boundary to 0.06 mmol of protons m −3 yr −1 downstream.The turnover time of protons can be defined as the local proton concentration divided by the proton cycling intensity.Using pH values of 7.55 upstream and 8.05 downstream, we find that the turnover time of the proton pool is 8 days upstream and 54 days downstream.This decrease in the intensity of proton cycling matches the decrease in the cycling intensity of total ammonium, dissolved inorganic carbon, oxygen, and nitrate as calculated in Hofmann et al. (2008b).This shows that high proton turnover is associated with high rates in biogeochemical processes, thus stressing the importance of such processes in the pH dynamics of estuaries.Moreover, the rates of proton production and consumption presented here (Fig. 6) are approximately the mirror images of the rates of oxygen production and consumption presented in Hofmann et al. (2008b).This suggests a tight connection between redox balance and proton cycling, which should be explored in future investigations.
Although the link between high biogeochemical rates and high proton turnover holds for all processes combined, it does not always apply to individual processes.The prime example is denitrification, which shows a negligible contribution to the overall proton cycling.Still, based on the baseline simulation, denitrification can be shown to be an important process in the overall cycling of NH + 4 and CO 2 (Hofmann et al., 2008b).Judging from the stoichiometry (Table 1), denitrification substantially consumes protons.Based on this, one would expect denitrification to be important in the proton cycling.The reason why this logic fails is that one has ignored acid-base equilibration: denitrification also produces CO 2 , a proton source.The protons subsequently released due to dissociation of produced CO 2 almost balance the direct proton consumption of denitrification.This illustrates the value of the explicit pH modelling approach: it is able to exactly quantify the contribution of a given process to proton cycling (properly accounting for the effects of acidbase equilibration).
Nitrification and oxic mineralisation are the important proton producers, and CO 2 degassing, and to a lesser extent primary production, are the important proton consuming processes.This is fully consistent with previous modelling studies on the Scheldt estuary (i.e., Regnier et al., 1997).Furthermore, these results are coherent with Abril and Frankignoulle (2001), who focussed on alkalinity variations and identified the nitrogen cycle and the process nitrification as important factor governing the acid-base regime of the Scheldt estuary.In addition to confirming this, the explicit pH modelling approach presented here allows to precisely calculate the contribution of nitrification (and other processes) in the overall proton cycling along the estuarine gradient.
The contributions of changes in the dissociation constants in the overall proton cycling are several orders of magnitude smaller than the influences of kinetic processes (Tables 7 and  8).Therefore, when describing the factors that govern the gross proton cycling intensity (P P , P C), these factors can be neglected.However, to model pH values in a system accurately, i.e. more accurate than 0.1 pH units, they should be included in the calculations.This is especially important for modelling the proton concentration explicitly over a longer period of time, since deviations in d[H + ] dt are likely to accumulate.

Proton cycling and CO 2 degassing
Estuaries are important sites of CO 2 degassing.Within the context of climate change and greenhouse gas budgets, it is important to understand the control on the rate of CO 2 degassing.The baseline simulation estimates a total CO 2 efflux of 3.3 Gmol yr −1 over the Scheldt estuary (averaged over the 2001-2004 year period, see Hofmann et al., 2008b, for details).A complex interplay between different factors controls this rate of CO 2 degassing.We performed simulations to examine the influence of the biogeochemical processes on CO 2 degassing.To this end, we toggled off all biological processes in the baseline simulation.In this scenario, the CO 2 efflux rates are reduced by approximately half over the whole estuary.While this shows that biological processes induce a large amount of CO 2 degassing, it also indicates that there is substantial CO 2 degassing driven by non-biological mechanisms.Most importantly, the incoming water is significantly oversaturated in CO 2 (around 300%, see Hofmann et al., 2008b), which induces degassing during transport along the estuary.
Among the biological processes, Frankignoulle et al. (1996), Abril et al. (2000), and Vanderborght et al. (2002) identified nitrification and oxic mineralisation as the important processes increasing CO 2 degassing in the Scheldt, and argued that these processes are interconnected via acid-base chemistry.Our analysis confirms this, but provides additional quantitative information: e.g. the proton production by nitrification is made explicit along the estuary, showing that its importance gradually decreases downstream (see Figs. 7  and 8).
Mechanistically, the control of acid base chemistry on CO 2 outgassing is complex: the kinetic rate of CO 2 outgassing is pH dependent, so processes that produce protons and lower the pH stimulate CO 2 outgassing.Conversely, processes that consume protons prevent CO 2 outgassing.For example, nitrification produces protons and stimulates CO 2 outgassing.However, some processes have mixed effects.As mentioned above, denitrification explicitly consumes protons (an effect that diminishes CO 2 outgassing) but also produces CO 2 .The latter effect increases CO 2 outgassing directly by increasing [ CO 2 ] and CO 2 oversaturation, and indirectly via the protons produced by CO 2 dissociation.This complexity warrants further research.The tools provided in this publication allow for a further investigation of the control of acid base chemistry on CO 2 outgassing.

Advective dispersive transport and proton cycling
The contribution of advective-dispersive transport to the rate of change of proton concentration changes sign around river km 32.Advective-dispersive transport is a proton consumer upstream and a proton producer downstream.This is due to the fact that the amount of protons produced by nitrification and oxic mineralisation in the upstream region of the estuary cannot be sufficiently compensated by proton consumption through CO 2 degassing.This causes a net downstream transport of protons, or equally, a net export of protons by advective-dispersive transport in the upstream model boxes.Downstream of river km 32, nitrification and oxic mineralisation reach lower levels.As a result, their proton production is now overcompensated by proton consumption by CO 2 degassing, which results in a net import of protons in the model boxes there.This import of protons declines towards the downstream boundary of the model.A thorough investigation of the role of advective-dispersive transport for the pH in the estuary, however, is beyond the scope of this publication.Changes in freshwater flow and boundary conditions influence the estuarine pH "directly" via influences on [ CO 2 ] and [TA].By this, we mean that a reduced freshwater flow imports less CO 2 , a proton source, into the estuary, resulting in a higher pH.However, there is also an "indirect" influence via [ NH + 4 ] and the nitrification rates in the estuary.This "indirect" pathway is about half as important as the "direct" influences via [ CO 2 ] and [TA].

Conclusions
In this publication, a novel 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 the conventional pH modelling approach.The presented method is generic and can be applied to any aquatic system, including the pore water of sediments.As a case study, the presented approach was applied to assess proton cycling in the Scheldt estuary.We showed that the pH increase along the Scheldt estuary is a result of a decreasing CO 2 /[TA] ratio due to mixing of riverine and marine end-member water masses.The distinct sigmoidal shape of the pH increase is the result of an interplay between CO 2 degassing and biology.The salinity dependency of acid-base dissociation constants only plays a minor role.Nitrification has been identified as the dominant biological process affecting pH and proton cycling in the Scheldt estuary.
The proton cycling intensity in the estuary significantly drops from upstream to downstream, which mirrors the decrease in cycling intensity for total ammonium, dissolved inorganic nitrogen, oxygen, and nitrate.This confirms the importance of biological processes in the pH chemistry in estuarine systems.The importance of different processes for total proton cycling changes along the estuary: upstream nitrification can be identified as most important, while midstream and downstream CO 2 degassing is dominant.In the whole estuary, nitrification and oxic mineralisation are the most important proton producers, while CO 2 degassing, and www.biogeosciences.net/6/1539/2009/Biogeosciences, 6, 1539-1561, 2009 to a lesser extent primary production, are the most important proton consuming processes.Advective-dispersive transport plays a special role as it changes its sign: it is a proton producer upstream and a proton consumer downstream.A clear inverse correlation between oxygen and proton turnover was found, consistent with theoretical considerations of redox chemistry in oxic waters.
The main driver of changes in the mean estuarine pH over the period from 2001 to 2004 is a changing freshwater flow.The pH is influenced "directly" via a changed loading of [ CO 2 ] and [TA] and also -to a significant amount -"indirectly" via [ NH + 4 ] and changing nitrification rates in the estuary.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.Note also that Eq. ( 12) 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. ( 5) and ( 11 [TA]= This means 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 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)

Fig. 1 .
Fig. 1.The Scheldt estuary.Gray dots represent positions in the river where longitudinal profiles of influences of processes on pH, presented in the Results section, show interesting features.

Fig. 2 .
Fig. 2. The model fit for pH for the modelled years 2001 through 2004.The dots represent NIOO monitoring data on the NBS scale (see Hofmann et al., 2008b), the lines represent modelled pH (modelled free scale pH values have been converted to the NBS scale using the Davies equation -cf.discussion of the use of the Davies equation in Hofmann et al., 2008b).

Fig. 3 .
Fig.3.Verification of the explicit pH modelling method using model runs for the year 2003.The black dots represent NIOO monitoring data (seeHofmann et al., 2008b), the black lines represent modelled pH calculated with the implicit approach and the blue lines represent modelled pH calculated with the explicit approach (modelled free scale pH values have been converted to the NBS scale using the Davies equation (cf.discussion of the use of the Davies equation inHofmann et al., 2008b)).(a): omitting all K * related terms from Eq. (14); (b): considering the terms describing the variations in the dissociation constants due to changes in S and T but without the pH scale conversion related terms; (c): considering all terms as described in Sect.(2.4.3).

Fig. 4 .
Fig. 4. Basic model simulation with no salinity dependence of the acid-base constants (values were calculated at T =15 • C and S=15), no gas exchange, and no biogeochemistry.[ CO 2 ], [TA] and pH (NBS) are plotted along the estuary as a function of salinity.The pH from upstream to downstream increases as a result of a decrease in the ratio between [ CO 2 ] and [TA].End-members for [ CO 2 ], [TA], and S are the upstream and downstream values for 2001 in the baseline model run.Due to the simplicity of the remaining model formulation the calculations are performed using an R script (R Development Core Team, 2005) and the acid-base modelling package AquaEnv (Hofmann et al., submitted) instead of the full FORTRAN implementation of the model.
estuarine pH from 2001 to2004Hofmann et al. (2008b) ) reported an upward trend in the mean estuarine pH over the years 2001 to 2004, but did not investigate the underlying causes of this trend.Potential factors are a different temperature forcing, changes in freshwater input (modifying the salinity gradient and the concentrations of all chemical species in the estuary), or temporal variations in the chemical composition of the water at the boundaries.Table

Fig. 6 .
Fig. 6.The influences of kinetically modelled processes on the pH -volumetrically a) and volume integrated b), averaged over the four modelled years.Note that the process abbreviations given in the legend represent the contribution of the respective process to d[H + ] dt , e.g., Tr signifies d[H + ] dt Tr .Process abbreviations: Tr=advective-dispersive transport, E CO 2 =CO 2 degassing, R Ox =oxic mineralisation, R Nit =nitrification, R PP =primary production, K * (T )=changes in dissociation constants via changes in temperature, K * (S)=changes in dissociation constants via changes in salinity, K * ([ HSO − 4 ])=changes in dissociation constants via changes in [ HSO − 4 ], K * ([ HF])=changes in dissociation constants via changes in [ HF].

Fig. 7 .
Fig. 7. Whole estuarine proton budget, averaged over the four modelled years.The error bars represent the standard deviations resulting from averaging over the four years.The process abbreviations in the legend denote the contribution of the respective biogeochemical process e.g., Tr refers to d[H + ] dt Tr .

Fig. 8 .
Fig. 8. Proton budget for three different zones in the estuary: (a) upstream region between river km 0 and 30, (b) midstream region between km 30 and 60, and (c) the downstream region, averaged over the four modelled years.The error bars represent the standard deviations resulting from averaging over the four years.The process abbreviations in the legend denote the contribution of the respective biogeochemical process e.g., Tr refers to d[H + ] dt Tr .

Fig. 9 .
Fig. 9. Left panel: Proton budget over the individual years 2001 to 2004.The process abbreviations in the legend denote the contribution of the respective biogeochemical process e.g., Tr refers to d[H + ] dt Tr .Right panel: trend in the mean estuarine pH over the period 2001 to 2004.Black line: model simulated pH (free scale values have been converted to the NBS scale).Blue line: pH data (NBS scale).NIOO monitoring data has been time-averaged using a center weighting scheme and results have been interpolated to every model box and then volume averaged.Red bars: indicator for the temporal and spatial variability of the modelled pH in the estuary.In every model box, the standard deviation of the temporal variability has been calculated, and these results are subsequently volume averaged over all model boxes.) ).A decreased freshwater input increases the www.biogeosciences.net/6 , changes in the loading the organic matter fractions, [O 2 ] and the rest of the state variables ([NO − 3 ], [ HSO − 4 ], [ B(OH) 3 ], and [ HF]) have a minor impact on the pH trend.
responsible for the change in the mean estuarine pH from 2001 to 2004 Interannual changes in freshwater flow and boundary concentrations (i.e.changes in the loading of chemical species) are the driving forces for changes in the whole estuarine mean pH over the years 2001 to 2004.The observed increase in mean estuarine pH from 2001 to 2004 can be attributed to a decrease in the freshwater flow Q, which is consistent with Hofmann et al. (2008b).Changes in the boundary concentration values further re-enforce this trend.
of [TA] with respect to equilibrium invariantsFor the work presented here, the partial derivatives of [TA] with respect to the equilibrium invariants (i.e. the terms∂[TA] ∂[X j ]) have been calculated analytically: ) are kept constant.That means, e.g. in the term ∂[TA]   ∂[ HSO − 4 ] the dissociation constants are considered constants, although they are also functions of[ HSO − 4 ].Likewise, in ∂[TA] ∂K * i , [ HSO − 4 ] is considered constant, while for ∂K * i ∂[ HSO −[ B(OH) 3 ] independently from the salinity S (although borate species contribute to S).Therefore, for∂[TA] ∂[ B(OH) 3 ] ,S is considered a constant, although, strictly speaking, changes in [ B(OH) 3 ] would also change S.This is done to mathematically separate influences of changes in S via the dissociation constants on [TA] and changes in the equilibrium invariant [ B(OH) 3 ] on [TA] directly.] can be evaluated analytically.Consider a system where total alkalinity equals carbonate alkalinity[TA]=[HCO − 3 ]+2[CO 2− 3 ] (A15)
data together with the pH boundary data in an inverse way to calibrate the [ CO 2 ] boundary conditions for the year 2003.Starting from an initial guess for the [ CO 2 ] boundary values, model predicted [TA] values were compared with the

Table 2 .
Mass balance equations for the state variables in the model.R Ox FastOM and R Ox SlowOM

Table 3 .
Left : acid-base equilibria in the model.Right: definition of dissociation constants (stoichiometric equilibrium constants).Values are calculated from temperature, salinity and pressure based on literature expressions (K *
Hofmann et al. (2008b)Scheldt estuary salinity gradient.The blue line represents the pH calculated with a closed system model without biology (comparable toMook and Koene, 1975); the red line represents the pH calculated with an open system model without biology (comparable to Whitfield and Turner (1986) but with realistic kinetic CO 2 air-water exchange instead of a fully equilibrated system); the magenta line represents a closed system model with biology; the black line represents the pH calculated with the full biogeochemical model as presented inHofmann et al. (2008b).All models are based on 2001 parameter values.The black circles are the observed pH data for 2001.

Table 5 .
Model forcings (upper part of the table) and a selection of mean estuarine model values resulting from the baseline simulation (lower part of the table).The subscript "up" denotes upstream boundary condition , the subscript "down" denotes downstream boundary condition.The boundary conditions for [TA] are calculated from boundary values for [ CO 2 ] in combination with pH up and pH down (which are otherwise not used as boundary conditions).Concentrations are given in mmol m −3 , the flow at the upstream boundary Q is given in m 3 s −1 .The pH output from the model simulations is converted to the NBS scale to be comparable with the data.

Table 6 .
Model scenarios to investigate the trend in the mean estuarine pH over the years 2001 to 2004.Changes in the total "loading" for particular chemicals either due to changes in the freshwater discharge (left column) or due to changes in the boundary concentrations (right column).The entries indicate the variables for which the loading has changed, while all other loadings have been fixed at 2001 values.Note that in all these scenarios, [TA] boundary conditions are calculated consistently from pH boundary forcing values.
Hofmann et al.:pH modelling with time variable acid-base constants at low salinities, the data pH profile shows less curvature at high salinities).

Table 7 .
Contributions of various biogeochemical processes to the proton cycling per unit of solution volume; values in mmol H + m −3 y −1 ; percentages are of total production (positive quantities) or consumption (negative quantities), respectively.

Table 8 .
Volume integrated proton budget; proton production/consumption rates in kmol H + y −1 per model box.