Mathematical model of the oxygen regime of Cartagena Bay

Ecological Modelling 165 (2003) 91–106

Mathematical model of the oxygen regime of Cartagena Bay Yuri S. Tuchkovenko a,∗ , Serguei A. Lonin b,1 a

b

Institute of Biology of Southern Seas, National Academy of Sciences of Ukraine, Odessa Branch, 37 Pushkinskaya St., 65011 Odessa, Ukraine Centro de Investigaciones Oceanográficas e Hidrográficas, Escuela Naval, A.A. 982, Cartagena de Indias, Colombia Received 14 November 2001; received in revised form 26 November 2002; accepted 7 January 2003

Abstract A mathematical model of eutrophication/oxygen regime for Cartagena Bay (Colombia) is described. The model consists of hydrodynamic and biogeochemical blocks to assess engineering alternatives for improving the water quality and oxygen regime in the bay. The hydrodynamic block of the model is based on the primitive three-dimensional equations for estuarine and coastal dynamics. The biogeochemical block includes the balance equations for the following components of the ecosystem: phytoplankton, bacteria, detritus, dissolved organic matter (DOM), phosphate, ammonium, nitrite, nitrate and dissolved oxygen. A non-traditional parameterization of bacteria role in biochemical oxidation of organic matter is used. Modelling of distinct ecological situations for Cartagena Bay has shown the principal strategic objectives to improve the water quality of the bay. The model has been tested and could be applied to other similar water bodies. © 2003 Elsevier Science B.V. All rights reserved. Keywords: Eutrophication; Dissolved oxygen; Water quality; Three-dimensional model; Hydrodynamics; Bay model

1. Introduction The natural biotic cycle of matter in marine ecosystems is modelled so that a nutrient balance in the “water-biota” phase exists. Primary production by phytoplankton associated with the consumption of mineral forms of nitrogen and phosphorus is regulated by nutrient regeneration due to the biochemical oxidation of metabolic excretions and dead organic matter in the presence of bacteria. Uncontrolled anthropogenic influx of nutrients increases primary production and a part of organic matter, which does not have sufficient time to be accumulated in higher trophic layers, and when no longer ∗ Corresponding author. Tel.: +38-0482-473992. E-mail addresses: [email protected] (Y.S. Tuchkovenko), [email protected] (S.A. Lonin). 1 Tel.: +57-56-694465; fax: +57-56-694297.

living, produces favorable conditions for bacteria growth, with the result that the nutrient cycle period is shortened, and the oxygen demand for biochemical oxidation of organic matter is increased. The consequences of the anthropogenic eutrophication process are: (a) chemical regime changes; (b) oxygen deficiency or its total disappearance in the bottom layer (anoxic conditions); (c) strong decline of habitat for aerobic aquatic organisms; (d) destruction of the equilibrium between the production and decomposition processes, stability, trophic structure and dynamics of the ecosystem. A key water quality indicator for eutrophic water bodies is seasonal oxygen dynamics (percentage of saturation). This characteristic is an integral indicator for the production–destruction relationship of organic matter and for the influence of hydrometeorological processes. Oxygen decreasing below some threshold value (about 1 ml l−1 ) causes mortality of various

0304-3800/03/$ – see front matter © 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0304-3800(03)00064-4

92

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

aerobic aquatic organisms and alters many chemical– biological processes. In this case, the recreation potential of the water is entirely lost and its quality becomes unacceptable according to sanitary-hygienic norms. Because the oxygen content depends on numerous chemical–biological processes with different intensities, the planning and forecasting of environmental protection and control measures and their efficacy in reducing eutrophication is impossible without modelling. In this work, such a model is considered. The aim of the present paper is to describe this eutrophication

model, which has been applied to the eutrophication problem in the bay in order to evaluate some recommendations for improving of its water quality. Before the applications the model has been calibrated with available data from ecological monitoring.

2. Study area Cartagena Bay is located on the Colombian coast of the Caribbean Sea (10◦ 16 –10◦ 26 N and 75◦ 30 –75◦ 35 W; Fig. 1). The bay has the following

Fig. 1. Location of the principal industrial and domestic sources of Cartagena Bay.

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

morphologic characteristics: maximum meridian length is about 16 km, latitudinal length is 9 km; water surface area is about 82.6 km2 , mean depth is 16 m and the maximum depth is 26 m. The bay consists of two parts: an external bay, which is connected to the Caribbean Sea through two straits, and an internal bay, situated on the north and without connection with the sea, except through the external bay. On the coast of the internal bay is situated the City of Cartagena de Indias, 40% of the sewerage water of the city directly enters the bay. The eastern coast of the external bay is an industrial zone of Cartagena, which includes 29 industries. Having a total water discharge of about 1.42 × 106 m3 per day, those anthropogenic sources transport into the bay 2.57 tons of mineral forms of nitrogen, 0.48 tons of mineral phosphorus and 22.2 tons (in units of biological demand of oxygen, BOD) of dead organic matter daily. Fresh waters enter the external bay from the Dique channel, which is connected with the Magdalena River. This channel is artificial, and its flow rate varies during the year from 55 m3 s−1 for the dry season (February–April) to 250 m3 s−1 for the rainy one (September–October). The channel influence on the hydrodynamic regime of the bay is two-fold. On the one hand, the channel provides mineral forms of ni-

93

