Biogeosciences A step-by-step procedure for pH model construction in aquatic systems

We present, by means of a simple example, a comprehensive step-by-step procedure to consistently derive a pH model of aquatic systems. As pH modelling is inher­ ently complex, we make every step of the model generation process explicit, thus ensuring conceptual, mathematical, and chemical correctness. Summed quantities, such as total inor­ ganic carbon and total alkalinity, and the influences of mod­ eled processes on them are consistently derived. The dif­ ferent time scales of processes involved in the pH problem (biological and physical reactions: days: aquatic chemical reactions: fractions of seconds) give rise to a stiff equation system. Subsequent reformulations of the system reduce its stiffness, accepting higher non-linear algebraic complexity. The model is reformulated until numerically and computa­ tionally simple dynamical solutions, like a variation of the operator splitting approach (OSA) and the direct substitution approach (DSA), are obtained. As several solution methods are pointed out, connections between previous pH modelling approaches are established. The final reformulation of the system according to the DSA allows for quantification of the influences of kinetic processes on the rate of change of proton concentration in models containing multiple biogeochemical processes. These influences are calculated including the ef­ fect of re-equilibration of the system due to a set of acid-base reactions in local equilibrium. This possibility of quantify­ ing influences of modeled processes on the pH makes the end-product of the described model generation procedure a powerful tool for understanding the internal pH dynamics of aquatic systems. Correspondence to: A. F. Hofmann (a.hofmann@nioo. knaw. nl)


Introduction
Human activities have increased atmospheric CO2 levels by 36% since pre-industrial times, and further increases are ex pected over the next decades (Prentice et al., 2001: Alley et al., 2007).Rising atmospheric CO2 levels lead to an input of CO2 into the oceans and to subsequent acidification of sur face waters (e.g.Orr et al., 2005).
Against this background, it is of high importance to ana lyze the impact of different biogeochemical processes onto alkalinity and the pH of natural waters (Sarmiento andGru ber, 2006: Stumm andMorgan, 1996).In recent years, var ious pH modeling approaches have been developed.These range from simple empirical correlations (Bjerknes and Tjomsland, 2001), over neural network approaches (Moatar et al., 1999), to mechanistic biogeochemical models that in clude reactive transport descriptions of varying complexity (e.g.Luff et al., 2001: Jourabchi et al., 2005).Mechanistic models have the advantage that they not only reproduce pH but also allow the prediction of future changes, and enable quantitative analysis of the processes that govern pH.As a result, they are a powerful tool to understand the pH dynam ics of aquatic systems.
However, there are still two pending problems with mech anistic pH models.The first issue relates to the apparent di versity of approaches.Most modeling approaches have been presented without cross linking to other methods.As a result, it is difficult to assess whether approaches are mutually con sistent, i.e., whether they would predict the same pH dynam ics for exactly the same underlying biogeochemical model.Moreover, it is not clear what the respective advantages of the different solution techniques are, and whether they yield the same amount of information with respect to pH dynam ics.Only some approaches are able to quantify the individual contribution of modelled processes on the pH.
Published by Copernicus Publications on behalf of the European Geosciences Union.The second issue relates to the complexity of the present approaches.The construction of pH models is inherently complex, involving many sequential steps and assumptions.Furthermore, the different time scales of processes involved in the pH problem (biological and physical reactions: days; aquatic chemical reactions: fractions of seconds) give rise to a stiff equation system.It is important to deal with this complexity by making every assumption explicit and ju sti fying every step.Even for a relatively simple biogeochem ical system, the model generation procedure becomes quite lengthy and intricate.A disadvantage of recent pH model ing approaches is that they have been typically applied to complex reaction sets, generating lengthy expressions.The illustration of a complex solution procedure with a complex model is not always optimal.Accordingly, there is a clear need to illustrate the various approaches to model pH with a simple biogeochemical application.
The objective of the present study is to provide a generic step-by-step procedure to construct and solve a pH model for an aquatic system.We will illustrate this step-by-step ap proach using an example, i.e., by constructing an example pH model for a simple estuarine system.This example is simple enough to facilitate understanding, yet complex enough to illustrate all features of the pH modeling approach.Accord ingly, the focus lies on concepts and principles rather than on mimicking the biogeochemical complexity of real aquatic systems.Models of more realistic and complex systems can be built by suitably changing the transport formulation or ex tending the reaction set.The feature of our analysis is that we carry out a number of sequential reformulations of the pH problem until elegant and efficient numerical solutions are possible.Along the way, we outline the implicit and explicit assumptions that are needed in every step of the procedure.This enables us to identify the weaknesses and strengths of past modeling procedures and solution methods.Our work therefore does not introduce a novel approach to pH mod elling, but gives a systematic framework which encompasses existing approaches.
2 pH model construction: a step-by-step procedure 2.1 Step 1: Formulation of the model questions Our example system is the upper Schelde estuary in north ern Belgium (Fig. la).The model domain includes 40 km of river ranging from the inflow of the Rupel tributary to the Belgian-Dutch border.A set of characteristic parameters is given in the Results section.Our principal goal is to examine the pH changes associated with some (drastic) perturbations in the biogeochemistry of this estuary.Two types of changes to the system are examined: 1) The estuary receives municipal water from the city of Brussels, which is one of the last major European cities to implement a coordinated waste water treatment pol icy.In 2007, a new sewage treatment plant for 1.1 mil lion inhabitants has started operating, and it is estimated that this will reduce the organic matter input to the estu ary by 50%.How will the pH of the estuary react to this abrupt change?W hich biogeochemical processes gov ern the pH steady state before and after the reduction?
2) Alongside the estuary lies the port of Antwerp, which concentrates one of the largest chemical industries in the world.The port harbours a large fertilizer industry with associated ship traffic of resources and products.Potential hazard scenarios include ship accidents with tankers carrying ammonia or ammonium-nitrate.What are the effects of such pulse-inputs on the estuarine pH and the influences of processes on it?Ammonia input and ammonium-nitrate input are examined as two sepa rate perturbation scenarios.

