Applications of a new ecosystem model to study the dynamics of phytoplankton and nutrients in the Ariake Sea, west coast of Kyushu, Japan

Journal of Marine Systems 75 (2009) 1–16

Contents lists available at ScienceDirect

Journal of Marine Systems j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j m a r s y s

Applications of a new ecosystem model to study the dynamics of phytoplankton and nutrients in the Ariake Sea, west coast of Kyushu, Japan Nguyen Thi Minh Hang a, Nguyen Cao Don b,c,⁎, Hiroyuki Araki b, Hiroyuki Yamanishi b, Kenichi Koga d a b c d

Department of the Environment, Water Resources University, 175 Tay Son, Dong Da, Hanoi, Vietnam Institute of Lowland Technology, Saga University, Saga 840-8502, Japan National Institute for Environmental Studies, 16-2 Onogawa, Tsukuba-City, Ibaraki, 305-8506 Japan Department of Civil Engineering, Saga University, Saga 840-8502, Japan

a r t i c l e

i n f o

Article history: Received 28 August 2007 Received in revised form 30 June 2008 Accepted 11 July 2008 Available online 23 July 2008 Keywords: Ariake Sea Estuary Tidal flat Flooding and drying New ecosystem model Hydrodynamic model Phytoplankton Nutrients Dissolved oxygen Suspended solids

a b s t r a c t In this study, we present the development and application of a new ecosystem model coupled with a hydrodynamic model to describe the important physical, chemical and biological processes of an ecosystem in the marine environment, the Ariake Sea in the west coast of Kyushu, Japan. The model was calibrated and validated using in-situ field measurements from various monitoring stations in the sea. The presented results covered the period from January 1991 to December 2000. The results showed that chlorophyll-a, nutrients and dissolved oxygen levels varied seasonally in response to weather and boundary condition. Through this study, the model was shown to be able to handle the flooding and drying processes that usually exist and play an essential role over the estuarine-tidal flats of the sea. © 2008 Elsevier B.V. All rights reserved.

1. Introduction Shallow waters in coastal zones are complex systems characterized by large fluctuations in the physical and chemical parameters. Recently, much attention has been paid to the functions and effects of shallow waters on the water quality of estuarine and tidal flat areas, a mixing region between rivers and ocean waters. The mixing of fresh and saline waters in relatively shallow depths in those areas can create not only an extremely productive ecosystem but also a region having high concentrations of suspended sediment that can limit light penetration and decrease algal production. Thus, suspended sediment can affect water quality, other productions and the aesthetics of coastal areas. ⁎ Corresponding author. Institute of Lowland Technology, Saga University, Saga 840-8502, Japan. Tel./fax: +81 952 288571. E-mail address: [email protected] (N. Cao Don). 0924-7963/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jmarsys.2008.07.006

In Japan, shallow waters in estuarine and tidal flat areas are the most valuable ecosystem related to lowland and urbanized areas that are often convenient and attractive locations for human settlement and their activities. The estuarine and tidal flat areas are habitation and spawning areas for various aquatic lives and are important in maintaining sustainable ecosystem and buffering against the impact of red tide and oxygen depletion (Suzuki et al., 1997). Low dissolved oxygen (DO) in the water body directly affects survivals of fishes and migrations of higher organisms, thus alters estuarine healthy ecological condition (Hang, 2007). In the west coast of Kyushu Island, the Ariake Sea is a typical shallow sea in Japan as schematically shown in Fig. 1. The sea is a complex system in terms of its biological ecosystem and hydrodynamics. The sea is a tidally-forced, semi-enclosed estuary encompassed by wide-spreading intertidal areas comprising mudflats and salt marshes, where flooding and draining are predominant. Its total area

2

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

Fig. 1. Map showing the Ariake Sea, west coast of Japan.

is about 1700 km2 and its length is about 100 km along the gulf axis. On average, its width is 16 km and water depth is about 20 m. The sea has the largest tidal range in Japan, about 6 m at a spring tide in the inner part. A tidal flat, which occurs there, is also the largest with an area of 207 km2 or about 40% of the total tidal flats in Japan. The Ariake Sea is particularly of a great importance in its uniqueness of ecological and biological characteristics. Its estuarine and tidal flat areas serve as conduits through which nutrients (i.e., nitrite, nitrate, ammonium, organic nitrogen and phosphorus) pass into productive coastal waters. They also provide important spawning and nursery grounds for many species of finfish and shellfish. The sea becomes famous for its rich fishery products and Porphyra Yezoensis sp. (Japanese Nori) cultivation. Previously, the sea used to provide a large amount of seaweeds and other fishery products and it was called “the sea of treasure”. In recent decades, however, due to environmental change and other influencing factors, there has been a rapid reduction in the fishery products. In particular, the product of Atrina Pectinata Japonica (Japanese Tairagi oyster), the principal product of the Ariake Sea, living in the upper 0.10–0.15 m of the mud bed, declined from 13.4 million tons in 1976 to 79.0 tons in 1999. That of Crassostrea gigas (Japanese Kaki oyster), living near the surface

mud, dropped from 0.80 million tons in 1979 to 0.12 million tons in 1999. Moreover, the product of Sinonovacula constricta (Japanese Agemaki oyster), living in mud or muddy-sand area of comparatively high ground, was even more worse with 17.0 tons harvested in 1976 but dropped to nil by 1992. Fig. 2 illustrates the fishery products harvested in the last three decades. It shows that, since 1984, the fishery product has been decreased dramatically. Moreover, the number of species living in the mudflats was also reduced. In the past (from 1960 to 1990), some species such as Asari, Tairagi, Agemaki, Sarubo could be found in the tidal flats, however, they all disappeared after the year 1990. Some are hard to find nowadays. So far, such habitat loss has occurred in within the sea as a result of direct and indirect impacts. The reason of such a reduction of the fishery products has been unclear. There are probable causes being argued for the declination of the fishery product, which are: (1) the aggravation of water quality due to the exchange of seawater with fresh water from surrounding rivers and land areas, resulting in abnormal amount of phytoplankton that thrives on sun, water and nutrients, leading to a seaweed scarcity; (2) sediment transport and resuspension that can affect planktonic primary productivity because it can limit light penetration and

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

3