trogen and phosphorus, and also suspended minerals that determine the water transparency of the bay (it varies from 0 to 1.5 m). On the other hand, because of the lower density of fresh water, a strong subsurface pycnocline (upper 4 m) is formed, which decreases vertical mixing. As a result, pollution, transported by the channel, is distributed into the surface freshwater layer, which lies entirely in the photic one. Depth of the straits between the bay and the sea varies from 0.5 to 3 m, excluding the navigation channel, situated in the southern strait, which is 100 m wide and 30 m deep. Because the bay depth reaches 26 m, hydrodynamic renewal of the bottom layer by the marine water is difficult. The horizontal advection of relatively pure marine waters, favourable for water renewal, is caused by tidal action at the marine boundary, but the tidal currents are weak due to the small amplitude of tides (20–40 cm). The marine water has a greater density (due to a high salinity) than the transformed fresh water in the bay; after entering through the shallow straits, the salt water is, therefore, passed to the deeper part of the bay, refreshing the bottom layer. The wind regime is characterised by strong (daily average velocity of 8 m s−1 ) north-easterly trade winds for the dry season and weak (up to 3 m s−1 ) winds for the wet season.

Fig. 2. Measured oxygen content (in ml l−1 ) in the bottom layer during dry (a) and wet (b) seasons.

94

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

During the wet season, when the channel flow rate is maximum and the vertical exchange is minimum due to the absence of wind, a strong pycnocline is formed and nutrients from the Dique channel and from other anthropogenic sources are distributed within the upper fresh water layer. The phytoplankton growth is limited in the major part of the bay by light because of turbidity. Nevertheless, the transparency increases in the vicinity of the straits due to deposition of the mineral particles and the mixing with the relatively pure sea water. Herein, local blooms of phytoplankton are observed. The autoctone organic matter, formed by phytoplankton photosynthesis, sinks into deeper layers, where it is subject to bacterial decomposition. A part of dissolved oxygen is consumed by organic matter oxidation and nitrification. The morphological features of the bay (shallow straits and deeper bay), a weak vertical turbulent exchange and absence of photosynthetic production beyond the pycnocline are the causes of the oxygen deficit for the bottom layer (Fig. 2). The aerobic hydrochemical conditions become anaerobic in a major part of the bottom layer, and this circumstance causes mortality of aerobic organisms. Biochemical decomposition of dead organic matter is decelerated and, as a consequence, the undecomposed organic matter is deposited into the bottom layer. After the rainy season (beginning in January), the Dique channel flow rate is minimal, the subsurface pycnocline is decreased in depth, and the water exchange and the vertical turbulent mixing are increased by wind action. As a result, the oxygen content in the bottom layer is increased by up to 2–3 ml l−1 (see Fig. 2a). Thus, a strong anthropogenic impact on the Cartagena Bay ecosystem from the industrial and domestic sources, pollution from Magdalena River together with the morphometric features of the basin cause the eutrophication in the bay and periodic oxygen deficit in its bottom layer.

3. Mathematical model 3.1. Hydrodynamic block The theoretical basis of the hydrodynamic model is the three-dimensional Model of Estuarine and Coastal Circulation Assessment (MECCA, Hess, 1985), modi-

fied by Lonin (1997) to the CODEGO model (Codego is the ancient indigenous name of Cartagena Bay). The basic equations are as follows: ∂uv ∂uw ∂u ∂u2 + + + ∂t ∂x ∂y ∂z
 ∂P ∂u ∂ + fv + 2Ah = −α0 ∂x ∂x ∂x
 
 ∂u ∂v ∂ ∂u ∂ + Ah + + Av , ∂y ∂x ∂y ∂z ∂z ∂v ∂uv ∂v 2 ∂vw + + + ∂t ∂x ∂y ∂z
 ∂ ∂P ∂v − fu + = −α0 2Ah ∂y ∂y ∂y
 
 ∂v ∂v ∂u ∂ ∂ Ah + + Av , + ∂x ∂x ∂y ∂z ∂z

(1)

(2)

∂P = −ρg, ∂z

(3)

∂u ∂v ∂w + + = 0, ∂x ∂y ∂z

(4)

ρ = ρ0 [1 + Fρ (S, T )],

(5)

∂uS ∂vS ∂wS ∂S + + + ∂t ∂x ∂y ∂z


 ∂ ∂S ∂S ∂S ∂ ∂ = Dh + Dh + Dv , ∂x ∂x ∂y ∂y ∂z ∂z (6) ∂T ∂uT ∂vT ∂wT + + + ∂t ∂x ∂y ∂z

 ∂ ∂T ∂T ∂ = Dh + Dh ∂x ∂x ∂y ∂y
 ∂ ∂T + Dv + R, ∂z ∂z

(7)

where u, v and w, velocity components for x, y and z (vertical), respectively; t, time; f, Coriolis parameter; P, pressure; g, gravity; ρ, water density; α 0 , specific volume; Ah and Av , eddy viscosity in horizontal and vertical directions, respectively; and Dh , Dv , turbulent diffusion on these directions; T, temperature; S, salinity; R, internal heat sources. The boundary conditions at the sea surface (z = 0) for (1)–(7) are defined as follows: dynamic conditions

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