2.2
Step 2: Formulation of the conceptual model In general, the concentration of a chemical species [X] in an aquatic system is influenced by a set of physical (trans port) processes P; , and a set of biogeochemical reactions R!.
The resulting mass conservation equation (MCE) (Morel and Hering, 1993) where is the stoichiometric coefficient of species X in the i-th reaction.Throughout this paper, all species concentra tions [X] are expressed as per kg of s o lu tio n (gravimetric units per mass of solution).A crucial step in the model de velopment is the decision which physical and biogeochemi cal processes to include in the model.This decision should be based on prior knowledge about the physics and biogeo chemistry of the system.For our problem, two physical pro cesses are of major relevance: advective-dispersive transport T x along the length axis of the river, and the exchange of volatile compounds with the atmosphere Ex-Although it Table 1.Estimated rates (/xmol-lV (kg • d)" 1) of biogeochemical processes in the example system ((a): Soetaert and Herman, 1995a: (b): Andersson et ah, 2006: (c) : Soetaert and Herman, 1995b: (d): Middelburg et al., 1996).pelagic primary production Rpri " « 0 .1 (a) pelagic nitrification Rnit ?« 7.5 (b) pelagic denitrification Rden ?« 6 .1 (c) pelagic oxic respiration Rox ' « 2.9 (c) benthic denitrification Rbden ?« 0.7 (c) benthic oxic respiration Rbox ?« 0.3 (d) would be possible here, no benthic exchange is taken into ac count to keep the mathematical expressions tractable.For the same reason, the estuary is modelled as a single box (Fig. lb).
Note that the implementation of a spatially explicit descrip tion would be entirely analogous in terms of pH modeling: the T x terms would simply give rise to partial rather than ordinary differential equations.
To assist in the selection of the biogeochemical reactions, Table 1 provides an overview of the relative importance of the various processes in the Schelde estuary.From the six biogeochemical reactions listed, we only retain pelagic oxic respiration R ox and pelagic nitrification Rnit.These are de scribed according to reaction stoichiometries: With y being the C/N ratio of organic matter (see Ta ble 14).
Pelagic primary production, benthic denitrification and benthic respiration can be justifiably neglected compared to pelagic nitrification and pelagic oxic mineralization (Ta ble 1).Pelagic denitrification was important in the 1970's, but, due to improved water quality, is now of minor signifi cance (Soetaert et al., 2006).For this reason and for didac tical purposes, we did not include it in the model (to avoid lengthy expressions in the mathematical derivations).How ever, we will include it a posteriori to check on the impor tance of denitrification in pH regulation of the model domain.
Since our aim is to model the pH, a number of acid-base reactions have to be accounted for.To select these reactions, we first have to constrain the set of chemical species that are modelled.For simplicity, we consider the estuary as an aque ous solution of the three most abundant seawater ions C1" , Na+ , and SO^" (DOE, 1994).For a realistic model applica tion other quantities like borate might be important, but we neglect these to keep the model as simple as possible.Fur thermore, we also incorporate organic matter, nitrate, oxygen and the ammonium and carbonate systems, as these species feature in the retained reactions (Rox and Rnit, Eq. 2).Ta ble 2 lists the set of acid-base dissociation reactions R(*ls that Table 2. Acid-base reactions in the example system, thermody namical pK ha's are infinite dilution values at 2 5° C as given in Stumm and Morgan (1996).According to the exclusion criterion given in Appendix A, reactions with an e below 0.5 % are ne glected.e has been calculated for a desired pH range of 6 to 9, with p K ^p K H A -and with [TA.] =5000 ¡x mol kg-1 (estimated from upstream and downstream boundary conditions given in Table 14), and with total concentrations for the given system as listed (total nitrate and ammonium are measured values for the example model system, total carbon dioxide has been estimated and all other to tal quantities have been calculated from salinity S=5 according to DOE, 1994).
After the selection of the pH scale, we can proceed to a formal delineation of the pH range of the model.This set ting of the pH range determines which dissociation reactions should be incorporated.Note that most pH modeling ap proaches do not explicitly mention this step.In these, the set of acid-base reactions is simply imposed without further consideration.However, models are simplified representa tions of reality, and they should be kept as simple as possible.This is particularly true for pH models, which are computa tionally demanding.Accordingly, one should avoid incorpo rating dissociation reactions that have no chance of affecting the pH dynamics.
Therefore, we propose a formal procedure for the selec tion of acid-base reactions which is based on prior knowl edge about the buffering capacity and the possible pH range of the specific system.In our case, we know that the part of the Schelde estuary which we model is strongly buffered, as are most estuarine and marine systems, with a total al kalinity [TA] of «50 0 0 /zmol kg-1 (estimated from upstrean and downstream boundary conditions given in Table 14).We furthermore know that the pH only fluctuates over a range from 7.5 to 8 .Nonetheless, we anticipate stronger excur sions because of the quite drastic perturbation scenarios out lined above.Allowing a suitable margin, we require that the model should represent the pH dynamics properly within a pH range of 6 to 9.This constraint enables us to reduce the reaction set in Table 2 considerably.
Whether or not a certain acid-base system has to be taken into account depends on a combination of 1) its pK value (s) which tell us whether the spéciation of the acid-base species will change within the pH range under consideration, 2) the total concentration of the acid [J^ A] which tells us how large theoretical changes in [H+ ] due to the spéci ation of the respective acid-base system would be in a completely unbuffered system, and 3) the mean [TA] of the system which tells us if these theo retical changes in [H+ ] will be appreciable or negligible in a buffered aquatic system.
Appendix A details a formal selection procedure which in tegrates these three criteria into a single quantity e for each acid-base system, e represents the amount of protons ignored (theoretical unbuffered proton concentration offset) by ne glecting the reaction in question, in percent of the average [TA] of the modeled system.
Finally, we exclude all reactions whose e value is smaller than 0.5%.This means the total amount of protons which could be taken up or released in the model, if the reaction in question would be included and the pH reaches the border of the pH range, is less than half a percent of typical alkalinity levels of the system.
Note that polyprotic acids are treated as a set of monoprotic acids considering each dissociation step independently.
Applying this rule (e values are given in Table 2), we do not need to incorporate the dissociation reactions of HC1, NaOH, H2 SO4 , H SO J, HNO3 and H2 O. Table 3 shows the reduced set of acid-base reactions considered in the model.Technically, it would not be "w rong" to include the other re actions.However, there is no reason to do so, provided that  Note that the auto-dissociation reaction of water is not included in Table 3. Effectively, this reaction has been treated in a rather arbitrary fashion in past models.The auto-dissociation of water is included in some models (e.g.Jourabchi et al., 2005), while excluded from others (e.g.Luff et al., 2001).Usually, the reasons for inclusion or exclusion are not mentioned.Here however, our formal selection pro cedure predicts that it is unimportant (we will check this a posteriori).Note that organic matter (CH2 0 ) y (NH3) has been ab breviated by OM and that the concentrations of Cl-, Na+ , H SO J, S 0 4-, NaOH and OH-are not simulated, since they are not affected by the modeled processes1.

Step
Although H20 does feature in the biogeochemical reac tions retained in the model (R ox and R nit; Eqs.(2 a) and (2 b)) and in the set of acid-base reactions (Table 3), its concentra tion is nevertheless considered constant (Morel and Hering, 1993).The resulting mass conservation equations for all 9 chemical species are given in Table 4, where again T x de notes advective-dispersive transport of chemical species X and Ex denotes the exchange of chemical species X with the atmosphere.
At this point, a first attempt to solve the system can be made.

Solution method [la]:
Together with suitable initial condi tions, the equation set in dt problem of ordinary differential equations (ODEs) (Fabian et al., 2001).Using suitable kinetic expressions for all mod eled process rates (i.e. for the forward and backward rates for acid-base reactions), this system is fully determined.In principle, it can be directly solved by common numerical integration techniques, such as Euler or Runge-Kutta in tegration (Press et al., 1992) or more complex integration schemes.This solution procedure is referred to as the Full Kinetic Approach (FKA) (Steefel andMacQuarrie, 1996: Meysman, 2001).
For the reaction f°r example, suitable kinetic ex pressions for the forward and backward reaction would be R dis forward p d is NH+ , and k f and forward V ' ''4 / backward k b being the forward and backward rate constants.Zeebe and Wolf-Gladrow (2001) give formulations for forward and backward rate constants of some acid-base systems relevant in seawater.However, problems arise when the values of the rate constants are not available.
Solution method [lb]: One way to avoid this problem is to adopt the principle of microscopic reversibility or detailed balancing (Morel and Hering, 1993), which requires the quotient of the kinetic forward and backward rate constants k f and k b to be equal to the equilibrium constant of the   (2001), and R ^ from Morel and Hering (1993).Values for the HCO3 remaining processes are estimated from Tables 1 and 14.For the exchange with the atmosphere, piston velocities Kp as given by Raymond and Cole (2001) were used.
(1) Rox Based on this principle, e.g., R ^f.+ can be written as Jmh4 which only features the forward rate constant kf.Assuming a sufficiently high value for kf this is an approximation of the local equilibrium assumption (Steefel and MacQuarrie, 1996) which will be discussed later.
Although it overcomes the problem of undetermined ki netic rate constants, solution method lb does not resolve a serious limitation of the FKA: it is bound to lead to very long computation times and numerical problems.The reason for this problem is that the transport and reaction processes that are included in the model occur on widely different time scales.Table 5 gives approximative values for the character istic time scale r for each process.
These characteristic time scales span several orders of magnitude, ranging from microseconds to days.This phe nomenon is called numerical stiffness (Boudreau, 1996b).Problems that are numerically stiff basically require special integration methods or rather small time steps in order to en sure accuracy.Effectively, the process with the smallest char acteristic time scale will set the pace of how the integration procedure progresses with time.Given the small time scales of the acid base reactions, pH models are very impractical or virtually impossible to solve with the FKA, even with in tegration methods that are specifically geared towards stiff problems (Chilakapati et al., 1998).A runtime comparison of all our presented approaches, including the FKA (solution method lb), is given at the end of the paper.
In conclusion, the FKA does not form a good choice for pH problems.As shown below, more refined alternatives to the FKA exist which do not depend on well constrained acid- base forward and backward rate constants2 and which dras tically reduce the computation time by reducing the stiffness of the system.