Fig. 2. Rapid reduction of fishery product during the last two decades.

nutrient availability; (3) other causes resulting in high sulphide content in the bottom sediment and DO depletion in seawater. All these factors are responsible for important ecosystem alterations characterized by algal blooms, oxygen depletion and hydrogen sulphide production. Additional problems arise from coastal erosion, land subsidence, groundwater flow and effects related to extreme meteorological events. Those causes have become the environmental issues for public concern and argument. According to salinity level, the sea is often divided into three regions, which are the inner part, the central part and the lower part, as shown in Fig. 1. The salinity, nutrients and chlorophyll-a levels are quite different in the three regions. According to observations, the salinity concentration is lowest in the inner part and highest in the lower part. In contrast, nutrients and phytoplankton concentrations in the inner part are much higher than those in the other regions. Moreover, it is observed that the spatial pattern of phytoplankton is closely related to the distribution of nutrients. Thus, much effort has been spent to study the water environment in the Ariake Sea, through field survey, experiment and data analysis, for example the work of Yamanishi et al. (2002, 2005). On the other hand, there have been a few attempts to model the dynamics of the sea, in such a way that the important physical, chemical and biological processes are included. Because those processes are very complicated, studies on water quality in the sea usually rely on water quality models, which become a powerful tool in determining and predicting the trend of dynamics of water quality constituents. As mentioned earlier, because the tidal difference in the sea is very great, currents and waves highly affect the material exchanging process in the water column. Moreover, because the flooding and drying processes exist and play an essential role over the estuarine-tidal flats of the sea, studies of the water quality in the sea must take this process into account. There are numbers of ecological models constructed and used to simulate natural ecological systems, for example the works of Fransz and Verhagen (1985), Franks et al. (1986), Frost (1987), Fasham et al. (1990), Taylor et al. (1993), Sharples and Tett (1994), McCreary et al. (1996), Drago et al. (2001), and others. However, such models were developed for describing eutrophication in rivers, lakes and reservoirs. On

the other hand, several models were developed for specific coastal waters, such as an eco-hydrodynamic model (Fennel, 1995; Neumann, 2000; Neumann, Fennel and Kremp, 2002; Fennel and Neumann, 2004) developed for the Baltic Sea, an eutrophication model (Cerco and Cole, 1993, 1995) developed for Chesapeake Bay. The chemical–biological model described in Fennel (1995) simulates a nitrogen cycle. This model was linked with a circulation model, the Modular Ocean Model MOM2.2 (Pacanowski et al., 1990) to estimate budgets and nitrogen transports. The model is considered an important step towards an ecosystem model of the Baltic Sea. CE-QUALICM (Cerco and Cole, 1995) is a finite volume eutrophication model that can be used to calculate 1-, 2-, 3-D water quality variables. This model incorporates 22 state variables that include multiple forms of algae, carbon, N, P, and silica; and DO. The model incorporates a predictive submodel of benthic processes including sediment oxygen demand and sedimentwater nutrient flux. Flows, diffusion coefficients, and volume must be specified externally and read into the model. On the other hand, in Japan, most of the numerical models for coastal marine ecosystem that have been used for environmental assessment are composed only of the pelagic system, while the benthic variables are given as boundary conditions. As a first step towards the management perspective and ecological restoration for the Ariake Sea, in this paper, we present a new ecosystem model that has been developed and validated for identifying and quantifying how the physical, biological and chemical processes and their interactions control the spatial distribution of water quality constituents in the sea. The model has served as a basic model for studies of water environment in the Ariake Sea (Don et al., 2005). 2. Model description We applied a hydrodynamic (HD) model (non-hydrostatics) to couple with our new ecosystem model configured for the Ariake Sea (Don et al., 2005, 2007; Hang et al., 2007). The HD model adapted in this study is MIKE-3 (DHI, 2003), which consists of two sub-modes: a HD module linked with an advection–dispersion (AD) model. The HD simulation is performed to produce HD data such as current velocities and

4

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

water elevations. To handle the transport of dissolved or suspended substances, nutrients and plankton, the AD simulation is performed based on the HD data. Thereafter, the ecosystem model (Don et al., 2005) facilitates integration of physical conditions with nutritional and biological dynamics, which describes pelagic dynamics with special consideration of phytoplankton primary production. 2.1. Hydrodynamic and advective-dispersion models The hydrodynamic model is a general numerical modeling system for the simulation of water levels and flows in estuaries, bays and coastal areas. It simulates unsteady three-dimensional flows and has been applied in a large number of studies. Hydrodynamic equations consist of the continuity equation (Eq. (1)) and equation of motion (Eq. (2)): @ρ @ρuj þ ¼0 @t @xj

ð1Þ


 @ui @ ui uj 1 @ρgi z þ 2Xij uj ¼ − þ þ gi ρ @xi @x @t     @ @ui @uj 2 þ νT þ − δij k þ ui R @xj 3 @xj @xi

ð2Þ

2.2. Ecosystem model

where ρ — the local density of the fluid, cs — the speed of sound in seawater, ui — the velocity in the xi-direction, Ωij — the Coriolis tensor, gi — the gravitational vector, νT — the turbulent eddy viscosity, δ — Kronecker's delta, k — the turbulent kinetic energy, R — source-sink terms, and t — time. The advection–dispersion model solves the advection– dispersion equation for dissolved or suspended substances, including salinity and temperature, in three dimensions:   @C @ ðui C Þ @ @C þ ¼ Di þ I0 þ E0 @t @xi @xi @xi

associated dispersion coefficients; I0 is the function that represents the internal source or sink of the water quality components; and E0 is the external loading from point sources (river discharge) and non-point sources of the water quality components. The information on water velocity components at each time step is provided by the HD model and is assumed constant during each time integration of the advection– dispersion equation. The advection–dispersion equation is solved at each time step following the time integration of the hydrodynamic equations (Eqs. (1) and (2)). The AD model offers a choice of turbulence models in terms of an eddy viscosity and a bed shear stress. The eddy viscosity can dynamically be specified by means of the mixed k–ε/Smagorinsky formulation with a standard k–ε model in the vertical and a Smagorinsky formulation in the horizontal. This HD model uses an alternating direction implicit technique to integrate these equations in the time–space domain. The equation matrices that result for each direction and each individual grid line are then resolved by a double sweep algorithm.