for the wind stresses (τ sx , τ sy ) and the atmospheric pressure Pa , heat QT and salinity QS fluxes and the kinematic condition for the free surface h, i.e.
 ∂u ∂v (τsx , τsy ) = ρAv , , P = Pa , ∂z ∂z

 ∂T ∂S QT dh Dv , = , (8) w= , QS , ∂z ∂z ρCw dt where Cw , specific heat capacity of water. For the bottom (z = H ), the drag law parameterization and zero-flux conditions for (6) and (7) are:
 ∂u ∂v (τbx , τby ) = ρAv , , ∂z ∂z
 ∂T ∂S , Dv = (0, 0). (9) ∂z ∂z At the open boundary of the bay, tidal oscillations h0 , patterns of temperature T∗ and salinity S∗ and the numerical outflow conditions are specified as:   (T , S) = (T∗ , S∗ ), if vn ≤ 0 , ∂(T , S)  ∂(T , S) + | v| = 0, if vn > 0 ∂t ∂n  H h = h0 + vn , (10) g where n, unit normal to the boundary with the calculated vector v of currents; H = h + d, total depth; d, local depth. The second condition in (10) expresses both the specified tidal input h0 and long-wave reflex ion in the linear approach. The flow rate Qr , temperature fT (z) and salinity fS (z) profiles are specified at the Dique channel entrance and a no-flow condition is used at the solid lines. The vertical turbulence is described by Prandtl’s concept of mixing length Λ = [κz(1 − z/H )]Θ(z), where κ is the Karman’s constant; Θ is a limiting profile function. Following Monin and Yaglom (1971) the stationary vertical turbulence in the presence of buoyancy would be expressed as

2 1/2 ∂u 2 ∂v 2 Az = Λ + (1 − CR1 Rf )1/2 , ∂z ∂z where CR 1 , the critical inverse Richardson number (1/Ricr ), i.e. turbulence exists when Ricr < 1/4

95

(CR1 = 4.0); Rf, the dynamic Richardson number defined as


g Rf = − ρ0




αρ

∂ρ ∂z




∂u ∂z




2 +

∂v ∂z

2 −1 ,

αρ ≈ 1.0. The transport equation for a substance (i) of concentration Ci representing one of the chemical–biological compounds of the ecosystem has the following form (similar to (6) and (7), but with a non-zero settling velocity wg in the case of solid particles such as phytoplankton or detritus): ∂HC ∂uHC ∂vHC ∂(w + wg )C + + + ∂t ∂x ∂y ∂σ

  ∂ ∂C ∂ ∂C Dh H Dh H − − ∂x ∂y ∂x ∂y
 ∂C ∂ Dv H −1 = F. − ∂σ ∂σ

(11)

Herein, the coordinate z has just been transformed to the non-dimensional coordinate σ = (z − h)/H. The right hand part of (11) represents the chemical–biological internal sources, described later (see Eqs. (12)–(19)). The numerical schemes used to solve Eqs. (6) and (7) were modified to be used in (11), because the original schemes (Hess, 1985) had conservative properties but had no transport (monotonous). Because of the horizontal discrete operator of (11) explicitly presented in time, the Flux Corrected Transport Scheme (Boris and Book, 1976) was applied, previously transforming (11) into the Lax–Vendroff’s finite-difference scheme. In contrast, the vertical discrete operator of (11) is in an implicit form, which may be presented as the Crank–Nicholson’s scheme and later transformed into the Total Variation Diminishing Scheme (Fletcher, 1991). Parasite negative oscillations in the solution tail are avoided when using those transport schemes, but a little dissipation effect is introduced so that the solution will not be completely conservative. Numerical details of the employed schemes are presented in Lonin (1997). The following block contains the detailed expressions for F in (11) and represents the local internal chemical–biological processes and the external fluxes.

96

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

Fig. 3. Diagram of chemical–biological structure of the eutrophication model and relationships among its components.

These fluxes are implicitly related to the boundary conditions for the last equation. 3.2. Biogeochemical block It includes the following ecosystem variables: Bf , DET , detritus, phytoplankton, Bb , bacterioplankton, Borg DIS = B ant + B nat , dissolved organic matter (DOM), Borg org org CPO4 , phosphate, CNH4 , ammonium, CNO2 , nitrite, CNO3 , nitrate, CO2 , dissolved oxygen. The diagram of relationships among the ecosystem components of the biogeochemical block is presented in Fig. 3. The transformation equations in the local interpretation are presented below.

tosynthesis (W m−2 ); Iz , luminosity at depth z (W m−2 ); α = α0 + αsus + αf , integral light extinction with depth (m−1 ), divided into components of the extinction of oceanic waters α 0 , extinction due to mineral suspension α sus and phytoplankton α f 0.542 (CIOH, 1997), (self-shading): αsus = 1.31Csusp 2 αf = 0.0088Bfchl + 0.054Bfchl (Parsons et al., 1984), where Csusp , mineral suspension concentration (mg l−1 ) and Bfchl , phytoplankton biomass, (mg ChlA m−3 ); CPO4 ,CNH4 , CNO3 , concentrations of phosphate (mg P l−1 ), ammonium (mg N l−1 ) and nitrate (mg N l−1 ), respectively; CkN , CkPO4 , half-saturation constants for the rate of utilization of the mineral forms of nitrogen and phosphorus by phytoplankton (mg l−1 ); γ f , part of phytoplankton

Phytoplankton:  dBf  1 z2 max = (1 − γ )σ (I , C , C )B − µ B , σ = V f (I )f (C , C ), f (I ) = fz (Iz )dz PO4 N f 1 2 N PO4 1 f f z f f f f dt local ,z z1 2.718 I0 = [exp(−Rz1 ) − exp(−Rz2 )], ,z = z2 − z1 , Rzi = exp(−αzi ), ,zα Iopt
 

Iz Iz CN CPO4 exp 1 − fz (Iz ) = , Iz = I0 exp(−αz), f2 (CN , CPO4 ) = min , , Iopt Iopt CkN + CN CkPO4 + CPO4 CN = CNH4 + CNO3 . Here, t, time (h); Bf , phytoplankton biomass (mg C m−3 ); σ f , specific rate of phytoplankton growth, determined by the illumination conditions (I) and the presence of nutrients (mineral forms of nitrogen CN and phosphorus CPO4 ); Vfmax , maximal specific rate of phytoplankton growth (h−1 ); I0 , photosynthetic active radiation flux at the sea surface (W m−2 ); Iopt , optimal luminosity for pho-