2.5
Step 5: Kinetic and equilibrium processes and species Table 5 shows that the characteristic time scales cluster in two groups.There is a group of comparatively slow pro cesses happening on a timescale of days (Processes 1 to 4) and a group of comparatively fast processes happening on timescales of fractions of a second to seconds (Processes 5 to 7).If the rate of one process is "sufficiently fast" compared to that of another process, this process can be assumed in lo cal equilibrium on the timescale of the "slower" process (e.g.Olander, I960; Aris and Mah, 1963;Otto and Quinn, 1971;DiToro, 1976;Saaltink et al., 1998).This allows to group the processes into slow kinetic processes, whose kinetics enter the model via suitable expressions, and fast equilibrium pro cesses, whose kinetics are neglected, i.e., local equilibrium is assumed to be reached instantaneously at any time.The designation of processes as "kinetic" or "equilibrium" also entails a corresponding classification of the species.K i netic species are those species whose concentrations are ex clusively influenced by kinetic processes, while equilibrium species are species that take part in at least one equilibrium reaction.
The grouping in kinetic and equilibrium processes de pends on the minimal time resolution of the model simu lations (cf.Saaltink et al., 1998) That means, the reference time to which to compare processes as "fast" and "slow" is the integration timestep of the model.In our model simula tions, the goal is to examine the pH changes over a period of days to weeks, resulting in an integration timestep of about one minute.
Accordingly, we assume the kinetics of reactions whose characteristic time scales are less than one minute to be neg ligible.Table 6 provides the resulting grouping of processes and species.
Note that local equilibrium is an assumption that is com monly used for systems to which our pH modelling approach Table 7. Kinetic process formulations.[X]up represents the up stream concentration of species X. [XHown its downstream concen tration.and [X]sat its saturation concentration.

Rox
is geared: macroscopic aquatic systems that contain reac tions on the timescale of several days as well as fast acid-base reactions.For models on this temporal and spatial scale, as suming local equilibrium does not change the model results.
However, the assignment of a process to the kinetic or equilibrium group is not absolute: it depends on the model time scale, and hence, on the questions addressed by the model.In our case, the exchange with the atmosphere is catalogued as a kinetic process.However, in a model that describes the pH evolution in the ocean over a million year time scale, exchange with the atmosphere can be considered an equilibrium process.Similarly, the dissociation Reac tions ( 5) to ( 7) are classified as equilibrium processes in our model.Yet, in a model that focuses on the fast relaxation of intracellular pH (model time scale of fractions of seconds), these same Reactions ( 5) to ( 7) would be considered kinetic reactions3.
The Processes (1) to (4) from Table 5 are modeled kinetically, and hence, we need to provide suitable constitutive expressions for their process rates .We describe oxic respi ration and nitrification as first order processes with respect to [OM] and [NH^~] respectively, and with a Monod depen dency on [ 0 2].The advective-dispersive transport is sim ply modeled as an exchange across the upstream and down stream system boundaries.The exchange with the atmo sphere is described by the classical reaeration mechanism (Thomann and Mueller, 1987).Table 7 lists the resulting kinetic expressions, parameters are given and explained in Table 14.

2.6
Step 6 : Mathematical closure of the system -the mass action laws Although the kinetics of the processes considered to be in lo cal equilibrium are neglected, these reactions are still part of the model.If the concentrations of species on either side of the reaction equation changes, the equilibrium shifts accord ing to Te Chatelier's principle.However, the rate of this shift is not governed by the kinetics of the equilibrium reaction themselves since they are assumed to be infinitely fast (local equilibrium is reached instantaneously).The net equilibrium  lco; reaction rates R dls are non-zero quantities and only depend on the supply rates of reactants and products of the equilib rium reaction in question due to slow kinetic processes in the model.In solution method lb we approximated the local equilibrium assumption by calculating, e.g., R*p+ by equa tion 6 b using a very high value for k f .However, under the true local equilibrium assumption the equilibrium is reached "instantaneously" which means

NH+
The latter expression is the equilibrium mass action law (Morel and Hering, 1993).This multiplication of an infi nite quantity with 0 renders the net equilibrium reaction rates R dls into mathematical unknowns.Chemical reasoning tells us that they are finite quantities.As a result, the system in Table 4 becomes an underde termined system with 9 equations (the M C E's) and 12 un knowns (9 species concentrations and the 3 equilibrium re action rates).To solve this system, it has to be mathemati cally closed by adding the equilibrium mass-action laws of the equilibrium reactions as additional algebraic constraints.
Including the mass action laws results in a fully deter mined initial-value differential algebraic equation (DAE) www.biogeosciences.net/5/227/2008/  Biogeosciences, 5, 227-251, 2008 Table 9. Model transformed into the canonical form: a fully de termined implicit initial-value DAE system.The combined mass conservation equations obtained by this transformation are equiv alent to the result of a series of linear combinations of the MCEs from Table 4: (4)+( 5)+( 6); ( 7)+(8); and (5)+2-( 6)+( 8)-( 9). (3) algebraic constraints (AEs) : m ass-action law s ( 7) 0 system (Fabian et al., 2001) (Table 8 ).The structure of this DAE system can be generalized as: where t is time.The DAE system is split into two parts: a differential part containing differential equations (Eq.8 a), and an algebraic part containing equations with no differen tials (Eq.8 b).It also contains two types of variables: the variables y whose differentials ^ are present (differential variables -the species concentrations) and the variables z whose differentials are absent (algebraic variables -the un known equilibrium reaction rates Rnh+> Rcc>2 > and R h c o -) -The algebraic part of the DAE system (Eq.8 b) only contains the differential variables y, and not the algebraic variables z (the equilibrium reaction rates).As can be seen from Ta ble 8 , the DAE system is fully determined (12 equations for 1 2 unknowns).
To our knowledge, this DAE system cannot be numerically integrated in the above form, despite being fully determined.To this end, the DAE system has to be reformulated first.For example, the DASSL routine (Differential Algebraic System Solver: Petzold, 1982) can solve implicit DAE systems (with suitable initial conditions given) of the form: This means that, if one wants to use the DASSL solver, the DAE equations may contain differentials of more than one variable (i.e., implicit differential equations), but the whole equation system can no longer contain the algebraic variables z.In the next step, we will discuss a suitable transformation that brings the DAE system in Table 8 in a form that can be solved by DASSL.

2.7
Step 7: Reformulation 1: transformation into the canonical form The system can be brought into a DASSL-solvable form by means of a transformation into the canonical form as discussed in, e.g., DiToro (1976), Steef el and MacQuarrie (1996), Lichtner (1996), Saaltink et al. (1998), Chilakapati et al. (1998) and Meysman (2001), based on an idea put forward by Aris and Mah (1963).During this transforma tion, the unknown equilibrium reaction rates are eliminated from the system.In a system with nes equilibrium species and n ep equilibrium reactions, the nes differential MCEs of the equilibrium species are then replaced by nei=nes-nep combined MCEs which do no longer contain the unknown equilibrium reaction rates.Appendix B details this proce dure for our problem.In our case the transformation of the system into the canonical form results in the reformulated DAE system as given in Table 9, which contains 9 variables and 9 equations.Note that the transformation procedure also provides explicit expressions for the unknown net and R*®+ (see equilibrium reaction terms R co 2> Appendix B).These can be used as output variables in the model and are sometimes of interest.

Solution m ethod [2]:
The model in Table 9 can be directly solved with the differential algebraic system solver DASSL (Petzold, 1982).This approach is referred to as the Full Numerical Approach (FNA).Still, this full numerical approach is not the most elegant way to approach the pH calculation.The equation set is sup plied "as it is" to an external numerical solver routine, which then performs the number crunching.A further reformula tion explicitly takes advantage of the chemical structure of the pH problem, thus allowing for less demanding numerical methods.

Step 8 : Introduction of equilibrium invariants
The Eqs. (4) to (6 ) of Table 9 contain differentials of more than one species on their left-hand sides.This means the differential part of the DAE system is not explicit and can not be solved by common integration methods such as Euler integration.To obtain a single differential on the left-hand side, one can introduce composite variables -as done in Ta ble 10 for our model.These composite concentration vari ables are referred to as equilibrium invariants.The reason for  9.
this nomenclature is straightforward.The right hand sides of Eqs. ( 4) to ( 6) do no longer contain the equilibrium reaction rates, and as a consequence, the rate of change of the equi librium invariants is not influenced by the equilibrium reac tions.Chemically, these equilibrium invariants thus can be seen as quantities that are conservative or invariant with re spect to the equilibrium reactions.Note that the definition of the equilibrium invariants introduces ne¡= 3 new variables.
To keep the DAE system determined, the definitions of the equilibrium invariants have to be added.The equilibrium invariants are in fact familiar quantities.We immediately recognize A and B as the dissolved in organic carbon (DIC) and total ammonium concentrations, which are denoted [ J ] C 0 2] and [J^N H ^"]. The third equi librium invariant is termed total alkalinity [TA].Again it has a familiar form: it is a subset of the total alkalinity [TA] as de fined by Dickson (1981).Still a number of subtleties should be stressed: 1) In our approach, the definition of total alkalinity fol lows naturally from the transformation into the canoni cal form and the elimination of the equilibrium process rates.It is not postulated a priori like in many previous pH modeling procedures (e.g.Regnier et al., 1997: Luff et al., 2001: Jourabchi et al., 2005).
2) The alkalinity definition is linked to a particular choice of kinetic and equilibrium reactions.Accordingly, when the reaction set is modified, the alkalinity definition might change as well.Also, even when keeping the same reaction set but choosing a different model time scale, one could arrive at a different alkalinity definition.
3) In our transformation procedure into canonical form we deliberately select suitable row operations during the Gauss-Jordan elimination (Appendix B) such that we obtain a subset of Dickson's total alkalinity as an equi librium invariant.Sticking to this practice, modifica tions in the reaction set and in the model time scale, as mentioned above, might result in different subsets of Dickson's total alkalinity.However, if this practice is abandoned, also different related quantities like, for ex ample, the "sum of excess negative charge" as used by Soetaert et al. (2007) can be obtained.
4) The right-hand side of the [TA] equation (Eq.( 6) in Ta ble 9) does not contain the rate of any equilibrium re action.This immediately shows that the alkalinity is a true equilibrium invariant, i.e., [TA] is not influenced by  (Morel and Hering, 1993) with corresponding mole balance equations including their equivalence to our equilibrium invariants.
species components equilibrium reactions, even though all its constituents are affected by these reactions4.
5) The influence of kinetic processes on [TA] can be di rectly inferred from the right-hand side of the [TA] equation (Eq.( 6) in Table 9).This implies that one does not need to invoke the electroneutrality of the solution or the notion of "explicit conservative total alkalinity" as advocated by Wolf-Gladrow et al. (2007) to obtain the influences of kinetic processes on [TA].
Note that the concept of equilibrium invariants is based on ideas put forward by amongst others DiToro (1976), Boudreau (1987), Boudreau and Canfield (1988), Boudreau (1991), Boudreau and Canfield (1993), and Morel and Her ing (1993).Furthermore, as illustrated in Table 11, trans forming the system into the canonical form and introduc ing equilibrium invariants is a formal mathematical way of finding suitable components in the tableau notation of Morel and Hering (1993) including the corresponding mole balance equations.
Note also that the sum of the transport terms on the right-hand side of the M C E's for the equilibrium invariants (Eqs.(4-6) in Table 9) can be directly calculated from the concentrations of the equilibrium invariants if the transport formulation for all species is the same, i.e., there is no dif ferential transport in the model.Mathematically that means that the transport formulation needs to be distributive over the sum.In our model, for example, this is realized by assum ing the same bulk dispersion coefficient E ' for all chemical species.
4Similarly the temperature invariance of [TA] can be inferred, since on the timescale of an integration timestep, temperature only influences the acid-base equilibrium reactions (via the .S'*'s), and these do not influence [TA].invariants !. System refomulated in terms of equilibrium invariants: explicit ODEs and equilibrium species as functions of [H+ ] and equilibrium :s.