ð3Þ

where C is the concentration of the water quality components; u i is the water velocity components; D i is the

The model structure is based on a conceptual model of the estuary nutrient cycling processes. This conceptual model is developed based on the results of the initial field investigation results and on the observed conditions in the sea estuaries. Refinement of the conceptual model will be an integral part of the model calibration and validation. The ecosystem model describes the main interactions among water, phytoplankton, zooplankton, nutrients, detritus and superficial sediment. The model has eleven (11) state variables in the water column, which are phytoplankton, zooplankton, detritus, nitrogen, phosphorus, dissolved oxygen, organic, inorganic sediment, etc., (as listed in Appendix A). A schematic review of the model is seen in Fig. 3 that depicts the state variables and primary producers that are

Fig. 3. Main interactions among the state variables in the ecosystem model.

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

represented by a single state variable, phytoplankton. In this model, phytoplankton was divided into three state variables: diatoms, green algae and blue-green algae. The growth of these three groups of phytoplankton is governed by photosynthesis, the uptake of inorganic nitrogen and phosphorus, their decay to detritus and inorganic compounds with oxygen consumption and mortality and is affected by the grazing processes. Phytoplankton production is controlled by a description of the intracellular concentration of nitrogen and phosphorus, light availability, temperatures and it is possible to switch between different dominating phytoplankton species (Zhang et al., 2004). In the model, the photosynthetic oxygen production and the autotrophic respiration vary with the water depth due to the light dependency of the autotrophs. The depth variation was modeled using the Lambert Beer Law, which requires a light extinction coefficient for the water to describe the dampening of the light irradiation through the water column. Lambert Beer Law or the light dampening function: F(H) = e− kH, where k is light extinction coefficient (m− 1), and H is water depth (m). The modified Lambert Beer function is listed in Table 2, Eqs. (16)–(18). Zooplankton grazing plays an important role in phytoplankton productivity in the estuarine system. A residual part of the phytoplankton in filter-feeding predation is directly transformed to detritus. Detritus settles in the sediment and decomposes to dissolved organics, which is then degraded to inorganic matter, thus consuming oxygen. Dissolved oxygen was modeled based on the re-aeration at the water–air interface, primary production and respiration, nitrification, degradation of detritus and sediment oxygen demand. The controlling factors for growth and death processes of benthic vegetation are incorporated in the model structure and the interaction with the other components is an integral description of the model. Sources and sinks of diatoms and other algae include production, respiration, zooplankton grazing, settling, resuspension and death. These are expressed in Eqs. (4)–(6) (see Table 1). The production of phytoplankton depends on the intensity of light and silica, the availability of nutrients and the ambient temperature. The gross production is computed with a multiplicative approach considering the maximal rate of production, the influence of light, temperature and the internal concentrations of nitrogen and phosphorus. Because no data on silica were available, the dependence on silica was given in the growth rate, k, of phytoplankton. In this study, although waves were not directly modeled, the effects of resuspension due to currents and waves to phytoplankton were considered by adding a “resuspension” term in each equation, from Eqs. (4)–(6) of Table 1. The detail of the conservative equations describing each biochemical process for phytoplankton production, nutrients cycling, dissolved oxygen and other water quality constituents is listed in Table 1, following Eqs. (4) –(15). The rate equation for each biochemical process is listed in Table 2. Those conservative equations of the ecosystem model are programmed and solved using an independent equation solver in MIKE-3 that allows the ecosystem model with tailor-made process descriptions to be incorporated directly into the HD and AD models.

5

3. Model setting The coupled simulation model was configured for the Ariake Sea (Don et al., 2005). A finite-difference grid was developed to adequately discretize the model domain. The area of interest was 53.1 × 86.4 km2 and was covered with a 3-D grid. The sizes of each cell were Δx = 900 m, Δy = 900 m. The vertical direction, z, was divided into 12 layers with Δz = 2 m interval. The 2-D grid contained 5664 cells: nx = 59, ny = 96, where ni denotes the number of cells in the i direction. Because the sea is composed of tidal flats that expose to the air after a tidal cycle, it is necessary to specify a drying water depth of 0.2 m and a flooding water depth of 0.3 m. These two depths are used to determine whether a computational point is dried and thus be taken out of the computations, or flooded (it means reentered into the calculations). Basic parameters and HD parameters include the most fundamental specifications such as bathymetries, simulation period, boundary locations, possible sources/sinks, turbulence model, flooding and drying depths, etc. HD boundary conditions, wind and bed friction, initial conditions also need to be included. Since the sea receives freshwater from eight surrounding rivers, namely the Rokkaku, Kase, Chikugo, Yabe, Kikuchi, Shira, Midori and Shiota River (with locations shown in Fig. 1), the nutrient supply from these rivers must be included. Concerning the terrestrial nutrient supply, a logarithmic correlation of the load (L) and the discharge of a river (Q) was analyzed to set up the ecosystem model. Based on monitoring data of the river water quality (Annual report by the River Bureau, Japan), the L–Q correlation was determined for each river. The simulation time covered the period from January 1991 to December 2000, with a one-minute time step. The model setup focused on estimating the different parameters required. The horizontal eddy viscosity was determined through the Smagorinsky concept, while for the vertical direction the eddy viscosity was determined from a onedimensional k–ε model. This model used transport equations for two quantities to describe the turbulent motion: the turbulent kinetic energy, k and the dissipation rate of turbulent kinetic energy, ε. Calibration of the model focused on choosing parameters was performed through trial and error for year 2000. During the calibration, the hydrodynamic and water quality parameters were adjusted until a satisfactory correspondence between model results and observed field data was obtained. Secondly, the parameters of the hydrodynamic and ecosystem models were validated for data from 1991–1999 by keeping the parameter sets obtained from the calibration step unchanged. 4. Results and discussion 4.1. Hydrodynamic model The hydrodynamic model was calibrated against water level and salinity data at three selected stations as marked on Fig. 1, namely S5 in the inner part, K14 in the central part and K4 in the lower part. The vector maps of measured and simulated current velocities for a flood tide were also

6

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

Table 1 Formulation of the biochemical process for the ecosystem model Conservation equation of the biochemical process

Equation