(12) production, demanded in respiration; µf , specific mortality rate of phytoplankton (h−1 ). Bacterioplankton:    DIS Borg dBb  Bb max =V ε − max Bb , DIS ing dt local b Bb Bkorg + Borg (13)

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106 DIS , where Bb , bacteria biomass (mg C m−3 ); Borg concentration of DOM in the water (mg O2 l−1 ), ant and natuwhich is the sum of anthropogenic Borg max nat ral Borg parts; Vb , maximal rate of growth for bacterioplankton (h−1 ); Bkorg , half-saturation constant for the growth of bacteria that is equal to concentration of the organic substrate when the real specific growth is equal to half of maximal one (mg O2 l−1 ); Bbmax , maximal possible biomass of bacteria (mg C m−3 ). The structure of (13) takes into account the influence of bacteria population density on the bacteria mortality. For magnitudes of Bb smaller than Bbmax , Eq. (13) is transformed into the Monod equation; DIS larger than Bk , it is for magnitudes of Borg org turned into an equation of logistic growth (INBUM, 1991). These asymptotes provide a necessary flexibility to (13). Dead organic matter:  DET  dBorg  DET = λf ηf µf Bf βO2 /C βm3 /l − δBorg ,  dt  local

 nat  dBorg   dt 

local

 ant  dBorg   dt 

local



(14)

= (γf σf (Iz , CPO4 , CN ) + (1 − λf )ηf µf )Bf + µb BB 

 nat Borg VBmax BB εing − DIS θ Bkorg + Borg DET × βO2 /C βm3 /l + δBorg , (14a)   ant Borg VBmax = Qant org − DIS θ Bkorg + Borg

× BB εing βO2 /C βm3 /l .

(14b)

DET , concentration of particulate dead Here, Borg organic matter (detritus) (mg O2 l−1 ); δ, specific rate of detritus autolysis (h−1 ); λf , fraction of detritus in phytoplankton cells (0 < λf < 1); ηf , labile fraction in dead organic matter of phytoplankton; µB = VBmax BB2 /BBmax , specific rate of bacterioplankton mortality (h−1 ) (INBUM, 1991); θ, non-dimensional constant expressing the effectiveness of substrate utilization by bacteria for its growth, which is determined as θ = Pb /(Pb + Rb )

97

where (Pb + Rb ) is energy accumulated in bacteria, which is summed from their production Pb and outgoing exchange Rb (IOAN, 1989; INBUM, 1991); εing = CO2 /(CO2 + CkO2 ), inhibition parameter for the organic matter biochemical oxidation and nitrification under conditions of oxygen deficit (0 < εing < 1), where CkO2 , half-saturation constant of the process (mg O2 l−1 ); Qant org , content variations for dead organic matter, due to anthropogenic sources; βm3 /l = 0.001, conversion coefficient for concentration from cubic meters to litres; βO2 /C , transition coefficient from carbon units of organic matter (mg C) to oxygen units (mg O2 ) (mg O2 mg−1 C). Phosphate:  dCPO4  dt local   

 1 1 nat ant ant = + Borg − 1 βP/C Borg − ωP βP/C θ θ max VB × B β 3 ε DIS b m /l ing Bkorg + Borg − σf (Iz , CPO4 , CN )Bf βP/C βm3 /l + Qant PO4 .

(15)

Ammoium:  dCNH4  dt local 

  1 1 nat ant ant = − 1 Borg βN/C + − ωN Borg βN/C θ θ VBmax × B β 3 ε DIS b m /l ing Bkorg + Borg − χ σf (Iz , CPO4 , CN )Bf βN/C βm3 /l − νN1 CNH4 εing + Qant NH4 .

(16)

ant , ω = β ant Here, ωP = βP/C /βP/C N N/C /βN/C , relationships between phosphorus and nitrogen contents in the natural and anthropogenic organic matter; ant , β nat , C:N:P ratios of natural βP/C , βN/C , βP/C N/C nat and anthropogenic B ant dead organic matter Borg org in (mg P mg−1 C) and (mg N mg−1 C), correspondingly; χ = CNH4 φ/[φCNH4 +(1−φ)CNO3 ], fraction of mineral nitrogen, consumed by phytoplankton in the ammonium form, where φ, preference coefficient of ammonium assimilation by phytoplankton ant ant ant with respect to nitrate; Qant PO4 , QNH4 , QNO2 , QNO3 ,

98

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

content variations for mineral forms of nitrogen and phosphorus, due to anthropogenic sources. Nitrite:  dCNO2  = νN1 CNH4 εing − νN2 CNO2 εing + Qant NO2 , dt local (17) where CNO2 , concentration of nitrite (mg N l−1 ); νN1 , νN2 , specific rates of the first and second stages of nitrification (h−1 ). Nitrate:  dCNO3  = νN2 CNO2 εing − (1 − χ ) dt local × σf (Iz , CPO4 , CN )Bf βN/C βm3 /l crit − νDN (CNO3 − CNO ) 3 crit − νfoto (CNO3 − CNO ) + Qant NO3 , 3

(18) where νDN , specific rate of denitrification (h−1 ) in the bottom layer, where oxygen concentration is less than 1 mg l−1 ; νfoto , specific rate of nitratereduction (h−1 ) due to physical–chemical processes in the surface layer (Raymond, 1983); crit , minimal nitrate concentration, when the CNO 3 nitrate-reduction processes are stopped (mg l−1 ). Dissolved oxygen:  dCO2  = σf (I, CPO4 , CN )Bf dt local 

 DIS Borg Vbmax Bb εing − DIS θ Bkorg + Borg × β O2 /C βm3 /l − (βO2 /N1 νN1 CNH4 + νN2 βO2 /N2 CNO2 )εing atm − Qbot O2 ± QO2 .