MCEs of kinetic species
(1) algebraic constraints (AEs) [TA] = 2.9 Step 9: Reformulation 2: Operator splitting The algebraic part of our DAE system now consists of the mass action relations (Eqs.(7-9) in Table 9) and the definitions of the equilibrium invariants (Table 10).These equations feature the equilibrium species.However, each of the equilibrium species concentrations ([C 0 2], [HCO^"], [CO2-] , [NH3] and [NH+]) can be readily expressed in terms of the proton concentration [H+ ] and the associated equilibrium invariants ( [ E C 0 2] and [ENH^j"])-Ap pendix C describes this reformulation of the algebraic part of the DAE system.As a result, we obtain a novel DAE system (Table 12) where both the DE part and the AE part are reformulated in terms of the equilibrium invariants.

Solution method [3a]:
Although it still can be directly solved with DASSL, the system given in Table 12 can be solved with less numerical effort using the Operator Splitting Approach (OSA).This two step approach decouples the DAE system into an ordinary differential equation (ODE) system describing the kinetic reactions and an algebraic equation (AE) system that governs the equilibrium part (Luff et al., 2001: Meysman, 2001).
At each time step, the differential equation system is nu merically integrated, e.g., with an Euler integration routine (Press et al., 1992), which provides values for the differen tial variables (kinetic species and equilibrium invariants) at the next time step.Subsequently, the algebraic equation sys tem is solved at each timestep using the values for the dif ferential variables provided by the numerical integration of the ODE system.Due to its nonlinearity in [H+ ], the alge braic equation system must be solved numerically (e.g.us ing the van Wijngaarden-Dekker-Brent or Newton-Raphson methods given by Press et al., 1992) to find the chemically meaningful root (ƒ ([H+ ])= 0) of the function: The classical OSA (solution method 3a) takes advantage of the specific structure of the model to solve it in a more ele gant fashion than the FNA using DASSL.Still it requires at each time step the iteration between a numerical integration solver and a numerical root-finding technique, which might be computationally demanding.
Solution method [3b]: Recently, a modified OSA has been proposed (Follows et al., 2006), which is computationally faster.Rather than solving Eq. ( 10) directly, it acknowledges that carbonate alkalinity ([CA]=[HC0 (("]+2 [C0 3 _ ]) con tributes most to total alkalinity.In our case, using the proton concentration of the previous timestep [H+ ]prev and modelled equilibrium invariants (here [TA] and NH]j"]), the modelled carbonate alkalinity can be estimated by: which allows a first guess for the [H+ ] at the current time step by analytically solving the quadratic equation: i jy't-jy't- This first guess for [H+ ] is then used to evaluate Eq. ( 10) and test if its root has been found (with sufficient accuracy).
If not, the first guess for [H+ ] is used to calculate a better estimate for [CA] and the procedure is iteratively repeated.Iteration is mostly not necessary for buffered systems.
Note that this method also works if there are more minor contribution terms to [TA] than in our simple example.Note further that this method is inspired by the classical pH calcu lation methods of Culberson (1980), who analytically solved a cubic equation for systems with total alkalinity consist ing of carbonate and borate alkalinity only, and Ben-Yaakov (1970), who iteratively solved an equation for [H+ ] by start ing with an initial guess and by subsequent uniform incre ment of [H+], Although this improved OSA approach (solution method 3b) is advantageous, it does not allow assessing the influ ences of modelled kinetic processes on the pH.A further re formulation of the system is possible, which avoids numeri cal root-finding as well as the iterative procedure according to Follows et al. (2006).This method allows for the assess ment of the influences of modelled kinetic processes (includ ing subsequent re-equilibration of the system) on the pH.

2.10
Step 10: Reformulation 3: Direct substitution The classical OSA needs a numerical root-finding procedure because the algebraic equation (AE) part is non-linear with respect to the unknown proton concentration [H+ ].There fore, if one could make [H+ ] a differential variable, its value would be known before the solution of the algebraic equa tion system.This way, the algebraic equation system could be solved analytically and the numerical root-finding proce dure would not be necessary.To achieve this goal, the dif ferential equation for [TA] in Table 12 should be substituted by a differential equation in [H+ ].Partially following the ideas developed by Jourabchi et al. (2005) and Soetaert et al. (2007), this can be done by starting with the total derivative of the equilibrium invariant [TA].
Equation ( 12) in Table 12 tells us, that if all the dissocia tion constants (K* 's) are constant, the equilibrium invariant [TA] can be written as a function of exclusively the proton concentration and the equilibrium invariants These variables are functions of time t.Consequently, the total derivative of [TA] can be written as The subscripts indicate which quantities are held constant upon differentiation, and the shorthand notation c = [J ^ CO2 ], w =[J]N H ]j"] and h= [H+ ] has been used.Eq. ( 14) can be readily solved for d^ ^, resulting in Equation ( 15) can replace the differential equation for [TA] in Table 12.Each of the quantities on the right hand-side of Eq. ( 15) is explicitly known.The time derivatives of the equilibrium invariants are given by expressions ( 4)-( 6) in Ta ble 12. Furthermore, Appendix D1 shows how the partial derivatives of total alkalinity can be analytically calculated.
Table 13 shows the reformulated D E 's/M CE's of the DAE system.The AE part is the same as given in Table 12 (ex cept for the equation for [TA] which is obsolete) and is not repeated.
The quantity is a central and important quantity for pH modelling, as it modulates the effect of changes in state variables on [H+ ].Soetaert et al. (2007) call a similar quantity the buffering capacity of the solution, and Frankignoulle (1994) refers to the inverse of a related quantity as the chemical buffer factor of the solution.
Solution method [4],The explicit ODE system in Table 13 can be numerically integrated.Subsequently, the AE system is used to analytically calculate the equilibrium concen trations for every timestep of the numerical integration.The resulting approach is referred to as the Direct Substitu tion Approach (DSA) (Saaltink et al., 1998: Meysman, 2001).
3 [H+] advantage is that it makes maximal use of the chemical struc ture of the pH problem, to gain understanding and insight and to reduce the numerical effort.However, depending on the application, the OSA improved according to Follows et al. (2006) might have about the same computational require ments.The second and major advantage is that Eq. ( 15) di rectly quantifies the influence of the various kinetic processes on [H+ ] and hence on pH.To show this, one can rearrange Eq. ( 15) (or rather the last equation in Table 13) to the form where the a coefficients and T can be calculated at each time step using the expressions given in Appendix D2.The a-coefficients are modulating factors that express the influ ence on pH for each of the four kinetic reactions/processes.Similarly, the term T lumps the influence of advectivedispersive transport processes on pH.
Splitting up a coefficients into process specific modulation factors and the buffering capacity of the solution, the influ ences of kinetic processes (except transport) on the rate of change of the proton concentration can be formalized as: where the ß coefficients represent the process specific mod ulation factors, which can also be found in Table D2 in Ap pendix D2.
The influence of transport on the rate of change of the pro ton concentration can be written as  This means that the influence of a modelled kinetic process (except transport) on the d^t ^ can be calculated by multiply ing the kinetic rate of the process in question by a modulat ing factor ß divided by the buffering capacity of the solution 9[H+] ' influence °f transport on , however, is an ex pression of the transport terms for the equilibrium invariants divided by, again, the buffering capacity of the solution.

Baseline simulation
In a first step, we performed a baseline steady state calcula tion for our model estuary with boundary conditions for the year 2004, which serves as a reference situation for the two perturbations scenarios outlined in the introduction.Table 14 provides an overview of the parameters and boundary condi tions that were used in this baseline simulation.
Using the set of parameter values in Table 14, the DSA approach (solution method 4) was implemented within the modeling environment FEMME (Soetaert et al., 2002).The FORTRAN model code can be obtained from the author or downloaded from the FEMME website: http://www.nioo.knaw.nl/ceme/femme/.
The upstream concentrations were used as initial condi tions, and a time-dependent simulation was performed until steady-state was reached.Table 15 compares the concentra tions in the baseline simulation with values averaged over the year 2004.There is a good agreement between measured and modeled values.Also, the steady state rates for oxic min eralisation (Rox= 2 .8/zmol-N kg-1 d-1 ) and nitrification (R"it=8.2/zmol-N kg-1 d-1 ) are in good agreement with values from Table 1.This correspondence between model and measurements was obtained without tuning of model pa rameters.This provides confidence that the baseline simula tion captures the essential features of the carbon and nitrogen   2004), and a Schmidt number for carbon dioxide at a temperature of 12 °C and a salinity of 5 from Wanninkhof (1992).rox has been obtained by dividing pelagic oxic mineralisation rates from Soetaert and Herman (1995b) by measured [O M ] values for 2004.rn¡t has been calculated in similar fashion using nitrification rates obtained from Andersson et al. (2006).[CC^lsat has been calculated according to a formula given in Weiss (1974) and atmospheric carbon dioxide levels from Borges et al. (2004).All dissociation constants are on the free hydrogen ion scale and for a temperature of T=\2 °C and salinity 5=5. 3 -1 n r s 1 (Heip, 1988 (Borges et al., 2004;Wanninkhof, 1992) First order oxic mineralisation rate rox 0 .1 d" 1 (Soetaert and Herman, 1995b) Mass balance closure was verified for carbon, ni trogen and oxygen.
As noted above, one of the major advantages of the DSA approach is that one can partition d^t    D2.
proton concentration) into contributions by different kinetic processes (Eq.16).At steady state, overall consumption of protons should match overall production.Figure 2 shows that in our baseline simulation, the pH steady state is dom inated by the interplay between oxic mineralisation, nitrifi cation and CO2 air-water exchange.Oxic mineralisation and nitrification respectively produce about 49% and 40% of the protons consumed by CO2 outgassing.The remaining 11% are the result of advective-dispersive transport T) • The NH 3 exchange with the atmosphere plays a negligible role, as it produces only 0.3% of the protons consumed by CO2 air-water exchange.

The influence of H2 O auto-dissociation and denitrifica tion
In the formulation of the model, we deliberately neglected the auto-dissociation of water (Rh2o) and denitrification (Rden) to keep the model analysis as simple as possible.A model including H2 O auto-dissociation does not show any differences in steady state results (Table 15, Fig. 2).Accord ingly, R h 20 can be safely omitted.

Three perturbation scenarios
In the perturbation scenarios, the baseline steady state values were imposed as initial conditions.
Scenario A: Decrease in the upstivam organic matter loading: It is estimated that the organic matter loading in the river Schelde will be halved by a new sewage-treatment plant for the city of Brussels, which became operational in 2007.To simulate the impact of this change, we started off from the baseline simulation (values for the year 2004), and de creased the upstream organic matter concentration [OM]up from 5 0 //m ol N kg-1 to 2 5 //m ol N kg-1 on the fifth day of a 40-day model run.Figure 3 shows the evolution of pH, [TA], [ £ C 0 2] and [O2] for this scenario.After about 35 days a new steady state is reached, in line with the 10 day response time scale of the dominant transport and reac tion processes (Table 5).The decrease in OM loading re duces the steady state concentration of organic matter [OM] by 38% (not shown), while oxygen levels increase by 10% and [J^ CO2 ] levels remain virtually unchanged (slight de crease by 0.3%).Note that the changes occur monotonically.This is different for the total alkalinity, which shows a slight "overshoot" response.TA decreases from 5928.9 // mol kg-1 to a minimum value of 5927.9 // mol kg-1 after 6 days, but then stabilizes at a higher level of 5928.1 // mol kg-1 .This dip in [TA] is explained by a different temporal response of the mineralization, nitrification and transport terms (Fig. 4a).The change in the upstream OM concentration leads to a sharp decline in [OM] in the system, causing R ox (which produces alkalinity) to drop sharply as well.The nitrifica tion rate Rnit (which consumes alkalinity) however drops less sharply.As a result, temporarily "too m uch" alkalinity TA is consumed, which results in the observed dip in [TA].Note that this discussion is only interesting in terms of model in ternals but is not relevant to the real system, since the changes in [TA] are near or below the measurement accuracy.Also note that the decrease in [TA] (0.8 ¡x mol kg-1 ) is much smaller than the corresponding decrease in the DIC concentration (16 n mol kg-1 ).This difference is due to the rising pH and the associated re-equilibration within the car bonate system.Although CO2 ] decreases, the CO2 sys tem dissociates more, due to the pH increase, increasing [CO3 -].Hence, alkalinity does not follow the decrease in CO2 ] to similar extent (see Table 16).
The new steady state pH of 7.734 is only 0.4% higher than the baseline pH of 7.705.Figure 5 shows that the abrupt decrease in organic matter loading has only a small influence on the pH steady state.The individual contributions of all processes decline, except for the small contribution of advective-diffusive transport.That means the model results are in agreement with intuitive expectations: In a system with less organic matter input, the influence of physical  processes on pH increases relative to that of biological processes.

Scenario B: Spill o f ammonium nitrate
Due to the presence of the harbour of Antwerp and the sur rounding chemical industry, there is potential danger of ship accidents and spills of chemicals into the Schelde estuary.As an example, we consider a spill of ten thousand tons of ammonium-nitrate fertilizer (NH^~ NO^"; in the molar ratio 1:1).Furthermore, we consider a slowly leaking ship, where the chemicals are released within a period of 1 0 days (be tween day 5 and 15 of the simulation).To model this release, we need to include an extra source term for ammonium and nitrate in the M C E's (cf.Table 4).
In a similar manner as for the other kinetic processes, one can derive the influence of the addition of substance X to the system Ax on the proton change , by including Ax in the M C E's from Table 4. W hereas the addition of nitrate has no effect on the pH, the contribution of A NH+ to d^t 1 has the form as given in Eq. ( 17), and results in The second row in Fig. 3 shows the profiles for pH, [TA], CO2 ] and [O2 ] for this scenario.Drastic perturbations in the geochemistry of the estuary are simulated during the 1 0 days of leakage, and during a small period of about 15 days afterwards.The leakage results in a distinct peak in [J^NH]j"] (not shown), with values rising by roughly 620% from 36 p mol kg-1 to 260 p mol kg-1 .This is accompanied by a peak in [NO^], rising by 130% from 340 p mol kg-1 to 778 p mol kg-1 , which is due to both the leakage ANO3 and increased nitrification.Total alkalinity and CO2 ] tem porarily drop by 4% and 1% respectively.Oxygen conditions drastically drop from 158 p mol kg-1 to hypoxic conditions at 43 p mol kg-1 , due to a short period of intense nitrifica tion.
pH levels drop by approximately two tenths of a unit from 7.71 to 7.49.Figure 5 shows that this is mainly due to an increase of nitrification, and that the contribution of A NH+ itself is negligible.After 10 days of leakage, dynamic pH equilibrium is almost re-installed, and the proton production of nitrification is compensated by the proton release due to outgassing of CO2 and from transport.The influence of oxic mineralisation on d d o e s not significantly change during the spill, compared to the dominant components.After 10 days, the end of the leakage imposes a new perturbation on the system.

Scenario C: Spill o f ammonia
In this sceneario, we investigate a similar ship accident, but now with a spill of ten thousand tons of ammonia (NH3).The leakage period is identical to the previous case, and the input term for ten thousand tons of ammonia within 1 0 days becomes The contribution of A nh3 to d^t ^ is Figure 3 shows the profiles for pH, [TA], [ J ] CO2 ] and [O2 ] for this scenario.Again a distinct peak in [J^ NH)j"] is ob served (not shown), with the baseline concentration rising by a factor of 37.This is again accompanied by a 50% increase in [NO^], which is now solely the result of increased nitri fication.Total alkalinity and [J^ CO2 ] temporarily rise by 20% and 1% respectively.Oxygen concentrations are again greatly reduced (by roughly 97%), now almost reaching full anoxia with a minimum of 5 p mol kg-1 .The oxic minerali sation rate is much lower than in the baseline-simulation due to low oxygen concentrations.The pH level increases by more than one pH unit from 7.71 to 8.78. Figure 5 shows that this is mainly due to the input of NH3 into the estuary by the leak.Nitrification ini tially counters the proton consumption of NH 3 input (via conversion to am monia), but this effect decreases drastically due to decreasing oxygen levels (cf. the initial steep spike in R nit shown in the right panel of Fig. 4).The effect of outgassing of ammonia E nh3 on ^ only becomes im portant towards the end of the 1 0 day spill period, when almost steady state conditions are reached.At this point, NH 3 outgassing and subsequent dissociation of ammonium (as the equilibrium state changes due to Le Chatelier's prin ciple) balance, together with the effects of nitrification and advective-dispersive transport, the proton consumption of A nh3 and subsequent association of ammonium.W hen the leakage is stopped, the system returns to the pre-leakage state within a matter of 15 days.There is however a dip in pH and alkalinity before baseline values are attained again.Immedi ately after the leakage stops, there is still a lot of NH|" in the system, which is further nitrified.The effects of CO2 out gassing and advective-dispersive transport (which changes  4 Writing down a MCE for all species whose concentrations are influenced by modeled processes.The system is now solvable with the full kinetic approach (FKA).
5 Partitioning the modeled process into kinetic and equilibrium processes according to their timescales and defining kinetic expressions for kinetic processes.
6 Mathematically closing the system by formulating the mass action laws of the equilibrium processes.
7 Transforming the system into the canonical form: reformu lating it into an implicit DAE system without any purely al gebraic variables.The system is now solvable with the full numerical approach (FNA).
Introducing the equilibrium invariants to convert the differen tial equations of the DAE into explicit ODEs.
9 Reformulating the algebraic part of the DAE to explicitly ex press all equilibrium species as functions of [H+] and equilib rium invariants.The system is now solvable with the operator splitting approach (OSA).
10 Reformulating the system according to the direct substitu tion approach (DSA): substitute the expression for by an expression for ^ to get rid of the AE systems nonlinearity in an unknown variable.The expression for ^ can be partitioned such that the influences of modeled kinetic processes onto ^ can be quantified.
sign again) compensate for the proton production associated with nitrification.However, this compensation occurs with a certain time lag, creating the dip in pH after the initial spike.
The net absolute values of proton consumption or produc tion of all processes decreases during the 1 0 day spill period due to an increase in the absolute value of the buffering ca pacity , which changes from -0.165 IO5 to -5.15 IO5 (Fig. 4).It can be noted that in our modelled pH range the absolute value of increases with increasing pH.This leads to the conclusion that the higher the pH of the system, the closer to zero the influences of processes on d .Along the course of the mathematical reformulations two approximations have been made.1) To make the transition from the FKA to the FNA (the transformation into the canonical form) the local equi librium assumption has been applied.As mentioned earlier, this approximation generally has no influence on the results of models of macroscopic systems.
2) To reformulate the system into a form solvable by the DSA, the K* 's of the system are assumed constant.This has been done for didactical reasons to keep the math ematical expressions simple.This, however, is no limi tation: variable K* 's can be integrated into the DSA as well.Although these approximations have been made, our four main solution methods are still different mathematical for mulations of the same model yielding the same results: both approximations can also be made from the very beginning.The local equilibrium assumption can be included into the FKA (solution method lb) and the K* 's can be assumed con stant for all approaches.W hat remains is a chain of mathe matical transformations, with no further approximations in volved.
As shown in Table 18 and in Fig. 6 , even if the local equi librium assumption and constant K* 's are applied to all ap proaches, there is a clear trade-off between reformulation effort and the numerical resources required.The more the pH problem is initially reformulated, the less computation time is spent on actual pH simulations afterwards.The re formulations transform the pH problem into a more elegant mathematical form, and only require a one-time investment during the model generation process.Accordingly, when do ing multiple simulations as in a sensitivity analysis, the ini tial time investment in reformulation is likely to pay off very rapidly.Although, in terms of computational performance, the improved OSA and the DSA are comparable, the DSA additionally allows for the quantification of the influences of kinetically modelled processes on the pH.These influences are calculated against the background of re-equilibration of the system due to a set of acid-base equilibria.
The DSA approach thus comes out as the most power ful procedure to tackle pH models.However, in a system where the dissociation constants (K * 's) cannot be assumed constant, the influence of temperature, salinity and pressure on the dissociation constants has to be incorporated into the DSA.This has been deliberately omitted from this paper for didactical reasons.

Comparison with previous approaches
Past pH modeling approaches can be equated to one of the four solution methods in Fig. 6 .The most basic approach, the Full Kinetic Approach (FKA) has only been implemented sporadically (Steefel and MacQuarrie, 1996;Zeebe, 2007), because of the numerical stiffness of the obtained equation systems (solution method la and lb), and the need to obtain parameters that are not very well constrained (forward and backward rates of acid-base reactions in solution method la).After one reformulation step termed the transformation into canonical form (DiToro, 1976;Lichtner, 1996;Steefel and MacQuarrie, 1996;Chilakapati et al., 1998;Saaltink et al., 1998;Meysman, 2001) based on an idea first put forward by Aris and Mah (1963), one can implement the Full Numeri cal Approach (FNA), which involves a direct numerical so lution of the resulting differential algebraic equation system.Steefel and MacQuarrie (1996) list a number of packages, in cluding DASSL (Petzold, 1982), capable of solving a system according to the FNA.Gehlen et al. (1999) applied this solu tion technique in a relatively simple pH problem (4 acid-base reactions) to study the distribution of stable carbon isotopes in the pore water of deep sea sediments.We are not aware of FNA applications with realistic "field-type" reaction sets (in cluding 10 or more acid-base reactions).In these situations, FNA simulations are expected to require significant compu tational resources.
The demanding computations of the FNA can be avoided by means of a second reformulation, via the introduction of equilibrium invariants.This reformulation allows uncou pling the differential and algebraic part of the DAE system and solving them independently.The resulting approach is termed operator splitting (OSA, steps 8 and 9).Regnier et al. (1997) used the OSA to model pH along an estuarine gradi ent, Marinelli and Boudreau (1996) used it to study the pH in irrigated anoxic coastal sediments, and Follows et al. (1996) used the OSA to investigate the carbonate system in the water column of the North Atlantic.Besides pointing out different varieties of the FNA, Chilakapati et al. (1998) also applied the OSA to simple groundwater problems.W hile not explic itly reformulating the system, Boudreau (1987), Boudreau and Canfield (1988), Boudreau (1991), Boudreau andCan field (1993), andBoudreau (1996a) (the CANDI model) used the notion of dividing the reaction set into kinetic reactions and equilibrium reactions.Imposed equilibrium invariants were then used to simulate steady state pH profiles of aquatic sediments.Therefore, these approaches can be viewed as predecessors of the OSA.Although equilibrium invariants were not explicitly defined, Wang and Van Cappellen (1996) (the STEADYSED1 model) uncoupled the DE and the AE part of the DAE system and solved them separately, mak ing their approach a quasi OSA.In a detailed methodological study on pH modeling, Luff et al. (2001) examined different variations of the OSA: three different possibilities for the al gebraic equation part of the DAE system were introduced.
As noted above, there are two major disadvantages asso ciated with the classical OSA approach (1) the equilibration step requires a numerical solution, which makes the OSA computationally intense, and (2) the OSA does not allow quantifying the influence of different processes on d^t ^.The numerical solution step can be eliminated using the im proved OSA put forward by Follows et al. (2006), but it still lacks the possibility of assessing influences of kinetically modelled processes on the pH.
The two problems of the OSA vanish after a third refor mulation, which leads to the Direct Substitution Approach (DSA).Therefore, we consider the DSA approach to be the most elegant and promising pH modeling procedure, espe cially if knowledge of the influences of modelled processes on the pH is desired.If this knowledge is not desired, the im proved OSA according to Follows et al. (2006) might be the method of choice, since the third reformulation of the system (step 10) is not necessary.In the DSA, the differential equa tion for total alkalinity is replaced by a differential equation for the proton concentration, which enables a direct analyt ical solution of the equilibration step.The most important advantage is that the change in [H+ ] can be partitioned into contributions by different processes, and hence, the influence of processes on pH can be directly assessed (as discussed fur ther below).
Although applying the DSA, Meysman et al. (2003) (the MEDIA modelling environment) did not make use of its ca pability of assessing influences of processes on the pH.
In recent years two other studies have employed DSATike approaches to assess influences of processes on the pH.Yet the way these methods were derived is not fully clear as the presentations are prone to internal inconsistencies.
The approach of Jourabchi et al. (2005) is situated some where between the DSA and the FNA.As a by-product in calculating stoichiometric coefficients for equilibrium species, Jourabchi et al. (2005) calculated a rate of change of protons over time for a given modelled process, starting from the total derivative of total alkalinity.However, these rates do not add up to a total rate of change since the effect of transport is not made explicit.Direct proton transport is even omitted as they remove the mass conservation equation for protons to cope with an overdetermined equation system, which introduces an error.Their equation system was sub sequently solved with a numerical solver that depended on steady state conditions of the system.This means dynamic pH simulations are not possible.Total quantities like total alkalinity were imposed and not consistently derived.Subse quently, [TA] was used in a way that in some points contra dicted D ickson's (Dickson, 1981) notion of [TA].Soetaert et al. (2007) also made a step towards a DSA, but fell short of deriving a total rate of change of protons.They needed to invoke several ad-hoc assumptions and concepts like the mean and total charge of postulated total quantities to arrive at formulae for the influence of modeled processes on the pH.These formulae did not add up to a total rate of change of protons over time, because no transport terms were included.This means that modeling the pH of a real system containing several processes at the same time was not possi ble.

Implicit assumptions
The subsequent reformulations of the system (Fig. 6 ) yield more insight into the physical, chemical, and mathematical structure of the pH problem.By delineating all steps of the model generation process explicitly, one achieves a high level of model transparency.Typically, past treatments do not ex plicitly list all assumptions and decisions made during the model generation process.This practice has led to the intro duction of unnecessary assumptions and constraints, as well as inconsistently employed concepts.
A first difference between our approach and past treat ments is that we do not need an a priori definition of alkalin ity.In other words, in our treatment, the way alkalinity is de fined in terms of the other chemical species follows directly from the model reformulation.As shown above, alkalinity is one of the equilibrium invariants (like total inorganic carbon and total ammonium).These equilibrium invariants emerge after the transformation of the pH problem into the canoni cal form and are equivalent to the mole balance equations of M orel's components (Morel and Hering, 1993) as well as to DiToro's reaction invariants (DiToro, 1976) of the system.Similar quantities appear in Boudreau (1987), Boudreau and Canfield (1988), Boudreau (1991), and Boudreau and Can field (1993).Equilibrium invariants, and hence alkalinity, are quantities that are conservative with respect to equilibrium reactions.The exact form of alkalinity depends on the cho sen set of equilibrium reactions, and hence, it is dependent on the chosen pH range of the model and the chosen time scale of the model.This practice ensures (1) consistency be tween the definition of total alkalinity and the model, and (2 ) the correct stoichiometric coefficients in the MCE for total alkalinity (cf.Eq. ( 6) in Table 12).
A second difference is that we do not need the assump tion of electroneutrality.Approaches like e.g.L uff's (Luff et al., 2001) charge balance approach, or the CANDI model (Boudreau, 1996a) implicitly assume electroneutrality of the solution.They use a measure of total charge (including con servative ions like Na+ ), which is assumed to be zero, to close their equation systems5.Although sometimes wrongly termed so (e.g.Boudreau, 1991;Follows et al., 1996Follows et al., , 2006), total alkalinity is not a charge balance, but a proton balance.It expresses the excess of proton equivalents (protons and proton donors) to proton acceptors (Dickson, 1981;Wolf-Gladrow et al., 2007).This means that if, for example, NO((" is assumed not to react with H+ in the pH range that is mod elled, the concentration of nitrate does not have any influence on total alkalinity, although it is an integral part of the total charge balance of the solution.By consistently using total alkalinity instead of a charge balance, concentrations of con servative ions do not enter the pH calculation. Furthermore, as mentioned above, our approach directly provides the stoichiometric coefficients for total alkalinity for all kinetic processes.Although these coefficients are not obvious from the definition of [TA], as all component concentrations are influenced by equilibrium reactions, our model generation procedure unambiguously provides them.To this end, a reformulation of the expression of [TA] into the "explicit conservative total alkalinity", which requires elec troneutrality of the solution, as put forward by Wolf-Gladrow et al. (2007), is not needed.

Assessing the influence of processes on pH
In a system that contains slow kinetic processes (slow ki netic biogeochemical reactions, but also transport) and fast equilibrium processes, the independent drivers of the system, that means the factors that determine the temporal evolution of the state of the system, are only the slow kinetic processes.W hen adopting local equilibrium, the net rates of the equilib rium processes become dependent on the rates of the kinetic processes.
Therefore it is of interest to assess the influences of the slow kinetic processes on the pH.To do so, one needs an explicit formulation for d w h i c h can be partitioned into explicit contributions by each of these different slow kinetic processes.
Two of our main solution approaches provide an explicit formulation for : the FKA and the DSA.However, the formulation for ^ as obtained by the FKA contains terms contributed by equilibrium reactions.These equilib rium terms implicitly contain the influences of all slow ki netic processes that influence the reactants and products of the equilibrium reaction in question.Therefore, the formula tion for as obtained by the FKA cannot be partitioned into explicit separate terms for the influences of all slow ki netic processes on the pH.
Exactly there lies the most important advantage of the DSA method as it provides an explicit partitioning of the for mulation for d^t ^ into the influences of all slow kinetic pro cesses on the pH against the background of buffering by an equilibrium reaction system.
Unlike Soetaert et al. (2007) and Jourabchi et al. (2005), in our DSA we obtain an explicit formulation for the contri bution of all kinetic processes to d^t ^, including transport.This enables a deeper understanding of how pH steady state is attained, and what processes exactly are responsible for a pH change upon disturbance of the system.This is clearly illustrated in our disturbance scenarios for a simple estuarine system.
Furthermore, the buffering capacity of the solution is identified as an important and central quantity, as it modu lates the influence of all processes on the pH.A process with the same rate, can have a different influence on the pH de pending on the state of the system, as represented by j • Our disturbance scenario C shows that it is possible that in certain circumstances, although process rates increase, the absolute values of influences of processes on can de crease, since the absolute value of increases due to an increased pH.Figuratively this can be explained by the fact that a higher pH means less free protons in solution.There fore, the amount of protons affected by a certain process is decreased, fjppj is a measure for this condition.

Conclusions
In the present paper, we systematically and consistently de rived a succession of methods to model pH, making every step of the model generation process explicit.The chemical structure of the model was used for sucessive reformulations until fast and elegant numerical solutions were possible.Ex isting pH modelling approaches were identified within this framework and advantages and drawbacks were pointed out.Definitions for summed quantities and the influence of all modelled processes on them where derived from the model.With the DSA the influence of modelled kinetic processes on the pH can be quantified.
A. F. Hofmann et al.: pH model construction in aquatic systems stoichiometric matrix for the equilibrium reactions, Vkin is the nkp>aies stoichiometric matrix for the influence of the ki netic reactions on the equilibrium species, Ru" is the vector of the kinetic reactions, and R eq is the vector of the equilib rium reactions.
The goal is to find a linear transformation P such that " d[C] " P X --= P X Vkin X Run at P X Veq X R eq (B2) P can be constructed by performing a Gauss-Jordan elimi nation applied to the matrix veq (By adequate selection of the row operations during the Gauss-Jordan elimination, a sub set of Dickson's total alkalinity (Dickson, 1981) as well as a subset of any other desired similar quantity like Soetaert's "sum of excess negative charges" (Soetaert et al., 2007) can be obtained as an equilibrium invariant).

Fig. 1 .
Fig. 1. (a): the example estuarine system: the model domain (the stretch of river between the blue lines) encompasses around 40 river kilometres; (b): The model domain is represented by a conceptual model scheme with biogeochemical processes; For explanations of symbols, see text.
mentioned chemical species.This finalizes the conceptual model formulation -see scheme in Fig. lb.2.3 Step 3: Constraining the model pH range -selection of acid-base reactions 4: A mass conservation equation (MCE) for each species Overall, our model set includes a set of n p = l processes en compassing 5 reactions (R ox, R nit, R ^+ , R c o 2 > R h c 0 -) a n d OM, 0 2, N O J, C 0 2, H C O J, C O |" , NH+ NH3, and H+ www.biogeosciences.net/5/227/2008/Biogeosciences, 5, 227-251, 2008 The DSA is the end result of three sequential reformulations of the pH problem.The DSA has two advantages.The first www.biogeosciences.net/5/227/2008/Biogeosciences, 5, 227-251, 2008Table 13.ODE part of the DAE system with direct substitution of ¿[TA] , rf[H+] o 2 + t h c o ^ + T c o l -+ e c o 2 + Y Rox d ÍE N H +]

E
Biogeosciences, 5, 227-251, 2008 www.biogeosciences.net/5/227/2008/ Boundary conditions o f the model domain'.Values for [J] CO2 ] have been obtained from Hellings et al. (2001).All other values are NIOO monitoring values for 2004, except for the values for [TA] which have been consistently calculated."NM 2004" refers to measured data from 2004 obtained by the NIOO monitoring program.
Fig. 5. Budget of 1 for the three model scenarios.The gray line indicates the net 1.
pH modeling in 10 steps 1 Formulation of the model question.2 Formulation of the conceptual model.3 Constraining the model pH range -selection of acid-base re actions.

Fig. 6 .
Fig. 6.The trade-off between numerical resource requirement and model reformulation effort.

Table 3 .
pH range adjusted set of acid-base reactions.

Table 4 .
Mass conservation equations (MCEs) for each chemical species.

Table 5 .
Characteristic time r of processes to be modeled.Values for and Rqq2 are obtained from Zeebe and Wolf-Gladrow

Table 6 .
Kinetic and equilibrium processes and species.n x denotes the number of respective species or processes.

Table 8 .
Zeebe and Wolf-Gladrow (2001) system.Note that the dis sociation constants used are stoichiometric constants (denoted by the star as superscript: in contrast to thermodynamic constants: seeZeebe and Wolf-Gladrow (2001)for a description of different dis sociation constants). rf[OM]

Table 10 .
Composite variables to create explicit ODEs in Table

Table 11 .
The model system written in tableau notation

Table 14 .
Characteristic parameters o f the model domain'.Kp has been calculated by using a CgOO value (piston velocity), normalized to a Schmidt number of 600 (the value for carbon dioxide in freshwater at 20°C), for the Schelde at Antwerp fromBorges et al. (

Table 15 .
Steady state baseline values compared with measured values from 2004 (NIOO monitoring data).All quantities except for pH have the unit // mol kg-1.

Table 16 .
The carbonate system before and after the change in up stream organic matter loading (all values in /+ mol kg-1 ).

Table 17 .
Summary of our pH modeling approach.

Table 18 .
CPU time in miliseconds for one model run of all men tioned solution approaches and scenarios.Values are averages of 1000 runs each.All approaches are integrated with DASSL to be comparable.The model output generated with the five ap proaches is identical.The FKA is implemented according to so lution method lb.The OSA (3a) has been implemented using the