Application of a multi-step computation scheme to an intratidal estuarine water quality model

Ecological Modelling 110 (1998) 281 – 292

Application of a multi-step computation scheme to an intratidal estuarine water quality model Kyeong Park a,*, Jian Shen b, Albert Y. Kuo b b

a Department of Oceanography, Inha Uni6ersity, Inchon, 402 -751, South Korea School of Marine Science/Virginia Institute of Marine Science, The College of William and Mary, Virginia, 23062, USA

Received 18 November 1997; accepted 30 March 1998

Abstract A computation scheme, which involves multi-step computation of the decoupled mass-balance equations in water quality models, was proposed by Park and Kuo (1996. Water Res. 30, 2255 – 2264). Using a simple hypothetical model, they explained the concept and demonstrated the accuracy and computational efficiency of the multi-step scheme. This paper describes the concept and the implementation of the multi-step computation scheme in a more generalized format for both intratidal and intertidal estuarine water quality models. The scheme is applied to a prototype estuary using a laterally-integrated two-dimensional intratidal model. The computational efficiency of the scheme is demonstrated by comparing the model results using the multi-step scheme with those using the traditional single-step explicit numerical solution of the kinetic processes. © 1998 Elsevier Science B.V. All rights reserved. Keywords: Multi-step computation scheme; Mass-balance equation; Estuarine water quality model; Intertidal model; Intratidal model

1. Introduction Water quality management in estuarine and coastal waters has received increased attention in recent years as human activities in these areas increase. Mathematical models based on physical * Corresponding author. Tel.: + 82 32 8607707; fax: + 82 32 8620988; e-mail: [email protected]

and biogeochemical principles have been extensively used to study the water quality conditions and to provide consistent, rational tools for water quality management in estuarine and coastal waters. Most mechanistic water quality models simulate the spatial and temporal distributions of water quality state variables by solving mass-balance equations, which in the general form may be expressed as:

0304-3800/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PII S0304-3800(98)00082-9

K. Park et al. / Ecological Modelling 110 (1998) 281–292

282