(19)

Here, CO2 , concentrations of oxygen (mg O2 l−1 ); βO2 /N1 , βO2 /N2 , oxygen equivalents for first and second stages of nitrification (mg O2 mg−1 N). Bottom absorption of oxygen may be estimated in first approximation according to Parsons et al. (1984) through the oxygen concentration in the water CO2 : b Qbot O2 = a f (x, y)(CO2 ) ,

(20)

−2 h−1 , C = ml l−1 ; a, where Qbot O2 O2  = mg m b, empirical constants; f(x, y), function that describes the spatial variations of the oxygen bottom act (x, y)/F act , flux, determined as f (x, y) = Forg med act (x, y), the organic matter flux through where Forg the sediments, spatially distributed, obtained in the act , the mean organic flux for the actual model;Fmed conditions. The bottom absorption flux of oxygen is due to organic matter sedimentation and oxidation. In this model, only two variables have non-zero gravity velocities: phytoplankton and detritus. The function f(x, y) in (20) takes into account the spatial non-homogeneity of this flux and may be corrected in the case when productivity changes by altering the anthropogenic sources. The flux of oxygen through the surface is described according to Lyakhin (1980) in the following equation: S Qatm O2 = γe,i nv nt (CO2 − CO2 ),

(21)

−2 where Qatm O2 , oxygen invasion or evasion, (mg m h−1 ); γ i,e , coefficient of invasion (evasion) (l m−2 h−1 ); nt , the temperature coefficient (for T = 30 ◦ C, S = 30–35, nt = 1.1), nv is the integral wind factor:  1.0 + 0.27W 2 , for W ≤ 8 m s−1 nv = , (22) −7.4 + 0.4W 2 , for W > 8 m s−1

where W, the wind speed (m s−1 ) at 10 m; COS 2 , saturated oxygen concentration (mg l−1 ) (at given temperature T and salinity S). Mathematical structure of the biogeochemical block is based on a traditional description (Ambrose et al., 1993; Sarmiento et al., 1993; Cerco and Cole, 1995; Lonin and Tuchkovenko, 2001) for phytoplankton dynamics and nutrient utilization during the photosynthesis and a non-traditional framework to the description of bacteria dynamics and parameterization of their role in the process of biochemical oxidation of organic matter and nutrient regeneration. The model code has been performed on Fortran-90; duration of 1-month simulation is about 10–12 h on “Alpha-500” Workstation. 4. Data and calibration Data of ecological monitoring at Cartagena Bay, measured since 1996 by the Centre for Oceano-

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

graphic and Hydrographic Research (CIOH, 1997), were used. The monitoring included from two to five seasonal measurements at 23 sites in the bay. The measured parameters were temperature and salinity, transparency, nutrient and oxygen contents, BOD and chlorophyll “a”. Some samples were obtained for phytoplankton primary production, bacterioplankton biomass and bottom demand of oxygen. Information about the sources of contamination was offered by the environment organisations of Cartagena. The scheme of source location is presented in Fig. 1, and the source characteristics are in Tables 1 and 2. A finite-difference grid of 250 m × 250 m, horizontal turbulence of 0.1 m2 s−1 , hydrodynamic time step of 12 s and biogeochemical time step of 1 h were

99

established for the model. Climatic winds, Dique channel flow rates, water transparency and meteorological parameters for calculation of the photoactive solar radiation were used from the observed data. Tidal levels on the marine boundaries were calculated by harmonic analysis. To select the grid parameters, the following requirements have been taken into account: (a) the necessity to include the adjacent part of the sea to calculate the water exchange between the sea and the bay; (b) the need to describe the spatial distributions of the measured elements; (c) numerical scheme stability; (d) optimal duration of calculations. Estimation of the chemical–biological parameters was carried out on the basis of the information from LGMI (1979), Parsons et al. (1984), LO GOIN (1987), Fasham et al. (1990), DHI (1994), Tufford and

Table 1 Industrial waters entering in Cartagena Bay Source Flow rate (m3 per day)

DBO5 (kg per day)

NH4 (kg per day)

NO2 (kg per day)

NO3 (kg per day)

Norg (kg per day)

PO4 (kg per day)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

10.25 3171.0 98.4 0.114 20.30 102.0 0.01 0 204.23 – 0.869 0.31 6.20 1393.1 1.10 1.43 0.01 0.02 0.02 130.34 13.15 564.05 10.21 0.09 11.45 – 60.31 29.35 213.93

8.47 225.0 7.09 0.01 – – – 9.90 – – 0.01 0.01 0.01 1483.5 0.13 0.10 0.01 – – 83.46 2.31 – 0.039 0.06 0.07 0.10 0.30 – 0.32

0.29 – 0.01 0.01 – – – – – – 0.01 0.02 0.01 – – 0 0.01 – – 0.01 – – 0.01 – – – – – 0.03

6.62 1.74 0 0.01 – – – 13.00 – – 0.01 0.022 0.13 – – 0.068 0.01 – – 0.01 – – 0.01 – – – – – 0.02

22.91 349.40 36.08 0.052 0.59 – – – – 0.01 0.025 0.082 0.074 1475.9 0.25 – 0.01 – – – – 14.80 0.47 0.07 3.04 4.30 0.37 0.59 –

– 34.11 3.80 – – – – – – – 0.013 0.013 0.102 – – – 0.01 – – – – – – 0.01 – – 0.39 – –

6042.2

1820.9

0.41

21.65

1909.0

