A reaction based rives tream water quality model, Model development and numerical schemes

A reaction-based river/stream water quality model:Model development and numerical schemesFan Zhanga,*,Gour-Tsyh Yeh b ,Jack C.Parker c ,Philip M.Jardineaa Environmental Science Division,Oak Ridge National Laboratory,Oak Ridge,TN 37831,USAb Department of Civil and Environmental Engineering,University of Central Florida,Orlando,FL 32816,USA cDepartment of Civil and Environmental Engineering,University of Tennessee,Knoxville,TN 37996,USAReceived 24June 2007;received in revised form 15October 2007;accepted 16October 2007KEYWORDS Rivers;Streams;Water quality;Sediment transport;Reactive transport;Simulation models;Chemical reactionsSummaryThis paper presents the conceptual and mathematical development of anumerical model of sediment and reactive chemical transport in rivers and streams.The distribution of mobile suspended sediments and immobile bed sediments is controlled by hydrologic transport as well as erosion and deposition processes.The fate and transport of water quality constituents involving a variety of chemical and physical processes is mathematically described by a system of reaction equations for immobile constituents and advective–dispersive–reactive transport equations for mobile constituents.To cir-cumvent stiffness associated with equilibrium reactions,matrix decomposition is per-formed via Gauss–Jordan column reduction.After matrix decomposition,the system of water quality constituent reactive transport equations is transformed into a set of ther-modynamic equations representing equilibrium reactions and a set of transport equations involving no equilibrium reactions.The decoupling of equilibrium and kinetic reactions enables robust numerical integration of the partial differential equations (PDEs)for non-equilibrium-variables.Solving non-equilibrium-variable transport equations instead of individual water quality constituent transport equations also reduces the number of PDEs.A variety of numerical methods are investigated for solving the mixed differential and algebraic equations.Two verification examples are compared with analytical solu-tions to demonstrate the correctness of the code and to illustrate the importance of employing application-dependent numerical methods to solve specific problems.ª2007Elsevier B.V.All rights reserved.IntroductionIn the past two decades,due to the rapid development of computer technology,water quality models have become popular tools for assessment of sediment and pollutant0022-1694/$-see front matter ª2007Elsevier B.V.All rights reserved.doi:10.1016/j.jhydrol.2007.10.020*Corresponding author.Tel.:+18655747314;fax:+18655740765.E-mail address:zhangf@ (F.Zhang).Journal of Hydrology (2008)348,496–509a v a i l ab l e a t w w w.sc i e n c ed i re c t.c o mjou rnal homep age:www.elsevier.c om/locate/jhydro lNomenclatureA river/stream cross-sectional area[L2]BP concentration of bed precipitates[M/M]BS concentration of bed sediment[M/L2]CB concentration of particulate sorbed on to bed sediment[M/M]C i concentration of water quality constituent i[M/M]CIMW concentration of dissolved chemical in the immobile water phase[M/M]CMW concentration of dissolved chemical in the mo-bile water phase[M/M]CS concentration of particulate sorbed on to sus-pended sediment[M/MD n deposition rate of the n th sediment[M/L2/T] DMA n amount of locally available dry matter of n th sediment[M/L2]E0n erodibility of the n th sediment[M/L2/T]E e concentration of the e th equilibrium-variable[M/L3]E m e concentration of mobile part of the e th equilib-rium-variable[M/L3]E i concentration of the i th non-equilibrium-vari-able[M/L3]E m i concentration of mobile part of the i th non-equilibrium-variable[M/L3]h water depth[L]h b river/stream bed depth[L]K e equilibrium constant of e th equilibrium reaction K x dispersion coefficient[L2/T]L transport operator incorporating source terms M total number of water quality constituentsME asi non-bank external source of the i th non-equilib-rium-variable[M/L/T]ME isi source of the i th non-equilibrium-variable from subsurface[M/L/T]ME os1i overland source of the i th non-equilibrium-vari-able from river banks1[M/L/T]ME os2i overland source of the i th non-equilibrium-vari-able from river banks2[M/L/T]ME rsi rainfall source of the i th non-equilibrium-vari-able[M/L/T]M asi non-bank external source of water quality con-stituent i[M/L/T]M isi source of water quality constituent i from sub-surface[M/L/T]M im number of immobile water quality constituentsM os1i overland sources of water quality constituent i from river bank1[M/L/T]M os2i overland sources of water quality constituent i from river bank2[M/L/T]M rsi rainfall source of water quality constituent i[M/ L/T]M m number of mobile water quality constituents M n concentration of the n th bed sediment[M/L2]MS asn non-overland source of the n th suspended sedi-ment[M/L/T]MS os1n overland source of the n th suspended sedimentfrom river bank1[M/L/T]MS os2noverland source of the n th suspended sedimentfrom river bank2[M/L/T]n*unit outward vectorN i weighting function at the i th node with the sameorder as base functionN j base function at the j th nodeNR number of biogeochemical reactionsN E number of equilibrium reactions/-variablesN K number of kinetic reactionsN NE number of non-equilibrium-variablesN S total number of sediment size fractionsp time step numberP river/stream cross-sectional wetted perimeter[L]Q river/streamflow rate[L3/T]r k rate of the k th reaction[M/L3/T]r bebackward rate of e th equilibrium reaction[M/L3/T]r feforward rate of e th equilibrium reaction[M/L3/T]r i j NR production rate of water quality constituent idue to all NR reactions[M/L3/T]R i production rate of i th non-equilibrium-variabledue to reactions[M/L3/T]R n erosion rate of the n th sediment[M/L2/T]S n concentration of the n th suspended sediment[M/L3]SP concentration of suspended precipitate[M/M]SS concentration of suspended sediment[M/L3]t time[T]V c Dn critical friction velocities for the onset of depo-sition of the n th sediment[L/T]V c Rn critical friction velocities for the onset of ero-sion of the n th sediment[L/T]V s n settling velocity of the n th sediment[L/T]w1time-weighting factor1w2time-weighting factor2W i weighting function at the i th node with sameorder as base function or one order higherx distance along the river/stream direction[L]h b porosity of the bed sediment[L3/L3]l ik reactant reaction stoichiometry of the i th waterquality constituent in the k th reactionm ik product reaction stoichiometry of the i th waterquality constituent in the k th reactionq i density of the phase associated with water qual-ity constituent i[M/L3]q w density of column water[M/L3]q wb density of bed water[M/L3]s b bottom shear stress or the bottom friction stress[M/L/T2]s cD n critical shear stress for the deposition of the n thsediment[M/L/T2]s c Rn critical shear stress for the erosion of the n thsediment[M/L/T2]D t simulation time step size[T]A reaction-based river/stream water quality model:Model development and numerical schemes497distributions to facilitate efficient resource management. Research on river/stream water quality modeling has in-volved studies of both sediment transport(Engelhardt et al.,1995;Zeng and Beck,2003;Rathburn and Wohl, 2003)and chemical transport(Park and Lee,2002;Boor-man,2003;Lopes et al.,2004).Most existing surface water quality models simulate either specific systems(Cerco and Cole,1995;Park and Lee,2002;Lopes et al.,2004)or are limited to certain chemicals and/or reactions(Brown and Barnwell,1987;Ambrose et al.,1993;Park et al.,2003). These models may serve as efficient monitoring and man-agement tools when applied to the specific conditions forwhich they have been developed and calibrated,but have limited applicability to other environmental conditions. With better understanding and improved mathematical for-mulations of complex biogeochemical interactions(Tho-mann,1998;Somlyody et al.,1998;Mann,2000;Yeh et al.,2001),generic models capable of simulating user-prescribed reaction networks have been developed that have broader applicability(Steefel and Cappellen,1998).A reactive system is completely defined by specifying all relevant biogeochemical reactions(Yeh et al.,2001),such as aqueous complexation,sorption,precipitation/dissolu-tion,volatilization,diffusion,sedimentation,etc.In trans-port simulations,it is useful to classify reactions as either ‘‘slow’’(kinetic)or‘‘fast’’(equilibrium)(Rubin,1983).A few reaction-based watershed models can handle contami-nant transport with kinetic reactions(Yeh et al.,1998; Cheng et al.,2000)using a primitive approach,in which the partial differential equations(PDEs)governing reactive transport are integrated directly to yield distributions of water quality constituents in a region of interest over time. However,when some reactions exhibit very fast kinetics (near equilibrium behavior),the PDEs becomes infinitely stiff making solution using the primitive approach intracta-ble(Fang et al.,2003).Therefore,a mixed differential and algebraic equations(DAEs)approach was introduced to overcome this problem by introducing DAEs for combina-tions of water quality constituents according to equilibrium reaction expressions and eliminating equilibrium reactions from the simultaneous solution of PDEs(Zhang et al., 2007).Although a systematic diagonalization approach ap-plied to the non-discrete system of PDEs involving decompo-sition by pivoting on equilibrium reactions and decoupling them from kinetic reactions has been developed for water quality modeling in subsurface porous media(Kra¨utle and Knabner,2005;Kra¨utle and Knabner,2007;Zhang et al., 2007),most surface water DAE models(DiToro,1976;Am-brose et al.,1993;Cerco and Cole,1995)have required manual identification of the DAEs.No existing surface water quality model,to our knowledge,has used a fully mechanis-tic approach to model chemical transport subject to both ki-netic and equilibrium reactions in river/streams.This paper presents a reaction-based numerical model simulating sediment and reactive chemical transport in riv-er/streams involving both kinetic and equilibrium reactions. In the model,sediments are categorized based on their physical and chemical properties.For each category of sed-iment,we include mobile suspended sediment particles within the water column and immobile sediment in the stream bed that is not subject to longitudinal transport. For chemical species,we define six phases and three forms.The six phases are suspended sediment,bed sediment,mo-bile water,immobile water,suspended precipitate,and bed precipitate phases,and the three forms are dissolved chem-icals,chemicals sorbed on sediment,and precipitates (Fig.1).Usually,water quality constituents of chemical spe-cies in the suspended sediment phase,the mobile water phase and the suspended precipitate phase are considered mobile.Constituents in the bed sediment phase,the immo-bile water phase and the bed precipitate phase are consid-ered immobile.Five numerical options are investigated to solve trans-port equations and three coupling strategies are considered to deal with reactive chemistry.For a given spatial and tem-poral discretization,some of these numerical options and coupling strategies are more accurate,while others sacri-fice accuracy to reduce CPU time.The main objective of this paper is to illustrate the applicability of various numer-ical options and coupling strategies to different types of problems for special application circumstances. Mathematical basisBed sedimentsThe balance equation for bed sediments stipulates that the rate of mass change is due to deposition/erosion asoðPM nÞo t¼PðD nÀR nÞ;n2N S;ð1Þwhere P is the river/stream cross-sectional wetted perime-ter[L],M n is wetted perimeter-averaged concentration of the n th bed sediment in mass per unit bed area[M/L2],t is the time[T],D n is the deposition rate of the n th sediment in mass per unit bed area per unit time[M/L2/T],R n is the erosion rate of the n th sediment in mass per unit bed area per unit time[M/L2/T],n2N S is used short for n2{1,...,N S},and N S is the total number of sediment size fractions.Concentrations of all bed sediments must be gi-ven initially for transient simulations.No boundary condi-tions are needed for bed sediments since they are regarded as immobile.The mechanics of sediment transport for cohesive and non-cohesive materials are different(Yang,1996).In gen-eral,sediment sizes smaller than2l m are generally consid-ered cohesive sediment;sediment of size greater than 60l m is coarse non-cohesive sediment;and silt(2–60l m) is considered to be between cohesive and non-cohesive sed-iment and usually determined with site-specific definition (Huang et al.,2006).We describe sediment deposition and erosion rates in Eq.(1)using different equationsfor Figure1Sediments and chemicals in river/streams.498 F.Zhang et al.cohesive sediments(Yeh et al.,1998;Gerritsen et al., 2000),such as Eqs.(2)and(3)D n¼minðV s n S n P Dn;S n h=D tÞwhere P Dn¼maxð0;1Às b=s c DnÞð2ÞR n¼minðE0n P Rn;DMA n=D tÞwhere P Rn¼maxð0;s b=s c RnÀ1Þð3Þand non-cohesive sediments(Yeh et al.,1998;Prandle et al.,2000),such as Eqs.(4)and(5)D n¼minðV s n S n N Dn;S n h=D tÞwhereN Dn¼max½0;1ÀðV c Dn=V c RnÞ2 ð4ÞR n¼minðE0n N Rn;DMA n=D tÞwhereN Rn¼maxð0;V c Dn=V c RnÀ1Þð5Þwhere V sn is the settling velocity of the n th sediment class [L/T],S n is the concentration of the n th suspended sedi-ment[M/L3],h is the water depth[L],D t is the simulation time step size[T],s b is the bottom shear stress or the bot-tom friction stress[M/L/T2],s c Dn is the critical shear stress for the deposition of the n th sediment[M/L/T2],E0n is the erodibility of the n th sediment[M/L2/T],DMA n is the amount of locally available dry matter of n th sediment,ex-pressed as dry weight per unit area[M/L2],s c Rn is the criti-cal shear stress for the erosion of the n th sediment[M/L/ T2],V c Dn and V c Rn represent the critical friction velocities for the onset of deposition and erosion,respectively[L/ T].It should be noted that any other phenomenological equations to estimate sediment deposition and erosion rate can be easily incorporated in the computer code. Suspended sedimentsThe continuity equation of suspended sediment can be de-rived based on the mass conservation law asoðAS nÞo t þoðQS nÞo xÀoo xAK xo S no x;¼MS os1n þMS os2nþMS asnþðR nÀD nÞP;n2N sð6Þwhere A is the river/stream cross-sectional area[L2],S n is the cross-sectional-averaged concentration of the n th sus-pended sediment in the unit of mass per unit column volume [M/L3],Q is the river/streamflow rate[L3/T],x is distance along the river/streamflow direction[L],K x is the disper-sion coefficient[L2/T],MS os1n and MS os2nare overland sourcesof the n th suspended sediment from river bank1and2,respectively[M/L/T],and MS asn is a non-overland source ofthe n th suspended sediment[M/L/T].Initial concentrations of all suspended sediments must be given for transient simulations.Four types of boundary con-ditions(Yeh et al.,1998)for suspended sediments are taken into account that may be stated as follows.Dirichlet boundary conditionThis condition is applied when concentration is given at the boundary asS n¼S nðx b;tÞ;ð7Þwhere x b is the location of the boundary node[L],and S n(x b,t)is a time-dependent concentration at the boundary [M/L3].Cauchy/Robin boundary conditionThis boundary condition is employed when the total mate-rialflow rate is given at the boundary asn*QS nÀAK xo S no x¼Q Snðx b;tÞ;ð8Þwhere n*is a unit outward vector,Q Sn(x b,t)is a time-depen-dent materialflow rate at the boundary[M/T].Usually,this boundary is an upstream boundary.Neumann boundary conditionThis boundary condition is used when the diffusive material flow rate is known at the boundary asÀn*AK xo S no x¼Q Snðx b;tÞ;ð9Þwhere Q Sn(x b,t)is a time-dependent diffusive materialflow rate at the boundary[M/T].Usually,this boundary is a downstream boundary.Variable boundary conditionThis boundary condition is utilized when the boundary loca-tion(upstream or downstream)is not predetermined or the flow direction changes with time during simulations.It re-fers to‘‘Neumann’’at the outflow part of the domain ðn*Q>0Þwith Q Sn(x b,t)=0and to‘‘Cauchy/Robin’’at the inflow part of the domainðn*Q<0Þwith Q Snðx b;tÞ¼n*QS nðx b;tÞ,where S n(x b,t)is a time-dependent concentration at the boundary that is associated with the incomingflow[M/L3].Immobile water quality constituentsThe balance equation for immobile water quality constitu-ents is simply the statement that the rate of mass change is due to biogeochemical reactions asoðPÁh bÁq wbÁh bÁCIMWÞo t¼PÁh bÁr CIMW jNR0;ð10ÞoðPÁh bÁq wbÁh bÁBPÞ¼PÁh bÁr BP jNR0;ð11ÞoðPÁCBÁBSÞo t¼PÁh bÁr CB jNR0;ð12Þwhere h b is the river/stream bed depth[L],q wb is the den-sity of bed water[M/L3],h b is the porosity of the bed sedi-ment[L3/L3],CIMW is the concentration of dissolved chemical in the immobile water phase in units of chemicalmass per bed-water mass[M/M],r CIMW jNR0is the production rate of CIMW due to all NR reactions in units of chemical mass per bed volume per time[M/L3/T],BP is the concen-tration of bed precipitate in units of chemical mass perbed-water mass[M/M],r BP jNR0is the production rate of BP due to all NR reactions as chemical mass per bed volume per time[M/L3/T],CB is the concentration of particulate sorbed on to bed sediment as chemical mass per bed sedi-ment mass[M/M],BS is the concentration of bed sedimentas sediment mass per bed area[M/L2],r CB jNR0is the produc-tion rate of CB due to all NR reactions as chemical mass per bed volume per time[M/L3/T].Defining r i jNR¼PÁh bÁr i jNR0=A for i=CIMW,BP and CB, Eqs.(10)–(12)can be modified asA reaction-based river/stream water quality model:Model development and numerical schemes499oðPÁh bÁq wbÁh bÁCIMWÞo t ¼AÁr CIMW jNR;ð13ÞoðPÁh bÁq wbÁh bÁBPÞo t ¼AÁr BP jNR;ð14ÞoðPÁCBÁBSÞo t ¼AÁr CB jNR:ð15ÞConsidering more than one variable CIMW,BP,BS, respectively,and defining q i=PÆh bÆq wÆh b/A for CIMWs and BPs and q i=PÆBS/A for CBs,Eqs.(13)–(15)can be sum-marized asoðA q i C iÞo t ¼Ar i jNR;i2M im;ð16Þwhere C i is the concentration of water quality constituent i, which is immobile,in units of chemical mass per unit phase mass[M/M],q i is the density of the phase associated with water quality constituent i[M/L3],r i j NR is the production rate of water quality constituent i due to all NR reactions in chemical mass per column volume per time[M/L3/T], and M im is the number of immobile water quality constitu-ents.Concentrations of all immobile water quality constitu-ents must be given initially for transient simulations.No boundary conditions are needed for immobile water quality constituents.Mobile water quality constituentsThe continuity equation of mobile water quality constitu-ents can be derived based on the conservation law of mate-rial mass stating that the rate of mass change is due to both advective-dispersive transport and biogeochemical reac-tions asoðAÁq wÁCMWÞo t þLðq wÁCMWÞ¼Ar CMW jNR;ð17ÞoðAÁq wÁSPÞo t þLðq wÁSPÞ¼Ar SP jNR;ð18ÞoðAÁCSÁSSÞþLðCSÁSSÞ¼Ar CS jNR;ð19Þwhere q w is the density of column water[M/L3],CMW is the concentration of dissolved chemical in the mobile water phase in units of chemical mass per column-water mass [M/M],r CMW j NR is the production rate of CMW due to all NR reactions in units of chemical mass per column volume per time[M/L3/T],SP is the concentration of suspended precipitate in chemical mass per column-water mass[M/ M],r SP j NR is the production rate of SP due to all NR reactions given as chemical mass per column volume per time[M/L3/ T],CS is the concentration of particulate sorbed on to sus-pended sediment in units of chemical mass per unit of sed-iment mass[M/M],SS is the concentration of suspended sediment as sediment mass per volume[M/L3],r CS j NR is the production rate of CS due to all NR reactions in units of chemical mass per volume per time[M/L3/T].Defining q i=q w for CMWs and SPs and q i=SS for CSs,Eqs.(17)–(19)can be summarized asoðA q i C iÞo t þLðq i C iÞ¼Ar i jNR;i2M m;ð20Þwhere C i is the concentration of water quality constituent i, which is mobile,in units of chemical mass per unit phase mass[M/M],q i is the density of the phase associated with water quality constituent i[M/L3],r i j NR is the production rate of water quality constituent i due to all NR reactions in chemical mass per volume per time[M/L3/T],M m is the number of mobile water quality constituents,and transport operator L(incorporating source terms)is defined asLðq i C iÞ¼oðQ q i C iÞo xÀoo x½AK xoðq i C iÞo xÀðM asiþM rsiþM os1iþM os2iþM isiÞ;ð21Þwhere M asiis the non-bank external source of water qualityconstituent i[M/L/T],M rsiis the rainfall source of waterquality constituent i[M/L/T],M os1iand M os2iare the overlandsources of water quality constituent i from river bank1and2,respectively[M/L/T],and M isiis the source of water qual-ity constituent i from subsurface[M/L/T].Concentrationsof all mobile water quality constituents must be given ini-tially for transient simulations.Similar to suspended sedi-ment transport,four types of boundary conditions aretaken into account for mobile water quality constituents,including Dirichlet,Variable,Cauchy/Robin,and Neumannboundary conditions.Decomposition of water quality constituent reactive transport systemsFrom a mathematical point of view,the temporal–spatialdistribution of water quality constituents is described by asystem of M im reaction equations(Eq.(16)),and M m reac-tive transport equations(Eq.(20)).These two equationscan be recast in the formoðA q i C iÞþa i Lðq i C iÞ¼Ar i jNR;i2M¼M imþM m;ð22Þwhere a i is0for immobile constituents and1for mobile con-stituents,and M is the total number of water quality constituents.The determination of r i j NR and associated parameters is a primary challenge in biogeochemical modeling.Instead of using adhoc methods to formulate r i j NR as empirical rates of water quality constituent concentration change due to all lumped reaction rates,we employ reaction-based formu-lations(Steefel and Cappellen,1998).In a reaction-based formulation,r i j NR[M/L3/T]is given by the summation of rates of all individual reactions that the i th water quality constituent participates,r i jNR¼X NRk¼1½ðm ikÀl ikÞr k :ð23ÞThis results in the transport equations of M water quality constituents described byoðA q i C iÞo tþa i Lðq i C iÞ¼AX NRk¼1½ðm ikÀl ikÞr k ;i2M;or Uo C Ao tþa LðCÞ¼A m r;ð24Þwhere m ik is the reaction stoichiometry of the i th water qual-ity constituent in the k th reaction associated with the prod-ucts,l ik is the reaction stoichiometry of the i th water quality constituent in the k th reaction associated with the500 F.Zhang et al.reactants,r k is the rate of the k th reaction[M/L3/T],U is a unit matrix,a is a diagonal matrix with a i as its diagonal component,m is the reaction stoichiometry matrix,and r is the reaction rate vector with NR reaction rates as its com-ponents.For the sake of simplicity,C A is a vector with com-ponents that represent M water quality constituent concentrations incorporating q i multiplied by the river cross-section area,and C is a vector with components rep-resenting M water quality constituent concentrations incor-porating q i.Eq.(24)is a representation of mass balance for any water quality constituent i in a reactive transport sys-tem,which states that the rate of change of any water qual-ity constituent mass is due to advection,dispersion and reactions that describe biogeochemical processes.In a primitive approach,Eq.(24)is integrated to yield the distributions and evolutions of water quality constituents in a region of interest.However,when some equilibrium reac-tions taking place in the system,this approach is not ade-quate(Fang et al.,2003).Here,we will take a diagonalization approach through decomposition.Eq.(24) written in matrix form can be decomposed based on the type of biogeochemical reactions via Gauss–Jordan column reduction of reaction matrix m(Chilakapati,1995).Among all the equilibrium and kinetic reactions,‘‘redundant reac-tions’’are defined as equilibrium reactions that can be de-rived from other equilibrium reactions.In order to avoid the singularity of the reaction matrix,redundant equilibrium reactions are automatically removed from the system prior to decomposition,if users inadvertently include them. Decomposition through Gauss–Jordan column reduction of the reaction matrix m is then performed by pivoting on the N E equilibrium reactions and automatically separating them from the N K kinetic reactions asA101 A2U1!o C A1d to C A2()þB101B2a1!LC1C2&'¼A D1K102K2!r1r2&';ð25Þwhere A1and A2are submatrixes of the reduced U matrix with size of N E·N E and(MÀN E)·N E respectively,01is zero submatrix of the reduced U and a matrixes with size of N E·(MÀN E),U1is unit submatrix of the reduced U ma-trix with size of(MÀN E)·MÀN E),C A1and C A2are subvec-tors of the vector C A with size of N E and MÀN E respectively,B1and B2are submatrixes of the reduced a matrix with size of N E·N E and(MÀN E)·N E respectively, a1is a diagonal submatrix of the reduced a matrix with size of(MÀN E)·(MÀN E),C1and C2are subvectors of the vec-tor C with size of N E and MÀN E respectively,D1is the diag-onal matrix representing a submatrix of the reduced m with size of N E·N E reflecting N E linearly independent equilib-rium reactions,K1and K2are submatrixes of the reduced m with size of N E·N K and(MÀN E)·N K respectively, reflecting the effects of N K kinetic reactions,02is a zero matrix representing a submatrix of the reduced m matrix with size of(MÀN E)·N E,and,r1and r2are subvectors of the equilibrium and the kinetic reaction rate vector r with size of N E and N K respectively.As a result,the decomposition of Eq.(24)to Eq.(25) effectively reduces a set of M simultaneous reactive trans-port equations into two sets of equations.Thefirst set con-tains N E equilibrium-variable transport equations involving equilibrium reactions and the second set contains N NE=MÀN E non-equilibrium-variable transport equations involving no equilibrium reactions.These two sets are de-scribed below.Thermodynamic equilibrium equations for equilibrium reactionsoðAE eÞþLðE meÞ¼AD1ee r1eþAX N Kj¼1K1ej r2j¼AD1eeðr f1eÀr b1eÞþAX N Kj¼1K1e j r2j;e2N E;e2N E)r f1e!1and r b1e!1)9a thermodynamically consistent equation:K e¼Yj2MA m j ej,Yj2MA l j ejor F eðC1;...;C M;p1;p2;...Þ¼0;ð26Þwhere E e¼P N Ej¼1A1ej C1j or E=A1C1and E me¼P N Ej¼1B1ej C1j or E m=B1C1.E e is an equilibrium-variable,which is a combina-tion of water quality constituent concentrations with reac-tion term involving equilibrium reactions;r feand r beare the forward and backward rates[M/L3/T],respectively,.Since all the equilibrium reactions in the system are linearly inde-pendent after redundant reactions are removed,one and only one equilibrium reaction appears in every one of the equilibrium-variable transport equations.For reversible reactions that are fast,equilibrium may be regarded as beingreached instantaneously(r f1e!1and r b1e!1)among all the relevant water quality constituents to assure a thermo-dynamically consistent expression.Therefore,E e is consid-ered to reach its equilibrium value instantaneously and an equilibrium-variable transport equation is replaced by the thermodynamically equilibrium equation representing the corresponding equilibrium reaction.In Eq.(26),K e¼P j2M A m j ej=P j2M A l j ejrepresents a mass action equilibrium equation where K e is the equilibrium constant of e th equilib-rium reaction and A j is the activity of j th water quality constituents.F e(C1,...,C M;p1,p2,...)=0symbolizes a user-defined equilibrium equation where F e is a nonlinear algebraic function of water quality constituent concentra-tions(C1,...,C M)and a number of parameters(p1,p2,...).Reactive transport equations for non-equilibrium-variablesoðAE iÞo tþLðE miÞ¼AX N Kj¼1K2ij r2j;i2N NE;where E i¼X N Ej¼1A2ij C1jþC2i or E¼½A2U1C1C2&';and E mi¼X N Ej¼1B2ij C1jþa1i C2i or E m¼½B2a1C1C2&';ð27ÞA reaction-based river/stream water quality model:Model development and numerical schemes501。