( (VC) (t = [physical transport processes] +[biogeochemical processes],

(1)

with V, volume of any given portion of a water body; C, average concentration of a dissolved or a suspended substance in the given volume; and t, time. Eq. (1) states that the time rate of change of mass within the given volume is determined by physical transport and biogeochemical processes. The time scale of the local derivative term in Eq. (1) can be either intratidal (real-time) or intertidal (tidal average), and the spatial scale is determined by the size of the given volume V. Physical transport processes in Eq. (1), consisting of advective and diffusive transports, are presumed to be identical for all water quality state variables in most water quality models. The time scale can be intratidal (minutes to an hour) or intertidal (tidal cycle to days), and the spatial scale can be zero-, one-, two- or three-dimensional. We often obtain the physical transport processes by solving a hydrodynamic model, which typically consists of continuity, momentum and salt-balance (sometimes heat-balance) equations. Biogeochemical processes in Eq. (1) consist of external and internal sources and sinks, the latter being primarily due to the kinetic processes for biogeochemical transformations. Different water quality state variables are affected by different kinetic processes. The fastest kinetic process in eutrophication models is usually the growth of algae, which has the rate coefficient of about 2 day − 1, and thus the biogeochemical processes have the time scale of the order of hours. As far as the numerical computation is concerned, biogeochemical processes have no spatial scale except the settling of particulate matters, since it is physical transport processes that move materials spatially. Materials go through biogeochemical transformations while being physically transported. Most numerical models have paid more attention to the solution scheme for the physical transport terms in Eq. (1) than to the kinetic terms. They have routinely treated the kinetic terms explicitly, which is first-order accurate, in order not

to introduce additional complications into the solution scheme, resulting in inconsistent treatment of physical transport and kinetic processes. Shanahan and Harleman (1982) have pointed out the necessity in numerical modeling to consider both hydrodynamics and biogeochemistry in a compatible and even-handed fashion. Intratidal real-time models must advance the computation of Eq. (1) with a time-step of the order of seconds to minutes, thus updating the kinetic processes unnecessarily too often. Intertidal tidal-average models must have a time-step of the order of hours because of the time scales of the kinetic processes, and thus update the physical transport processes, which can have a time-step of one tidal cycle, unnecessarily too often. Hence, differences in the time scales of the physical transport and kinetic processes impose computational inefficiency in numerical solutions of Eq. (1). A great deal of inconsistency and computational inefficiency may be avoided in the numerical solution of Eq. (1) if the physical transport and the kinetic processes are decoupled. The decoupling method has been employed in the water quality modeling of surface water (Ambrose et al., 1988; Kuo and Neilson, 1988) and groundwater (Valocchi and Malmstead, 1992). The solution schemes of these models involved two-step computation, in which substances are physically transported and then followed by the application of kinetic processes. The decoupling of Eq. (1) not only simplified the solution scheme but also made the models more flexible with respect to the addition of new state variables and to the modification of the kinetic formulations. Park and Kuo (1996) proposed a new computation scheme that also involves the decoupling of Eq. (1). Their proposed scheme solves the decoupled equations by employing multi-step computation, in which the kinetic processes are applied for the first half time-step, followed by the physical transport for a full time-step, and the kinetic processes are applied again for the remaining second half time-step. The kinetic equations may be solved numerically or analytically, which may require prior linearization of the equations. The numerical solution scheme of the kinetic equa-

K. Park et al. / Ecological Modelling 110 (1998) 281–292

tions alternates between explicit and implicit schemes, and is equivalent to the second-order accurate Crank–Nicholson scheme. Park and Kuo (1996), using a simple hypothetical model, explained the concept of the multi-step computation scheme and demonstrated the enhanced accuracy and computational efficiency comparing it with the single-step explicit or implicit solution of the kinetic terms and with the two-step computation scheme of the decoupled equations. They also discussed the usage of the multi-step scheme for intertidal water quality models. This paper describes the concept and the implementation of the multi-step computation scheme in a more generalized format for both intratidal and intertidal estuarine water quality models. The multi-step scheme is applied to a prototype estuary using a laterally-integrated two-dimensional intratidal model. The model described in Park and Kuo (1993) consists of hydrodynamic and water quality models, and the original version solves the governing massbalance equations for eight water quality state variables by treating the kinetic terms explicitly in a single-step computation. The model results using the multi-step scheme are compared with those of the original version to demonstrate the computational efficiency of the multi-step scheme. 2. Concept of the multi-step computation scheme The concept and the implementation of the multi-step computation scheme, proposed by Park and Kuo (1996), are explained in a generalized format for both intratidal and intertidal estuarine water quality models using a general model consisting of a hydrodynamic model (e.g. solving continuity, momentum, salt-balance and heat-balance equations) and a water quality model (solving mass-balance equations for water quality state variables). The governing mass-balance equation for a water quality state variable may be expressed as: (C ((uC) ((6C) ((wC) + + + (y (z (x (t

=



 

  

283

( ( ( (C (C (C Kx + Ky + Kz + Si (x (y (z (x (y (z + Se,

(2)

where C, concentration of a water quality state variable; u, 6 and w, velocity components in the x-, y- and z-directions, respectively; Kx, Ky and Kz, turbulent diffusivities in the x-, y- and z-directions, respectively; Si and Se, time rates of change resulting from internal and external sources and sinks respectively. The last three terms on the left-hand-side of Eq. (2), the advective terms, and the first three terms on the right-hand-side, the diffusive terms, constitute the physical transport, of which the information is provided by the hydrodynamic model. The last two terms in Eq. (2) represent the kinetic processes and external loads for each of the state variables. The kinetic terms in Eq. (2), while being solved, are decoupled from the physical transport terms. The mass-balance equation for physical transport only, which will be referred to as physical transport equation, is: (C ((uC) ((6C) ((wC) + + + (t (x (y (z =



 

  

( (C ( (C ( (C Kx + Ky + Kz . (x (y (z (x (y (z

(3)

The equation for kinetic processes only, which will be referred to as kinetic equation, is: (C = Si + Se, (t

(4)

which may be expressed as: (C =K · C+R, (t

(5)

where K, kinetic rate (time − 1) and R, source–sink term (mass volume − 1 time − 1). Eq. (5) is obtained by linearizing some terms in the kinetic equations, mostly Monod type expressions. In reality, most biogeochemical processes are non-linear, especially when viewed over a long time period. However over a relatively short time increment, at which numerical models advance their computa-

284

K. Park et al. / Ecological Modelling 110 (1998) 281–292

Fig. 1. The multi-step computation scheme for the governing mass-balance equation, employing an alternate solution of physical transport and kinetic processes. The subscripts K and P indicate kinetic processes and physical transport respectively. The subscript + indicates surplus of processes over Dt, while the subscript – indicates lacking of processes over Dt. Solid-line arrows indicate either kinetic processes or physical transport, while dotted-line arrows displace the same intermediate concentrations to a different notation.

tion, these processes may be approximated by linear representation. Therefore, K and R in Eq. (5) are known quantities in each step of numerical computation, though they may vary with time. Eq. (3) is identical to the salt-balance equation in the hydrodynamic model; therefore, its numerical method of solution is also the same. This paper concerns the method of solution for Eq. (2) in terms of interfacing Eqs. (3) and (4) (or Eq. (5)). The hydrodynamic model may employ either a two- or a three-time level scheme for the physical transport equation (Eq. (3)). We will describe the multi-step computation scheme using a threetime level as an example where the time-step is 2 · Dt = tn + 1 −tn − 1. For models using a two-time level scheme for Eq. (3), the only difference is that the time-step is Dt. In the multi-step scheme, the solution scheme of the kinetic equation (Eq. (5)) is derived by dividing the solution procedure over a time period of 2 · Dt into two steps, alternating between explicit and implicit schemes. Fig. 1(a) illustrates the solution procedure over the time period from tn − 1

to tn + 1. The first step, S1, solves Eq. (5) over Dt from tn − 1 to tn by the explicit scheme: n−1 n C− = Dt · K n − 1 · C n − 1 + Dt · R n − 1, P−C

(6)

which subjects the concentration at t=tn−1, C n − 1, n to the kinetic processes alone to give C − P. The superscript designates the time-step. The subscript − P designates an intermediate concentration that lacks the physical transport over Dt, whereas the subscript + P will designate one with surplus physical transport over Dt. In Fig. 1, the subscript − K designates an intermediate concentration that lacks the kinetic update over Dt, whereas the subscript +K designates one with surplus kinetic n n−1 update over Dt, hence, C − P = C + K . Next, the n−1 intermediate concentration, C + K , is physically transported over 2 · Dt from tn − 1 to tn + 1 in step S2 (Fig. 1(a)) by the finite difference form of Eq. (3): n n n+1 n−1 C− K − C + K = 2 · Dt · PT(C , u ),

(7)

where PT is a physical transport operator over 2 · Dt from tn − 1 to tn + 1, which transports the

K. Park et al. / Ecological Modelling 110 (1998) 281–292

concentration field C n with the velocity field u n, n+1 and C + K is another intermediate concentration at t = tn + 1 lacking the kinetic update over Dt from tn to tn + 1. Eq. (7) is identical with the procedure of the solution of the salt-balance equation. Since the concentration field C n is not available for substance subjected to kinetic processes, we may approximate Eq. (7) by:





n−1 n+1 C+ K +C − K n+1 n−1 C− , un , K −C + K =2 · Dt · PT 2

(8) by assuming: Cn=

n−1 n+1 C+ K +C − K . 2

(9)

Note that the approximation of the concentration field C n (Eq. (9)) results in an implicit solution of n+1 Eq. (8) for C − K . The solution scheme of Eq. (8) may be further simplified to become an explicit or semi-implicit one by further approximation of concentration field. Two-time level solution of the physical transport equation would not require the approximation in Eq. (9). Finally, the step S3 solves Eq. (5) over Dt from tn to tn + 1 by the implicit scheme: n−1 n C n + 1 −C + · C n + 1 +Dt · R n + 1, P =Dt · K

(10) n+1

is the concentration at t =tn + 1. In where C the linearized kinetic equations (Eq. (5)), the kinetic rates, which may depend on other state variables, are evaluated using old conditions, i.e. n n+1 K n − 1. Also note in Eq. (10), C + P =C − K (Fig. 1(a)). The explicit (Eq. (6)) and implicit (Eq. (10)) solution steps of the kinetic equation (Eq. (5)) may be replaced by analytical solution of Eq. (5). In principle, the same three-step procedure, S1, S2 and S3, may be repeated for the next time period from tn + 1 to tn + 3, with the equation for the step S4 given by: n+1 n+2 C− =Dt · K n + 1 · C n + 1 +Dt · R n + 1, P −C

(11) n+2 −P

where C is an intermediate concentration at t =tn + 2 lacking the physical transport over Dt from tn + 1 to tn + 2. In practice, the computation steps S3 and S4 may be combined by adding Eq. (10) and Eq. (11) to arrive at (Fig. 1(b)):

285

n−1 n+2 n C− + K n + 1)C n + 1 P − C + P = Dt(K

+ 2 · Dt · R n + 1,

(12)

which may be approximated by: n n n+2 n n n+2 C− P − C + P = Dt · K (C + P + C − P )+ 2 · Dt · R , (13)

by assuming: K n − 1 + K n + 1 # 2 · K n,

(14)

1 n n+2 C n + 1 # (C + P + C−P ) 2

(15)

and R n + 1 # R n.

(16)

Eq. (13) is a second-order accurate trapezoidal Crank–Nicholson solution of Eq. (5) over 2 · Dt from tn to tn + 2, with the concentration at t=tn + 1 given by Eq. (15). In eutrophication models, the kinetic rate term K has the time scale on the order of hours or longer, and the source–sink term R usually consists of external loads and sediment– water exchange fluxes. In model application, external loads are usually specified as a daily input and sediment–water exchange fluxes have the time scales of days and months. Since Dt in an intratidal model is on the order of minutes, the assumption in Eq. (16) does not significantly affect the accuracy of the solution. With the combination of computational steps S3 and S4 in Fig. 1(a), the solution scheme for Eq. (3) and Eq. (5) becomes an alternate solution of the physical transport (Eq. (7)) and the kinetic processes (Eq. (13)), as illustrated in Fig. 1(b). It should be noted that the intermediate concentrations with the subscripts, either 9 K or 9 P, are not real concentran+1 tions at their time-steps. For example, C − K and n+1 C + K are imaginary concentrations at t=tn + 1, with the former lacking the kinetic update from tn to tn + 1 and the latter with surplus kinetic update from tn + 1 to tn + 2. The real concentrations at t=tn + 1, C n + 1, may be evaluated by the average of these two. The multi-step scheme described above may be generalized for intratidal and intertidal models (Fig. 2). In intratidal (real-time) models, the physical transport processes have shorter time scales

286

K. Park et al. / Ecological Modelling 110 (1998) 281–292

Fig. 2. Generalized usage of the multi-step computation scheme for intratidal and intertidal models. See Fig. 1 for the meaning of subscripts and arrows.

than the water quality kinetic processes. Hence, the kinetic equation (Eq. (5)) may not be solved as often as the physical transport equation (Eq. (3)). In general then, u =m · (2 · Dt) where m is a positive integer (Fig. 2(a)). The kinetic equation is solved once over a time interval of u from tn to tn + 2m for every m steps of computation of physical transport, with the new concentration at t = n n + 2m tn + m given by the average of C + mP and C − mP (Eq. (15)). However, it should be cautioned that u should not be too large to cause instability in the solution of kinetic equation by consuming all materials within a cell over a time period of u. The solution scheme in Fig. 2(a) is implemented in a laterally-integrated two-dimensional intratidal model described in Park and Kuo (1993), and is compared in the next section with the traditional explicit one-step solution of the kinetic terms. In intertidal (tidal-average) models, the physical transport may be updated over a time interval of one tidal cycle, which is too long a time-step for the kinetic processes. In Fig. 2(b), the kinetic equation (Eq. (5)) is solved every Dt for m times over the first half tidal cycle, and then the physical

transport equation (Eq. (3)) is solved once over a time interval of one tidal cycle (u= tn + 2m − tn ), and finally the kinetic equation is again solved every Dt for m times over the second half tidal cycle to give the new concentration C n + 2m.

3. Application of the multi-step computation scheme to an intratidal model The laterally-integrated two-dimensional intratidal model described in Park and Kuo (1993) consists of hydrodynamic and water quality models. The hydrodynamic model solves the laterally integrated continuity, momentum and salt-balance equations using a two-time level finite difference method. The water quality model, using the physical transport information from the hydrodynamic model, solves the laterally integrated massbalance equations for eight water quality state variables (CHL, chlorophyll-a; ON, organic nitrogen; NH4, ammonium nitrogen; NO3, nitrite+nitrate nitrogen; OP, organic phosphorus; PO4, inorganic phosphorus; CBOD, carbonaceous biochemical oxygen demand; DO, dissolved oxygen):

K. Park et al. / Ecological Modelling 110 (1998) 281–292

 

( ( ( (CB)+ (uBC)+ (wBC) (t (x (z =



 

u [C]N [K1]O k = [I] − k 2



(C ( (C ( KxB + KzB +B(Si +Se), (x (z (z (x (17)

with C, laterally averaged concentration of a dissolved or suspended water quality state variable; B, estuarine width. The water quality model in Park and Kuo (1993) solved Eq. (17) in one-step and treated the kinetic terms explicitly. The multi-step computation scheme in Fig. 2(a) is now implemented in this water quality model. The kinetic formulations used for the eight state variables are presented in Park and Kuo (1993). The kinetic equations may be expressed in the format of Eq. (5) by linearizing some terms in the kinetic equations, mostly Monod type expressions. The linearized kinetic equations for the eight state variables can be expressed using an 8×8 matrix: ( [C]=[K] · [C]+ [R], (t

(18)

where [C] is in mass volume − 1, [K] is in time − 1 and [R] is in mass volume − 1 time − 1. Since the settling of particulate matters from the overlying cell acts as an input for a given cell, when Eq. (18) is applied to a cell of finite volume, it may be expressed as: ( [C]k = [K1]k · [C]k +l1[K2]k · [C]k − 1 +[R]k, (t (19) where the four matrices [C]k, [K1]k, [K2]k and [R]k are defined in Appendix A. The subscript k designates a cell at the kth vertical layer. The layer index k increases downward, k =1 is the surface layer and k =KC is the bottom layer. Then, l1 = 0 for k= 1 (surface layer), otherwise l1 = 1. The matrix [K2] is a diagonal matrix, and the non-zero elements (Appendix A) account for the settling of particulate matter from the overlying cell. Eq. (19) is solved using a second-order accurate trapezoidal scheme over a time-step of u, which may be expressed as:



287

−1

u O · [C]O {[K1]O k + k · [C]k 2



A O + l1[K2]O k · [C]k − 1}+ u[R]k ,

(20)

where u= m · Dt is the time-step for the kinetic equations (Fig. 2(a)) and Dt, time-step for a twotime level solution of the physical transport equation, i.e. salt-balance equation in hydrodynamic model; [I] is a unit matrix; [C]A = [C]N + [C]O. The superscripts O and N designate the variables before and after being adjusted for the relevant kinetic processes over time-step. Since Eq. (20) is solved from the surface layer downward, the term with [C]A k − 1 is known for the kth layer and thus placed on the right-hand-side. In Eq. (20), inversion of a matrix can be avoided if the eight state variables are solved in a proper order. The kinetic equations are solved in the order of the variables in the matrix [C] defined in Appendix A. The final forms of Eq. (20) are also listed for each of the state variables in Appendix A. The model described in Park and Kuo (1993) was applied to the tidal Rappahannock River, a western shore tributary of Chesapeake Bay. The results of model calibration and verification are presented in Park and Kuo (1994) for the hydrodynamic model and in Park et al. (1996) for the water quality model. The model calibration results from the original version, which solves the kinetic terms explicitly, are compared in Fig. 3 with those using the multi-step scheme (Fig. 2(a)). Comparisons are made in Fig. 3 for CHL, DO, CBOD, total nitrogen (total organic nitrogen plus total inorganic nitrogen), total organic nitrogen (ON plus the portion in algal biomass), NH4, NO3, total phosphorus (total organic phosphorus plus PO4), total organic phosphorus (OP plus the portion in algal biomass) and PO4. The original version uses a Dt of 108 s for both the hydrodynamic and the water quality models. The modified model with the multi-step scheme uses the same Dt of 108 s for the physical transport equation and u of 108 min (i.e. m= 60 when u= m · Dt) for the kinetic equations. The computation time when solving Eq. (17) for eight state variables using multi-step scheme is

288

K. Park et al. / Ecological Modelling 110 (1998) 281–292

Fig. 3. Comparison of the model results as daily averages between the explicit treatment of the kinetic terms (solid and dashed lines at the surface and bottom layers respectively) and the multi-step scheme ( × and + at the surface and bottom layers respectively: for clarity, only every other points are marked).  and  are field data at the surface and bottom layers respectively.

: 53% of that of original version of the model. The multi-step scheme therefore is computationally much more efficient, since it does not solves the kinetic equations as often as the original

version. The multi-step scheme, using a larger time-step for the kinetic equations, may have larger truncation error. Fig. 3 however shows that virtually no difference exists between the original

K. Park et al. / Ecological Modelling 110 (1998) 281–292

289

Fig. 3. (Continued)

version and the multi-step scheme, indicating that the usage of a second-order accurate trapezoidal scheme over the time interval of 108 min (Eq. (20)) gives as accurate results as the explicit scheme with the time interval of 108 s.

4. Discussion The time-step in intratidal models is governed more by the physical transport processes than by the kinetic processes. To relieve the constraint in time-step and to achieve high accuracy, secondorder accurate finite difference schemes, either three-time level or two-time level (Smolarkiewicz and Margolin, 1993), have been used in some hydrodynamic models (Wang, 1982; Blumberg and Mellor, 1987; Moustafa and Hamrick, 1994). When coupling a biogeochemical model to these

hydrodynamic models, the kinetic processes have been routinely treated explicitly, which is first-order accurate, so that additional complications are not introduced into the solution scheme, which results in inconsistent treatment of physical transport and kinetic processes. The multi-step scheme for intratidal models (Fig. 2(a)) is second-order accurate in treating the kinetic processes, which is consistent with second-order accuracy in many hydrodynamic models. The scheme in Fig. 2(a) employs a larger time-step for the kinetic terms. It has thus enhanced computational efficiency, compared with the traditional explicit or implicit solution of the mass-balance equations, which solves the kinetic processes unnecessarily too often with the same time-step as the physical transport processes. The truncation error in the numerical solutions increases with a larger time-step, which is compensated for by the increased accuracy of the second-order solution, as shown in Fig. 3.

290

K. Park et al. / Ecological Modelling 110 (1998) 281–292

The time-step in intertidal models is governed more by the kinetic processes than by the physical transport processes. In the application of the Chesapeake Bay three-dimensional water quality model, e.g. physical transport and kinetic processes were computed on a 2-h time-step, even though the hydrodynamic model provided information of intertidal (12-h) physical transport (Cerco and Cole, 1993). The multi-step scheme for intertidal models (Fig. 2(b)) requires the computation of the physical transport processes only once over a tidal cycle with multiple updates of kinetic processes, which again saves the computation time. The multi-step scheme for intratidal models (Fig. 2(a)) was implemented in a three-dimensional water quality model (Park et al., 1995), which is being directly coupled to the three-dimensional hydrodynamic model described in Moustafa and Hamrick (1994). The scheme for intertidal models (Fig. 2(b)) was implemented in a tidal prism water quality model (Kuo and Park, 1995), and is being implemented in a three-dimensional water quality model, which will be indirectly coupled to an intratidal hydrodynamic model (KORDI, 1997). The same scheme may also be applied to an ecological model if it is to be coupled with physical transport processes.

Acknowledgements This paper is contribution number 2120 of the School of Marine Science/Virginia Institute of Marine Science, The College of William and Mary, Virginia, USA. The study was supported partly by the Virginia Department of Environmental Quality’s Coastal Resources Management Program through Grant cNA47OZ0287-01 of the National Oceanic and Atmospheric Administration, Office of Ocean and Coastal Resource Management, under the Coastal Zone Management Act of 1972, as amended, and partly by the G-7 Environmental Research Program, Ministry of Environment and Ministry of Science and Technology, Republic of Korea (Phase II).

Appendix A

The four matrices in Eq. (19) are defined as:

Á CHL Â Ã Ã ON Ã Ã NH4 Ã Ã NO3 [C]k = Ã Ã; OP Ã Ã PO4 Ã Ã CBOD Ã Ã Ä DO Åk Á Ãa 1 0 Ãa 2 b 2 Ãa 3 b 3 Ã [K1]k= Ãa 4 0 Ãa 5 0 Ãa 6 0 Ãa 7 0 Ãa 8 0 Ä

0 0 0 0 0 0 c3 0 0 c4 d4 0 0 0 e5 0 0 e6 0 0 0 c8 0 0

Át 1 à à t à 2à 0 à à 0 entii([K2]k )= à à , t5 à à t à 6à t à 7Ã Ä 0 Åk

0 0 0 0 0 f6 0 0

0 0Â Ã 0 0 Ã 0 0 Ã 0 0 Ã; 0 0 Ã 0 0 Ã g7 0 Ã g 8 h 8 Åk

Ár 1 à à r à 2à r à 3à r [R]k = à 4à , r à 5à r à 6à r à 7à Är 8Åk

where the non-zero elements in [K1]k, [K2]k and [R]k are given below. Hereinafter, the subscript k, which designates the kth layer, is omitted. WS O CHL a1 = G O − R O − P O − , Dz O t1 =

(WS O CHL)k − 1 , Dz O

r1 =

WO CHL , VO

a2 = FON(R O + P O)AN, b2 = − K O ON −

WS O ON , Dz O

t2 =

(WS O ON)k − 1 , Dz O

K. Park et al. / Ecological Modelling 110 (1998) 281–292

r2 =

WO BFO ON ON Bk −l2Bk + 1 + , O V Dz O Bk

a3 = [(1− FON)(R O +P O) − PR OG O]AN, b3 = K O ON, r3 =

c3 = −K O NIT,

WO BFO NH4 NH4 Bk −l2Bk + 1 + , O V Dz O Bk

a4 = − (1− PR O)G OAN,

c4 =K O NIT,

d4 = − K O DENIT, r4 =

WO BFO NO3 NO3 Bk −l2Bk + 1 + , O V Dz O Bk

a5 = FOP(R O +P O)AP,

e5 = −K O OP −

t5 =

(WS O OP)k − 1 , Dz O

r5 =

WO BFO OP OP Bk −l2Bk + 1 + , VO Dz O Bk

a6 = [(1− FOP)(R O + P O) −G O]AP, f6 = − r6 =

WS O PO4 , Dz O

t6 =

e6 =K O OP,

(WS O PO4)k − 1 , Dz O

WO BFO PO4 PO4 Bk −l2Bk + 1 + , O V Dz O Bk

a7 = ACOP OAC,

g7 = −K O CBOD −



c8 = − ANOK

r8 = h8DOs +



u CHLN = CHLO + {a1CHLO + l1t1CHLA k − 1} 2



, + ur1





RO AC, RQ

O NIT

,

g8 = − K O CBOD,

with Dz, layer thickness; G, R and P, algal growth, respiration and predation rates respectively; WSCHL, WSON, WSOP, WSPO4 and WScbod settling rates for CHL, ON, OP, PO4 and CBOD respectively; WCHL, WON, WNH4, WNO3, WOP, WPO4, WCBOD and WDO, external loads for eight state variables respectively; FON and FOP, fractions of metabolically produced nitrogen and phosphorus, respectively, returned to the organic pool; AN, AP and AC, ratios of nitrogen, phosphorus and carbon, respectively, to CHL in algae; KON and KOP, mineralization rates of ON and OP respectively; BFON, BFNH4, BFNO3, BFOP and BFPO4, benthic fluxes of ON, NH4, NO3, OP and PO4 respectively; PR, algal preference for NH4 uptake; KNIT, nitrification rate; KDENIT, denitrification rate; KCBOD, CBOD decay rate; SOD, sediment oxygen demand; PQ and RQ, photosynthesis and respiration quotients respectively; Kr, surface reaeration rate; DOs, saturated DO concentration. In the sediment–water exchange terms in the matrix [R]k, l2 = 0 for k= KC (bottom layer), otherwise l2 = 1. The terms for dissolved oxygen reaeration in h8 and r8 are applied only when k=1 (surface layer), hence l3 = 1 for k= 1, otherwise l3 = 0. The kinetic formulations are presented in detail in Park and Kuo (1993). Eq. (20) for each of the eight state variables are:



u 1− a1 2

−1

,

ONN = ONO

WO SODO Bk −l2Bk + 1 CBOD + , O V Dz O Bk

a8 = 2.67 PQ · G O −

O CBOD O

WS Dz

(WS O CBOD)k − 1 , t7 = Dz O r7 =

WS O OP , Dz O

291

h8 = −l3K O r ,

WO SODO Bk −l2Bk + 1 DO − , VO Dz O Bk

u + {a2CHLA + b2ONO + l1t2ONA k − 1} 2



+ur2



u 1− b2 2



−1

,

u NH4 = NH4O + {a3CHLA + b3ONA + c3NH4O} 2



+ ur3



u 1− c3 2

−1

,

292

K. Park et al. / Ecological Modelling 110 (1998) 281–292



NO3N = NO3O u + {a4CHLA +c4NH4A +d4NO3O} 2





+ur4



u 1− d4 2

−1

,

u OPN= OPO + {a5CHLA +e5OPO +l1t5OPA k − 1} 2



+ ur5





u 1− e5 2

−1

,

PO4N = PO4O u + {a6CHLA +e6OPA 2 A +f6PO4Ol1t6PO4k−1 }+ur6







u 1 − f6 2

−1

,

CBODN = CBODO u + {a7CHLA +g7CBODO 2





+l1t7CBODA k − 1} +ur7



u 1 − g7 2

−1

DON = DOO u + {a8CHLA +c8NH4A +g8CBODA 2



+ h8DOO}+ ur8



u 1 − h8 2

−1

where C A = C O +C N and l1 =0 for k =1 (surface layer), otherwise l1 =1 (see Eq. (20)). References Ambrose Jr., R.B., Wool, T.A., Connolly, J.P., Schanz, R.W., 1988. WASP4, a hydrodynamic and water quality model: model theory, user’s manual and programmer’s guide. EPA/600/3-87/039, U.S. Environmental Protection Agency, Athens, GA, 297 pp. Blumberg, A.F., Mellor, G.L., 1987. A description of a threedimensional coastal ocean circulation model. In: Heaps, N.S. (Ed.), Three-dimensional Coastal Ocean Models. American Geophysical Union, Washington, DC, pp. 1–16.

.

Cerco, C.F., Cole, T.M., 1993. Three-dimensional eutrophication model of Chesapeake Bay. J. Environ. Eng. ASCE 119 (6), 1006 – 1025. Korea Ocean Research and Development Institute (KORDI), 1997. Development of coastal water quality assessment and prediction technology. BSPN00303-945-2, Annual Report of Marine Environmental Monitoring and Assessment Technology (Phase II), Ministry of Environment and Ministry of Science and Technology, Korea, 356 pp (in Korean with English abstract). Kuo, A.Y., Neilson, B.J., 1988. A modified tidal prism model for water quality in small coastal embayments. Water Sci. Technol. 20 (6/7), 133 – 142. Kuo, A.Y., Park, K., 1995. A PC-based tidal prism water quality model for small coastal basins. In: Brebbia, C.A., Traversoni, L., Wrobel, L.C. (Eds.), Computer Modelling of Seas and Coastal Regions II. Computational Mechanics Publications, Southhampton, UK, pp. 371 – 378. Moustafa, M.Z., Hamrick, J.M., 1994. Modeling circulation and salinity transport in the Indian River Lagoon. In: Spaulding, M.L., Bedford, K.W., Blumberg, A.F., Cheng, R.T., Swanson, J.C. (Eds.), Estuarine and Coastal Modeling III. ASCE, New York, pp. 381 – 395. Park, K., Kuo, A.Y., 1993. A vertical two-dimensional model of hydrodynamics and water quality in estuaries. Special Report in Applied Marine Science and Ocean Engineering (SRAMSOE) No. 321, School of Marine Science(SMS)/ Virginia Institute of Marine Science (VIMS), The College of William and Mary (CWM), VA, 47 pp. Park, K., Kuo, A.Y., 1994. Numerical modeling of advective and diffusive transport in the Rappahannock Estuary, Virginia. In: Spaulding, M.L., Bedford, K.W., Blumberg, A.F., Cheng, R.T., Swanson, J.C. (Eds.), Estuarine and Coastal Modeling III. ASCE, New York, pp. 461 – 474. Park, K., Kuo, A.Y., 1996. A multi-step computation scheme decoupling kinetic processes from physical transport in water quality models. Water Res. 30 (10), 2255 – 2264. Park, K., Kuo, A.Y., Shen, J., Hamrick, J.M., 1995. A threedimensional hydrodynamic-eutrophication model (HEM3D): description of water quality and sediment process submodels. SRAMSOE No. 327, SMS/VIMS, CWM, VA, 102 pp. Park, K., Kuo, A.Y., Neilson, B.J., 1996. A numerical model study of hypoxia in the tidal Rappahannock River of Chesapeake Bay. Estuar. Coast. Shelf Sci. 42 (5), 563 – 581. Shanahan, P., Harleman, D.R.F., 1982. Linked hydrodynamic and biochemical models of water quality in shallow lakes. Report No. 268, Ralph M. Parsons Laboratory, Massachusetts Institute of Technology, MA, 279 pp. Smolarkiewicz, P.K., Margolin, L.G., 1993. On forward-intime differencing for fluids: extension to a curvilinear framework. Mon. Weather Rev. 121, 1847 – 1859. Valocchi, A.J., Malmstead, M., 1992. Accuracy of operator splitting for advection – diffusion – reaction problems. Water Resour. Res. 28 (5), 1471 – 1476. Wang, D.-P., 1982. Development of a three-dimensional, limited-area (island) shelf circulation model. J. Phys. Oceanogr. 12, 605 – 617.