38.45

Total

559.60 400.00 800.00 20.00 55.00 87.40 2.80 984285.0 44.80 0.014 8.33 35.55 68.00 121089.0 71.30 930.0 253200.0 17.30 6.50 404.80 52.60 755.13 19.20 33.20 565.92 172.80 57.70 54.00 316.00 1364112.0

Information was obtained from Corporation of the Dique channel (CARDIQUE, 1996).

PTot (kg per day) – 73.36 5.84 0.019 0.24 – – – – 0.01 0.084 0.022 0.13 48.44 0.15 – 0.01 – – – – 0.54 0.02 0.19 5.60 – 0.50 0.17 – 135.33

100

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

Table 2 Domestic waters entering in Cartagena Bay Source

Flow rate (m3 per day)

1 2 3 4 5 6 7 10 11 12 13 14 16 17 18 Open Channels

40 3500 300 20 200 150 40 25 150 100 25000 1000 11000 2000 2000 15000

Total

60525

DBO5 (kg per day) 13.25 1159.2 99.36 6.62 66.24 49.68 13.25 8.28 49.68 33.12 8280.0 331.2 3643.2 662.4 662.4 1129.5 16207.0

NH4 (kg per day)

NO2 (kg per day)

PO4 (kg per day)

0.54 47.14 4.04 0.27 2.69 2.02 0.54 0.34 2.02 1.35 336.75 13.47 148.17 26.94 26.94 110.85

0.001 0.098 0.008 0.001 0.006 0.004 0.001 0.001 0.004 0.003 0.70 0.028 0.308 0.056 0.056 0.420

0.37 32.44 2.78 0.18 1.85 1.39 0.37 0.23 1.39 0.93 231.75 9.27 101.97 18.54 18.54 17.4

724.1

1.70

439.4

Information was obtained from Corporation of the Dique channel (CARDIQUE, 1996).

McKellar(1999), laboratory and field experiments (CIOH, 1997). Some parameters (Vfmax ,Vbmax , δ, νN1 , νN2 ) were adapted to the temperature conditions in situ (water temperature varies during the year within the limits of 28–30 ◦ C) utilising empirical relationships from the above-mentioned literature. The selection of those parameters was based, on one side, on the best agreement between the model dynamics and the measured data and, on the other side, on reasonable agreement with values obtained by other investigators. The first step of model calibration was a simplified one, i.e. it was carried out on the one-dimensional version (model of annual dynamics of the vertically distributed components of the ecosystem, taking into account the external fluxes of matter and energy), obtained from (1) to (7), cancelling the horizontal terms as a result of spatial integration. Numerical simulations with the one-dimensional model were carried out in two stages. The first stage aimed to achieve a stationary vertical distribution corresponding to the measured data and obtained for fixed external forces (mean annual conditions). This aim was achieved through some adjustments of the original values of biogeochemical constants (within permissible limits), known from the bibliographic sources, mentioned before. In the second stage, annual behavior of radiation, wind, channel fluxes and salinity distribution in the

bay were specified. A secondary correction of the constants was carried out, comparing the model and real seasonal dynamics of the ecosystem variables. These values of the constants were used later in the three-dimensional version of the model. Preliminary use of one-dimensional model was related with the following circumstances: the onedimensional variant requires less computer time than the three-dimensional to adjust the parameters and calibrate the model; the amount of numerical experiments for different sets of parameters and for several years of the model time may be much more for the one-dimensional model. It is simpler to achieve a production–destruction balance through the parameter values selection in their permissible limits, established in the available literature. Calibration of the three-dimensional version was also carried out. Typical distributions of the elements for the wet and dry seasons, obtained from the model, were compared with the observed data. Herein, the “typical” (characteristic) distribution of the chemical– biological elements is one that corresponds to a stationary state of the ecosystem without external flux variability (anthropogenic fluxes and hydrometeorological conditions) for given parameters of the system. Analysis of the biogeochemical block sensitivity over the variation of the internal parameters has shown that the principal (influential) ones are the specific

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

101

Fig. 4. Phytoplankton biomass (I) (in mg ChlA m−3 ) and phosphates (II) (in mg P l−1 ) at the surface for the dry season (January), obtained by field measurements (a) and by the model (b). The mapping was carried out in Surfer (Golden Software Inc.) using the inverse distance to a power gridding method.

102

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

Table 3 Chemical–biological constants in the eutrophication model for Cartagena Bay, established as a result of its calibration Symbol

Value

Units

Vfmax VBmax CkPO4

4.0 2.0 (4.0)a 0.01 0.073 0.1 0.5 110.0 1.0 0.9 0.1 1.0 0.4 1.0 0.16 3.0 0.33 0.4 (0.8)a 10.0 40.0 0.024 0.176 0.034 0.51 2.67 3.4 1.1 0.1 0.1 128.0 0.66 22.0 11.5

per day per day mg P l−1 mg N l−1 per day per day W m−2 m per day – per day m per day – mg O2 l−1 per day per day – – mg O2 l−1 mg C mg−1 ChlA mg P mg−1 C mg N mg−1 C mg P mg−1 C mg N mg−1 C mg O2 mg−1 C mg O2 mg−1 N mg O2 mg−1 N per day per day mg O2 m−2 h−1 – l m−2 h−1 l m−2 h−1

CkN γf µf Iopt w gf ηf δ wgd φ CkO2 νN 1 νN2 θ ␭f Bkorg β C/ChlA β P/C β N/C ant βP/C ant βN/C βO2 /C βO2 /N1 βO2 /N2 νDN νfoto a b γe γi a

Values used for one-dimensional version of the model.