Phytoplankton diatoms: PPDi dPPDi ¼ PrPPDi −DePP Di −GrZP Di −SePP Di þ RePBDi ¼ kPrDi  fPrDi ðT Þ  fPrDi ðIÞ  f ðDINÞ  f ðPO4 Þ  PP Di dt PP 1 kSe kRe PP  PP Di þ  f ðReÞ  Di  PB −kDeDi  fDeDi ðT Þ  PP Di −kGrZ  fZ ðT Þ  Di   ZP − CHL YZ=Chl dH dH CHL

(4)

Phytoplankton green algae: PPGr dPPGr ¼ PrPPGr −DePP Gr −GrZP Gr −SePP Gr þ RePBGr ¼ kPrGr  fPrGr ðT Þ  fPrGr ðIÞ  f ðDIN Þ  f ðPO4 Þ  PP Gr dt PP Gr 1 kSe kRe PP Gr   PP Gr þ  f ðReÞ   PB −kDeGr  fDeGr ðT Þ  PP Gr −kGrZ  fZ ðT Þ   ZP− CHL YZ=Chl dH dH CHL

(5)

Phytoplankton blue-green algae: PPBl dPPBl ¼ PrPPBl −DePP Bl −GrZP Bl −SePP Bl þ RePBBl ¼ kPrBl  fPrBl ðT Þ  fPrBl ðIÞ  f ðDIN Þ  f ðPO4 Þ  PP Bl dt PP 1 kSe kRe PP  PP Bl þ  f ðReÞ  Bl  PB −kDeBl  fDeBl ðT Þ  PP Bl −kGrZ  fZ ðT Þ  Bl   ZP− CHL YZ=Chl dH dH CHL

(6)

Chlorophyll-a: CHL CHL ¼ PPDi þ PPGr þ PPBl

(7)

Zooplankton: ZP dZP ¼ YZ=Chl  ðGrZPDi þ GrZPGr þ GrZPBl Þ−kDeZ  fZ ðT Þ  ZP dt

(8)

Organic detritus: D dD kSe kRe ¼ YD=Chl ðDePPDi þ DePPGr þ DePPBl Þ þ kDeZ  fZ ðT Þ  ZP−kDeg  fDeg ðT Þ  D− Dþ  f ðReÞ  SD dt dH dH

(9)

Ammonia nitrogen: NH4 kf luxNH dNH4 NH 4 ¼ YN=D  kDeg  fDeg ðT Þ  D−YN=Chl ðPrPP Di þ PrPP Gr þ PrPP Bl Þ −k  fN ðT Þ  NH 4 þ ðSNH4 −NH 4 Þ dt NH 4 þ NO3 NiNH dH

(10)

Nitrate nitrogen: NO3 kf luxNO3 dNO3 NO3 ¼ kNiNH  fN ðT Þ  NH4 −kDeNO3  fN ðT Þ  NO3 −YN=Chl ðPrPP Di þ PrPP Gr þ PrPP Bl Þ þ ðSNO3 −NO3 Þ dt NH 4 þ NO3 dH

(11)

Dissolved Inorganic Nitrogen: DIN DIN ¼ NH4 þ NO3

(12)

Phosphoric acid phosphorus: PO4 kf luxPO4 dPO4 ¼ YP=D  kDeg  fDeg ðT Þ  D−YP=Chl ðPrPP Di þ PrPP Gr þ PrPP Bl Þ−kabs  fabs  PO4 þ ðSPO4 −PO4 Þ dt dH

(13)

Dissolved Oxygen: DO dDO ¼ YDO=Chl ðPrPPDi þ PrPPGr þ PrPPBl Þ−YDO=Chl ½kRes  fPrDi ðT Þ  PP Di þ kRes  fPrGr ðT Þ  PP Gr þ kRes  fPrBl ðT Þ  PP Bl  dt kf luxDO ðDO−SDOÞ þ f ðkRa Þ  ð f ðDOCr Þ−DOÞ −YDO=Z  kRes  fZ ðT Þ  ZP−YDO=D  kDeg  fDeg ðT Þ  D− dH

(14)

Inorganic Suspended solid: ISS dISS kSe kRe ¼−  ISS þ  f ðReÞ  ISS dt dH dH

(15)

Suspended solid: SS SS ¼ YD=Chl ðPPDi þ PPGr þ PPBl Þ þ ZP þ D þ ISS

(16)

Chemical oxygen demand: COD dCOD ¼ YDO=Chl ðPPDi þ PPGr þ PPBl Þ þ YDO=Z  ZP þ YDO=D  D dt Note: Notation of state variables is listed in Appendix A.

(17)

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

7

Table 2 Rate equation of each process Rate equation of each process Effect of solar radiation to production of – Diatoms

– Green algae

– Blue-green algae

Expression  fPrDi ðIÞ ¼

   I I exp 1− expð−kI  SS  H Þ IDi IDi

 fPrGr ðIÞ ¼

Equation (18)

   I I exp 1− expð−kI  SS  H Þ IGr IGr

(19)

    I I exp 1− expð−kI  SS  H Þ IBl IBl

fPrBl ðIÞ ¼

(20)

Effect of water temperature to production of – Diatoms

– Green algae

– Blue-green algae

Effect of nutrient to algae production

 fPrDi ðT Þ ¼

   T T exp 1− TDi TDi



   T T exp 1− TGr TGr

(22)

   T T exp 1− TBl TBl

(23)

fPrGr ðT Þ ¼  fPrBl ðT Þ ¼

f ðDINÞ ¼

(21)

DIN ; HSN þ DIN

f ðPO4 Þ ¼

PO4 HSP þ PO4

(24)

Effect of water temperature to death of θDeDi temperature coefficient for death of diatoms

θDeGr temperature coefficient for death of green algae

θDeBl temperature coefficient for death of blue-green algae

  T ðT−T Þ fDeDi ðT Þ ¼ θDeDiDi ¼ exp 1− Di T

  T ðT−T Þ fDeBl ðT Þ ¼ θDeBlBl ¼ exp 1− Bl T

(27)

  ToGr ¼ exp 1− T

(28)

ðT−ToGr Þ

θZ temperature coefficient for grazing rate

fZ ðT Þ ¼ θZ