rates of biological production and mortality, and also the gravitational deposition velocity of phytoplankton and detritus. Some results of the verification are presented in Fig. 4. The constants of the chemical–biological block, obtained from the experimental and bibliographic sources and during the one- and three-dimensional model calibration, are presented in Table 3. 5. Model results The final version of the model was used for an environmental improvement efficiency study, especially improvement of the oxygen regime of the

bay. Oxygen absorption in the bottom sediments has been corrected, taking into consideration the organic new (x, y)/F act , modifyflux changes f (x, y) = Forg med new ing (20), where Forg is organic matter flux for the prognostic situation, obtained by the model as a result of gravitational deposition of phytoplankton and detritus with respective velocities wgf and wgd (see Table 3). For this purpose, each variant of simulations was divided into two steps. In the first step, the anthropogenic source characteristics were altered and a new new was calculated field of the organic matter flux Forg bot at given fixed oxygen flux QO2 , corresponding to the actual conditions. Then, a new field of the function

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

f(x, y) was calculated and the simulation was repeated again for the new Qbot O2 . According to the plan of Cartagena development, it is supposed that by 2025 all urban discharges will be conveyed through a central channelisation system

103

to the open sea. It is also planned to reduce the nutrient discharge from the industrial sources by 80% of the current level. The simulation of this alternative has shown that the oxygen deficit on the bottom layer of the internal bay will disappear (see Fig. 5a and b).

Fig. 5. Oxygen distribution (in ml l−1 ) at the bottom layer for the wet season, obtained by the model for (a) the actual state, (b) elimination of 100% of domestic and 80% of industrial contamination, (c) the last variant under limitation of the Dique channel discharge down to 50 m3 s−1 , (d) the same variant under limitation by one-half of nutrient transport by the channel.

104

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

For the external bay, the deficit will be diminished, but the situation will not be changed qualitatively, because the major area of this layer will still contain less than 1.5 ml l−1 of oxygen. The explanation of this simulation result is as follows: diminishing the anthropogenic source intensity does not eliminate the main reason of the oxygen deficit—the Dique channel input of fresh water with a high nutrient content, a strong pycnocline formation in the bay and a weak ventilation of the bottom layer. Therefore, in addition to the one mentioned before, the following three alternatives were considered: 1. Reducing the water discharge from the channel for the wet season from 150 to 50 m3 s−1 (the last value is a flow rate for the dry season); 2. The same alternative, but with additional reduction of mineral nitrogen and phosphorus by one-half in the channel water; 3. Channel closure. By reducing the channel flux to 50 m3 s−1 , the suspended particle discharge is decreased and the transparency in the bay is increased. This effect in the model is taken into account by applying the dry season’s transparency field. This causes an increase

of the phytoplankton primary production and the flux of dead organic matter to the sediments in the central part of the bay, where the oxygen minimum is increased (see Fig. 5c). Thus, the partial elimination of the channel water discharge is not a favourable alternative to improve the oxygen regime of the bay, because the channel is still a strong source of nutrients, and the zones of maximum phytoplankton production simply move from the bounds to the center of the bay. Nutrient reduction in the channel with a water flux of 50 m3 s−1 (Alternative 2), gives a better result. The phytoplankton and bacterioplankton biomass and the organic matter flux are decreased significantly in the major part of the bay. As a consequence, the oxygen deficit will be less in the bottom layer (see Fig. 5d). For the entire channel closure, the eutrophic status of the bay tends to the open sea regime. The oxygen deficit disappears (Fig. 6b). In the southern part of the external bay, which has a slow oceanic water renewal and where the major group of the industrial sources is situated, the oxygen in the bottom layer for the wet season will be more than 1.5 ml l−1 , while in the internal bay it will be more than 2.5 ml l−1 . The above-mentioned alternatives of the Dique channel limitation, or its entire closure, are rather real

Fig. 6. Oxygen distribution (in ml l−1 ) at the bottom layer for the wet season, obtained by the model under entire closure of the channel without anthropogenic source limitation (a) and 100% domestic and 80% industrial source elimination (b).

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

105

ones. At this moment the channel has lost its navigational importance because of other connections with the sea apart of the bay. Moreover, there is a primary technical decision to construct a dam to limit the water discharge from the channel to Cartagena Bay and redistribute this water through the other branches. Also, it is planned to unify the fluvial system and the channel in one network. The contaminated water of the Magdalena River will be filtered from nutrients passing through this system. Finally, the main condition of the oxygen regime improvement in the bay is the control of anthropogenic influence within the above-mentioned limits. Simulation of the channel closure with conservation of the actual level of anthropogenic sources has shown that the situation will unexpectedly be worse (see Fig. 6a). The explanation is that the channel closure will cause the transparency to increase in the bay up to 6 m, and this factor will not be a limitation of photosynthesis, while the amount of nutrients from the anthropogenic sources will be still sufficient to maintain a high level of phytoplankton production.

tial distribution of bacteria biomass affects biochemical oxidation of organic matter, nutrient regeneration and oxygen consumption. In this respect, the present structure of biogeochemical block is more complete and systematic than that usually used in engineering practice in eutrophication models (Ambrose et al., 1993; DHI, 1994). In the cited models the organic matter oxidation rate is given as a constant, independent of functional characteristics and available biomass of bacteria. The model presented in this paper also takes into account the difference between the chemical compounds of the natural and anthropogenic organic matter in an explicit form. Mathematical modelling of distinct ecological situations for Cartagena Bay has shown that an effective strategic objective must be to reduce water and nutrient discharge of the Dique channel as much as possible. Simultaneously, anthropogenic impacts from industrial and domestic sources should be regulated. The model has been verified and could be used for analogous problems in other deep marine bays and gulfs.

6. Conclusions

Acknowledgements

The oxygen regime of Cartagena Bay is formed by a complex of natural and anthropogenic factors, some of which are: basin morphology, Dique channel discharge, water exchange with sea, wind regime and anthropogenic impact. Therefore, it is impossible to define the efficiency of environment-improvement strategies without mathematical modelling. Taking MECCA (Hess, 1985) as a base of the presented hydrodynamic model, an original three-dimensional eutrophication model has been developed, including both transport and biogeochemical blocks. The transport block is based on transport quasi-conservative finite-difference numerical schemes (Lonin, 1997) in agreement with those used for the hydrodynamic equations. The biogeochemical block of the model describes local chemical–biological processes, responsible to the non-conservativity of the modelled components of the marine ecosystem. An important feature of the last block is that the role of bacteria in organic matter mineralisation is explicitly taken into account. This explicitness allows to consider how non-uniform spa-

The paper is based on works supported by DIMAR, Colombian Navy and Colciencias. The authors are grateful to their colleagues at the CIOH; special thanks to Mr. L. Calero for scientific support and Mr. C. Tejada for help in English translation. References Ambrose, R.B., Wool, T.A., Martin, J.L., 1993. The water quality analysis simulation program. In: WASP5. Part A: Model Documentation. Environmental Research Laboratory, Athens, GA, USA, 210 pp. Boris, J.P., Book, D.L., 1976. Flux-corrected transport. III. Minimal-error FCT algorithms. J. Comp. Phys. , 397–431. Corporation of the Dique channel (CARDIQUE), 1996. Informe Preliminar sobre Contaminación de fuentes industriales a la Bah´ıa de Cartagena. Cartagena, Colombia (in Spanish). Cerco, C.F., Cole, T., 1995. User’s Guide to the CE-QUAL-ICM 1995 Three-Dimensional Eutrophication Model. US Army Corps of Engineers, Waterways Experiment Station, Washington, DC, 3-1-3.54. Centre for Oceanographic and Hydrographic Research (CIOH), 1997. Estudios que identifican la condición ambiental del ecosistema, grado de impacto sobre los componentes biológicos, potencial de la capacidad de recuperación y acciones de

106

Y.S. Tuchkovenko, S.A. Lonin / Ecological Modelling 165 (2003) 91–106

reabilitación del sistema de la Bah´ıa de Cartagena. Final Report. Project GEF/RLA/93/G41, UNOPS, PNUD, Cartagena, Colombia, 315 pp. (in Spanish). Danish Hydraulic Institute (DHI), 1994. User Guide and Reference Manual Water Quality Module, Release 2–4 MIKE21, Denmark, 14 pp. Fasham, M.J.R., Ducklow, H.W., McKelvie, S.M., 1990. A nitrogen-based model of plankton dynamics in the oceanic mixed layer. J. Mar. Res. , 591–639. Fletcher, C.A.J., 1991. Computational Techniques for Fluid Dynamics. 2. Specific Techniques for Differential Flow Categories. Springer, Berlin. Hess, K.W., 1985. Assessment model for estuarine circulation and salinity. NOAA Technical Memorandum NESDIS AISC 3—National Environmental Satellite, Data, and Information Service, USA, Washington, DC, 39 pp. Institute of Biology of the Southern Seas (INBUM), 1991. Modeling of the processes of water purification in the shelf zone. Ukraine, Sevastopol, 227 pp. (in Russian). Intitute of Oceanology of Russian Academy of Sciences (IOAN), 1989. Models of Oceanic Processes. Moscow, 366 pp. (in Russian). Leningrad Hydrometeorological Institute (LGMI), 1979. Modelling of transport processes and transformation of matter at sea. Leningrad, Gidrometeoizdat, 290 pp. (in Russian). Leningrad Branch of the State Oceanographic Institute (LO GOIN), 1987. Modelling of ecosystem components. In: Problems of Studying and Mathematical Modelling of the Baltic Sea

Ecosystem, vol. 3. Leningrad, Gidrometeoizdat, 255 pp. (in Russian). Lonin, S.A., 1997. Cálculo de la Transparencia del agua en la Bah´ıa de Cartagena. Bolet´ın Cient´ıfico CIOH—Centro de Investigaciones Oceanográficas e Hidrográficas, Cartagena, Colombia, N.18, pp. 85–92 (in Spanish). Lonin, S.A., Tuchkovenko, Y.S., 2001. Water quality modelling for the ecosystem of the Cienaga de Tesca coastal lagoon. Ecol. Model. , 279–293. Lyakhin, Y.I., 1980. About oxygen exchange intensity between the ocean and the atmosphere. Oceanologia (6), 1014–1021 (in Russian). Monin, A.S., Yaglom, A.M., 1971. Statistical Fluid Mechanics, vol. 1. MIT Press, Cambridge, MA. Parsons, T.R., Takahashi M., Hargrave, H., 1984. Biological oceanographic processes. Pergamon Press, Oxford, New York, 324 pp. Raymond, J., 1983. Plankton and Ocean Productivity: Phytoplankton, vol. 1. Moscow, 567 pp. (in Russian). Sarmiento, J.L., Slater, R.D., Fashman, M.J.R., Ducklow, H.W., Toggweiler, J.R., Evans, G.T., 1993. A seasonal three-dimensional ecosystem model of nitrogen cycling in the North Atlantic euphotic zone. Global Biogeochem. Cycles (2), 417–450. Tufford, D.L., McKellar, H.N., 1999. Spatial and temporal hydrodynamic and water quality modeling analysis of large reservoir on the South Carolina (USA) coastal plain. Ecol. Model. , 137–173.