θDeg temperature coefficient for degradation rate

  ToDeg ðT−TODeg Þ fDeg ðT Þ ¼ θDeg ¼ exp 1− T

Effect of temperature to absorption of phosphorus

Effect of re-aeration to resuspension

ðT−TON Þ

fN ðT Þ ¼ θN  fabs ¼

 

Dissolved oxygen saturation concentration (g/m3)

Note: Notation of state variables is listed in Appendix A.

(29)

  ToN ¼ exp 1− T ðTb19Þ ðTN19Þ

SISS 0

f ðReÞ ¼ f ðVelÞ ¼

Re-aeration rate constant (1/day)

(26)

  TGr ðT−T Þ fDeGr ðT Þ ¼ θDeGrGr ¼ exp 1− T

Temperature function

θN temperature coefficient for nitrification rate

(25)

(30)



(31)

 ðHbHCr Þ ðHNHCr Þ  ðVELNVELCr Þ ðVELbVELCr Þ

(32)

1 f ðVelÞ

1 0

f ðkRa Þ ¼ 3:93 

VEL0:5 0:728  WIND0:5 −0:371  WIND þ 0:0372  WIND2 þ H H 0:5

f ðDOCr Þ ¼ 14:652−0:0841  SAL ( 0:00256  SAL− 0:41022 þ T  ð0:007991 þT  −0:0000374  SAL−0:000077774  TÞ

(33)

(34)

g

8

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

Fig. 4. Simulated and measured current velocity vectors.

constructed and are shown in Fig. 4. The measured one was assumed to be representative and was used to qualitatively evaluate the flow simulation. Although the measured data points were not dense enough for a direct comparison, the flow direction and magnitude for each time period were found to be similar. The simulated tidal levels also matched well with the observed ones as graphically shown in Fig. 5 for June 2000. The mean absolute error between predictions and observations (MAE) was 4 cm and the root mean square error (RMSE) was 6 cm. These errors were approximately 0.9% of the water depth at the location ST1. A comparison between measured and simulated salinity levels at the three aforementioned stations in Fig. 6 for both calibration and validation shows reasonable agreement

although salinity levels in some periods having high river inflows were underestimated. This figure verifies the ability of the numerical model to reproduce salinity throughout the system. The salinity variation is considerably affected by the amount of fresh water from inland rivers, especially the overall dilution process by freshwater during the rainy seasons. Near the river mouths of the Shiota, Rokkaku, Kase and Chikugo Rivers, salinity levels at a surface layer drop in response to freshwater inflows. Low salinity levels appear in July corresponding to a large amount of fresh water discharge from those rivers. In the tidal flats, there is a small difference in salinity levels in the surface and bottom layer because water in these areas is strongly mixed under currents and waves.

Fig. 5. Simulated and measured tidal levels.

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

9

4.2. Ecosystem model

Fig. 6. Seasonal variation of simulated (lines) and measured (dots) salinity.

Fig. 7 particularly demonstrates the salinity contours during periods of high freshwater inflow from the Kase River. The distribution of freshwater flow and the circulation pattern of the sea result in a gradient of increasing salinity with increasing alongshore distances from the inner to the lower part. Near the river mouths, water column stratification periodically occurs when less dense (warmer or fresher) surface water overlies denser (colder or more saline) bottom waters and predominately exists during the summer months. Such a phenomenon is found to be similar in the other river mouths.

Comparisons between measured and simulated water quality parameters at the three aforementioned stations for both calibration and validation are plotted on the same graph for convenience. The ecosystem model simulates the water quality with respect to nutrients, phytoplankton and dissolved oxygen and other water quality constituents in the water column. The values of rate parameters employed in this study are summarized in Table 3. In this Table, some parameters were fine-tuned during the model calibration to fit observed variations in water quality and some values were adopted from literature. Fig. 8 shows the time variation of COD (Chemical Oxygen Demand), Chl-a (Chlorophyll-a), PO4–P (Phosphoric acid phosphorus) and DIN (Dissolved Inorganic Nitrogen) at three aforementioned stations (S5, K14 and K4). Overall, the simulated results have good agreement with the observed data. COD is supplied from land areas and produced by algal production. There is a small difference between the measured COD and the computed ones. On average, COD is about 4.2 mg/l in the inner part, 1.5 mg/l in the central part and 0.8 mg/l in the lower part. In the inner part, high concentrations of COD in summers are caused by the discharge loading from land areas. For Chl-a, the model produced a good tendency variation of Chl-a, although the surveyed Chl-a data were available only for the year 2000. Chl-a concentration is very low in the lower part but becomes much higher in the inner parts. The highest Chl-a levels appear in the summer months and local maxima are observed in the inner part during July of a year. These results show that, in the inner part, the variation of water quality constituents such as COD, Chl-a, PO4–P and DIN follows a seasonal variation, increasing during the rainy seasons while decreasing during winter. It is because runoff from land areas during the rainy seasons contributes nutrients load and runoff water dilutes seawater in areas near the river mouths. As seen in Fig. 8, nutrients and COD levels are much higher in the inner part compared to those in the other parts, because of resuspension and release processes from the bottom mud. There is a distinct seasonal pattern in the change of nutrients. The concentration of

Fig. 7. Salinity contours during high freshwater inflow (in a vertical section).

10

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

Table 3 Parameters and constants for the ecosystem model Description Phytoplankton and Zooplankton Growth rate of diatoms Growth rate of green algae Growth rate of blue-green algae Optimal solar radiation of diatoms Optimal solar radiation of green algae Optimal solar radiation of blue-green algae Optimal temperature of diatoms Optimal temperature of green algae Optimal temperature of blue-green algae Nitrogen half-saturation constant Phosphorus half-saturation constant Diatoms death rate Green algae death rate Blue-green algae death rate Optimal temperature for grazing Optimal temperature for degradation Optimal temperature for nitrification Grazing rate by zooplankton Mass zooplankton to Chl-a ratio Settling rate Resuspension rate of phytoplankton Phytobenthos in sediment Critical velocity Critical water depth Death rate zooplankton Optical attenuation coefficient Organic detrituss Mass detritus to Chl-a ratio Degradation rate detritus Detritus mass in sediment Dissolved Oxygen DO to Chl-a ratio DO to mass detritus ratio DO to mass zooplankton ratio Respiration rate of phytoplankton Respiration rate of zooplankton Flux rate of DO DO sediment exchange rate NH4–N and NO3–N N to mass detritus ratio N to Chl-a ratio Nitrification rate of NH4–N Flux rate of NH4–N Flux rate of NO3–N NH4–N in sediment Denitrification rate of NO3–N PO4–P P to Chl-a ratio P to mass detritus ratio P to mass zooplankton ratio Absorption rate of PO4–P Flux rate of PO4–P Inorganic SS Mass in sediment PO4–P in sediment

Nomenclature

Value

Unit

Reference

kPrDi kPrGr kPrBl IDi IGr IBl TDi TGr TBl HSN HSP kDeDi kDeGr kDeBl ToGr ToDeg ToN kGrz yz/Chl kSe kRe PB VELCr HCr kDeZ kI

0.8 0.6 0.6 300 350 400 15 25 30 0.006 0.003 0.005 0.005 0.005 20 20 20 0.2 0.1 0.1 0.04–2 0–1000 0.2 0.25 0.025 0.01

day− 1 day− 1 day− 1 E/m2/day E/m2/day E/m2/day °C °C °C mg/l mg/l day− 1 day− 1 day− 1 °C °C °C day− 1 g Z/mg Chl m/day m/day mg/m3 m/s m day− 1 m2/g

Eppley (1972) Eppley and Renger (1974) Calibration Rhee and Gotham (1981) Parson et al. (1990) Sakshaug and Slagstad (1991) Calibration Calibration Calibration Jorgensen (1979) Nakata (1989) Calibration Calibration Calibration Umgiesser et al., 2003 Umgiesser et al. (2003) Umgiesser et al. (2003) Drago et al. (2001) Calibration Jorgensen (1979); Nakata (1989) Calibration Calibration Calibration Calibration Calibration Calibration

yD/Chl kDeg SD

0.1 0.005 200–50000

g D/mg Chl day− 1 mg/m3

Calibration Calibration Calibration

yDO/Chl yDO/D yDO/Z kRes kResZ kfluxDO SDO

0.035 0.35 0.35 0.05 0.05 0.03 0–80

g DO/mg Chl mg DO/mg D mg DO/mg Z day− 1 day− 1 m/day mg/l

Calibration Calibration Calibration Drago et al. (2001) Calibration Calibration Calibration

yN/D yN/Chl kNiNH kfluxNH kfluxNO3 SNH4 kDeNO3

0.15 0.014 0.05 0.02 0.02 10–400 0.09

mg N/mg D g N/mg Chl day− 1 m/day m/day mg/m3 day− 1

Calibration Aminot et al. (1997); Cugier et al. (2005) Umgiesser et al. (2003) Calibration Calibration Calibration Ambrose et al. (1993, 2001); Zheng et al. (2004)

yP/Chl yP/D yP/Z kabs kfluxPO SISS SPO4

0.0012 0.012 0.012 0.1 0.02 0–1000 10–500

g P/mg Chl mg P/mg D mg P/mg Z day− 1 m/day mg/m3 mg/m3

Calibration Calibration Calibration Calibration Calibration Calibration Calibration

nutrients becomes low during the winter–summer periods when the influence of aerobic area in the tidal flat areas increases. It then elevates to a high level during the summer–fall period. This is because the uptake and decomposition of nutrients take place under the aerobic and anaerobic conditions in the mud bed (Araki et al., 2001). The DIN calculated as a combination of ammonium, nitrate and nitrite was also modeled well. Fig. 9 indicates that the concentration of phytoplankton at a specific time is very high in areas near the river estuaries

but becomes low in the other locations. Comparisons of model predictions and observations at particular locations reveal interesting information about the temporal evolution in water quality conditions. Surface water concentrations of phytoplankton are much higher and less variable at the inner part. For dissolved oxygen (DO), it is an important indicator of the tolerance level of the sea to the level of nitrogen and phytoplankton. DO concentration represents the balance between inputs from photosynthesis and from the overlying

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

11

Fig. 8. Seasonal variation of COD, Chl-a, PO4–P and DIN (lines) with the measured data (dots).

atmosphere and outputs due to respiration of animal and plant communities and decaying organic matter (Howes et al., 1999). Good agreement between measured and simulated DO levels at the three aforementioned stations

for both calibration and validation is illustrated in Fig. 10, although some peak values were underestimated. DO concentration decreases during the summer–fall period, then gradually increases during the winter season. In these

12

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

Fig. 9. Distribution of total phytoplankton, at time 14 h:00, in 27 June.

periods, when seawater shows low oxygen levels, it clearly indicates a disruption of the balance due to an overabundance of respiration and decay relative to the amount

Fig. 10. Seasonal variation DO (lines) and measured data (dots).

of oxygen input. In addition, low oxygen levels are directly stressful to living animal and plant communities in the sea (Howes et al., 1999). In the inner part, oxygen content is about 5.5 mg/l in September. On average, the concentration of DO in the sea is between 5.5 to 8.5 mg/l. The spatial distribution of the DO (Fig. 11) indicates that the inner part is sometimes in anoxic or hypoxic conditions as it is close to the point sources (river discharges) of pollutant loads (Fig. 12). In a short time scale, daily dissolved oxygen concentrations are typically lower at night and early morning because of the absence of photosynthesis at night to meet plant and animal respiration. Its concentrations in the inner part show considerable diurnal variability due to changing temperature, winds and daily light levels. In a year, the lowest bottom DO level appears in September. The highest surface DO level appears in February and during this time, DO elevates to around 10.0 mg/l. The salinity and DO stratification generally follow the same pattern of correlating stratification by salinity, temperature and DO. In all cases, bottom levels are badly depleted of oxygen, although not with the same severity as occurred in the summer months (Hang, 2007). For suspended solids, Fig. 13 indicates reasonable agreement between simulated and observed data. The

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

13

Fig. 11. Distribution of dissolved oxygen at time 14h:00, 27 June (surface layer).

suspended solids (SS) concentration is higher in areas near the river mouths but becomes lower in the other locations. Sediment concentration is also rather high near the bottom during initially accelerating tidal flows and

during the final ebb retreat (Yamanishi et al., 2002). The amplitude of the concentration variation increases during spring tides because the larger energy sources input to the system.

Fig. 12. DO contours during a flood tide (in a vertical section).

14

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

inwards as a result of the estuarine circulation processes. This phenomenon results from the tidal phase difference between the offshore and the near-shore waters; the flow direction differs in these two areas, which occurs twice during a tidal cycle. Besides advection and dispersion, other natural processes such as sedimentation, release and biomass production, as well as discharged loads from land areas can contribute effectively to the water quality and mud distribution in the Ariake Sea (Don et al., 2007). 5. Conclusion Fig. 13. Computed and observed SS concentration, year 2000.

Water quality in the sea is closely tied to sediment transport processes because resuspension of sediment increases turbidity and releases stored nutrients (Don et al., 2005). Fig. 14 shows the distribution of total sediment, which includes suspended organic matter and detritus, in the sea after the simulation period. It is apparent that high SS levels appear in the western part of the inner part as a result of current circulation and other estuarine processes. Moreover, in the river mouths, tides, currents and waves can influence some distance from the mouths. In these regions, a thick layer of sediment is deposited during the flood phase of the tide, which is then resuspended and transported towards the sea during the ebbs. During periods of strong winds or during spring tides, this sediment is probably resuspended and transported further

The current study provides an attempt to develop a new ecosystem model describing the important physical, chemical and biological processes of an ecosystem in the marine environment. The model has been considered to be an important step towards an ecosystem model of the Ariake Sea. The model was successfully been calibrated and moderately validated for the Ariake Sea, a semi-closed sea in the west coast of Kyushu, Japan. The model is shown to be able to simulate the spatial and temporal variation of phytoplankton, nutrients and other water quality constituents in response to the variation of boundary and weather condition. The model is found to be able to handle flooding and drying processes that are predominant over the intertidal flats of such a bay. The results show that there is a distinct seasonal pattern in change of nutrients. High concentrations of nutrients, chlorophyll-a and suspended solid appear in the inner bay during summers when high discharge loadings from inland rivers and

Fig. 14. Distribution of total sediment (surface layer).

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

high release and resuspension from the bottom mud are predominant. On the other hand, nutrient levels become low during winter to early summer, which coincide with the period when the influence of aerobic area in the tidal flat increases, and become high during summer to fall because nutrients consumption is accelerated under the condition of aerobic and anoxic/ anaerobic in the mud bed. Acknowledgements The revised version of this paper was prepared during the professorship of the second author at the Institute of Lowland Technology (ILT), Saga University, Japan in 2008. Our sincere thanks go to the editor of this journal, Prof. Wolfgang Fennel, who read the entire manuscript and offered valuable suggestions that led to many improvements in this paper. Special thanks to the anonymous reviewers who suggested important improvements to an earlier version of the manuscript.

Appendix A Notation PrPPDi , DePPDi , GrPPDi , Production, death, grazing, settling, resuspension SePPDi, RePPDi of Phytoplankton Diatom, respectively PrPPGr, DePPGr, GrPPGr, Production, death, grazing, settling, resuspension SePPGr, RePPGr of Phytoplankton Green algae, respectively PrPP Bl , DePP Bl , BlPP Bl , Production, death, grazing, settling, resuspension SePPBl, RePPBl of Phytoplankton Blue-green algae, respectively PhDO Photosynthesis ResDO Respiration of DO DegDO Degradation of DO FluxDO Flux of DO Re-air Re-aeration of DO NiNH4 Nitrification DeNO3 Denitrification UpNO3 Uptake of NO3 FluxNO3 Flux of NO3 MiPO4 Mineralization UpPO4 Uptake of PO4 AbsPO4 Absorption of PO4 FluxPO4 Flux of PO4 State variables PPDi (μg/l) PPGr (μg/l) PPBl (mg/l)

Phytoplankton Diatoms Phytoplankton Green algae Phytoplankton Bluegreen algae

NO3 (mg/l)

Nitrate nitrogen

DCOD

Dissolved COD

ISS (mg/l)

Inorganic Suspended Solid Phosphoric acid phosphorus Dissolved Oxygen

ZP (mg/l)

Zooplankton

PO4 (mg/l)

D (mg/l)

Organic detritus

DO (mg/l)

NH4 (mg/l)

Ammonia nitrogen

Weather and flow condition variables Solar radiation I (J/cm2/day)

VEL (m/s)

T (°C)

Water temperature

DIR (degree)

dH (m)

Water depth in model layer Water depth

WIND (m/s) SAL (psu)

H (m)

Horizontal current speed Horizontal current direction Wind speed Salinity

15

References Ambrose Jr., R.B., Wool, T.A., Martin, J.L., 1993. The Water Quality Analysis Simulation Program, WASP5, Part A: Model Documentation. U.S. Environmental Protection Agency, Athens, Georgia. 202 pp. Ambrose, B., Wool, T.A., Martin, J.L. (2001). The Water Quality Analysis Simulation Program, WASP6, User Manual, US EPA, Athens, GA. Aminot, A., Guillaud, J.F., Kerouel, R., 1997. La baie de Seine: hydrologie, nutriments et chlorophylle 1978–1994. Edition IFREMER, Reperes Ocean, vol. 14. 148 pp. Araki, H., Yamanishi, H., Koga, K., Sato, K., 2001. Study on environmental change and peculiarity of the Ariake Sea, Japan. In: Brebbia, C.A., Anagnostopoulos, P., Katisfarakis, K. (Eds.), Water Resources Management. WIT Press, pp. 341–350. Cerco, C.F., Cole, T., 1993. Three-dimensional, eutrophication model of Chesapeake, Bay. Journal of Environmental Engineering 119, 1007. Cerco, C.F., Cole, T. (1995). User's Guide to the CE-QUAL-ICM ThreeDimensional Eutrophication Model, Release Version 1.0, Technical Report EL-95-15, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Cugier, P., Menesguen, A., Guillaud, J.F., 2005. Three-dimensional 3D ecological modelling of the Bay of Seine English Channel, France. Journal of Sea Research 54, 104–124. DHI, 2003. MIKE 3: Scientific Documentation. Estuarine and Coastal Hydraulics and Oceanography Hydrodynamic Module. Don, N.C., Hang, N.T.M., Araki, H., Yamanishi, H., 2005. Modeling of water quality and sediment transport in Ariake Bay, Kyushu, Japan. The 8th International Conference on Cohesive Sediment Transport. Don, N.C., Araki, H., Hang, N.T.M., Yamanishi, H., Koga, K., 2007. Sediment transport and short-term sedimentation processes in the tidal flats of the Ariake Sea, west coast of Kyushu, Japan. Journal of Coastal Research, Special Issue 50, 837–841. Drago, M., Cescon, B., Iovenitti, L., 2001. A three-dimensional numerical model for eutrophication and pollutant transport. Ecological Modelling 145, 17–34. Eppley, R.W., 1972. Temperature and phytoplankton growth in the sea. Fishery Bulletin 70, 1063–1085. Eppley, R.W., Renger, E.M., 1974. Nitrogen assimilation of an oceanic diatom in nitrogen-limited continuous culture. Journal of Phycology 18, 534–551. Fasham, M.J.R., Ducklow, H.W., McKelvie, S.M., 1990. A nitrogen-based model of plankton dynamics in the oceanic mixed layer. Journal of Marine Research 48, 591–639. Fennel, W., 1995. Model of the yearly cycle of nutrients and plankton in the Baltic Sea. Journal of Marine Systems 6, 313–329. Fennel, W., Neumann, T., 2004. Introduction to the Modeling of Marine Ecosystems. Elsevier Oceanographic Series, vol. 72. Amsterdam, 297 pp. Franks, P.J.S., Wroblewski, J.S., Flierl, G.R., 1986. Behavior of a simple plankton model with food-level acclimation by herbivores. Marine Biology 91, 121–129. Fransz, H.G., Verhagen, J.H.G., 1985. Modelling research on the production cycle of phytoplankton in the Southern Bight of the North Sea in relation to riverborne nutrient loads. Netherlands Journal of Sea Research 19, 241–250. Frost, B.W., 1987. Grazing control of phytoplankton stock in the open subarctic Pacific Ocean: a model assessing the role of mesozooplankton, particularly the large calanoid copepods, Neocalanus spp. Marine Ecology. Progress Series 39, 49–68. Hang, N.T.M., 2007. Development and Applications of an Ecosystem Model for Water Environmental Management in the Ariake Sea of Japan. Ph.D. Dissertation, Saga University, Japan, 181 pages. Hang, N.T.M., Araki, H., Don, N.C., Yamanishi, H., Koga, K., 2007. Hydrodynamics and water quality modeling for the ecosystem of the Ariake Sea, Kyushu, Japan. Journal of Coastal Research, Special Issue 50, 800–804. Howes, B., Williams, T., Rasmussen, M., 1999. Baywatchers II. Nutrient related water quality of Buzzards Bay embayments: a synthesis of Baywatchers monitoring 1992–1998. The Coalition for Buzzards Bay. Jorgensen, S.E., 1979. In: Vaerlose, D.K. (Ed.), Handbook of Environmental Data and Ecological Parameters. International Society for Ecological Modelling — [ISEM]. 1162 pp. McCreary, J.P., Kohler, K.E., Hood, R.R., Olson, D.B., 1996. A four-component ecosystem model of biological activity in the Arabian Sea. Progress in Oceanography 37, 117–165. Nakata, K., 1989. Circulation of water of Tokyo Bay. Journal of Oceanographic Society of Japan 42, 5 (in Japanese). Neumann, T., 2000. Towards a 3-D ecosystem model of the Baltic Sea. Journal of Marine Systems 25, 405–419. Neumann, T., Fennel, W., Kremp, C., 2002. Experimental simulations with an ecosystem model of the Baltic Sea: a nutrient load reduction experiment. Global Biogeochemical Cycles 16, 7_1-7_19.

16

N.T. Minh Hang et al. / Journal of Marine Systems 75 (2009) 1–16

Pacanowski, R.C., Dixon, K. and Rosati, A. (1990). The GFDL Modular Ocean Model Users Guide Version 1.0, GFDL Technical Report No. 2, Geophysical Fluid Dynamic Laboratory, NOAA, Princeton University, Princeton. Parson, T.R., Takahashi, M., Hargrave, B., 1990. Biological Oceanographic Processes. Pergamon Press, Oxford. Rhee, G.-Y., Gotham, I.J., 1981. The effect of environmental factors on the phytoplankton growth: Light and the interactions of light with nitrate limitation. Limnology and Oceanography 26, 649–659. Sakshaug, E., Slagstad, D., 1991. Light and productivity of phytoplankton in polar marine ecosystems: a physiological view. Polar Research 10, 69–86. Sharples, J., Tett, P., 1994. Modelling the effect of physical variability on the midwater chlorophyll maximum. Journal of Marine Research 52, 219–238. Suzuki, T., Aoyama, H., Hata, K., 1997. The quantification of nitrogen cycle on tidal flat by ecosystem model—in the case of Isshiki tidal flat in Mikawa bay. Journal of Advanced Marine Science and Technology Society 3 (1), 63–80. Taylor, A.H., Harbour, D.S., Harris, R.P., 1993. Seasonal succession in the pelagic ecosystem of the North Atlantic and the utilization of nitrogen. Journal of Plankton Research 15, 875–891.

Umgiesser, G., Melaku Canu, D., Solidoro, C., Ambrose, R., 2003. A finite element ecological model: a first application to the Venice Lagoon. Environmental Modelling Software 18, 131–145. Yamanishi, H., Araki, H., Koga, K., Sato, K., 2002. Study on water quality and mud property in the inner part of Ariake Bay. Environmental Engineering Research 39, 219–227 (in Japanese, with English Abstr.). Yamanishi, H., Araki, H., Koga, Y., Himura, K., Ohishi, K., 2005. Field survey on suspended matters transport and variations of water quality on an intertidal mudflat in the gulf of Ariake Sea using an automatic profiler sensor. Environmental Engineering Research 42, 297–304 (in Japanese, with English Abstr.). Zhang, J., Jorgensen, S.K., Mahler, H., 2004. Examination of structurally dynamic eutrophication model. Ecological Modeling 173, 313–333. Zheng, L., Chen, C., Zhang, F.Y., 2004. Development of water quality model in the Satilla River Estuary. Ecological Modeling 178, 457–482.