Assesment of the response of a shallow macrotidal estuary to changes in hydrological and wastewater inputs through numerical modelling

Ecological Modelling 221 (2010) 1194–1208

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Assesment of the response of a shallow macrotidal estuary to changes in hydrological and wastewater inputs through numerical modelling Andrés García a,∗ , José A. Juanes a,1 , César Álvarez a,1 , José A. Revilla a,1 , Raúl Medina b,2 a b

Submarine Outfall and Environmental Hydraulics Group, Institute of Environmental Hydraulics, University of Cantabria, Avda. de los Castros s/n, 39005 Santander, Spain Ocean and Coastal Research Group, Institute of Environmental Hydraulics, University of Cantabria, Avda. de los Castros s/n, 39005 Santander, Spain

a r t i c l e

i n f o

Article history: Received 16 June 2009 Received in revised form 20 November 2009 Accepted 31 December 2009 Available online 4 February 2010 Keywords: Eutrophication Shallow estuary Modelling Phytoplankton Wastewater discharge Urdaibai

a b s t r a c t The aim of this study was to investigate the response to short-term changes in river freshwater discharges and in nutrients loadings (mainly from the treatment of urban wastewater), of the shallow macrotidal Urdaibai estuary (north of Spain), by using numerical tools. A two-dimensional hydrodynamic model and a water quality model were applied to the estuary, in order to better use it as a prediction tool in the study of the effects of variations in hydrodynamic conditions and in waste water inputs. The model was calibrated and verified using data measured under different hydrological conditions (spring and summer). A model calibration was carried out with field data measured during the summer, while the model validation was conducted for spring conditions. The calibration process allowed the model parameter definition, while the model validation permitted the verification of the calibrated parameters under different environmental conditions. The model results were in reasonable agreement with field measurements, in both the calibration and the validation phases. The model showed a significant decrease in phytoplankton concentration with river input increase. A study on the effects of nutrient input reduction from the Gernika Waste Water Treatment Plant (WWTP) was conducted. It showed a decrease in phytoplankton concentration with decreasing levels of nutrient discharge. This reduction was more pronounced in conjunction with the highest river discharge. In that case, a 50% decrease for the elimination of the WWTP discharge was observed. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Estuaries receive inorganic nutrients and organic matter from the land and represent important systems where terrestrial nutrients and organic matter are processed before entering the ocean (Ortega et al., 2005). As a result, the ecology and biodiversity of estuarine waters in many parts of the world are under the threat of increasing anthropogenic nutrient inputs (Nixon, 1995). Phytoplankton growth is limited by several environmental factors, including light, temperature and nutrients. Among these factors, only nutrients can be controlled; hence, they have been the focus of most efforts to control algal blooms responsible of water quality deterioration (Na and Park, 2006). Knowledge of the nutrient assimilation capacity of estuarine ecosystems is essential for water quality management and rehabilitation. The main adverse

∗ Corresponding author. Tel.: +34 942201704; fax: +34 942201714. E-mail addresses: [email protected] (A. García), [email protected] (J.A. Juanes), [email protected] (C. Álvarez), [email protected] (J.A. Revilla), [email protected] (R. Medina). 1 Tel.: +34 942201704. 2 Tel.: +34 942201810. 0304-3800/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2009.12.027

effects of uncontrolled eutrophication are depression of oxygen levels and the massive occurrence of harmful algal species. All this has negative effects on water quality and food webs. Thus, management of aquatic ecosystems has traditionally focused on reducing nutrient input. However, besides nutrient control strategies, hydrological conditions may have an impact on the response of estuarine ecosystems to eutrophication (Chan et al., 2002). Water mass circulation in estuaries is predominantly controlled by fresh water discharges, tidal circulation and atmospheric forces. The time scale of physical processes ranges from hours to seasons and large spatial gradients affect the hydrographic conditions and nutrient distribution (Dyer, 1997). The high number of factors influencing eutrophication processes has led to consider the use of complex mathematical tools for studying and predicting their evolution. Water quality models are essential tools to evaluate the impact of human activities in estuarine ecosystems. Shallow estuaries are characterized by the presence of wide tidal flats which yield a considerable variability in the water domain. This feature increases the difficulty of eutrophication modelling and restricts the application of many water quality models that have been successfully validated in environments in which water

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

domain changes have no substantial effects, such as deep estuaries and bays, coastal areas or enclosed seas. The objective of this study was to develop a quantitative understanding of the Urdaibai estuary response to hydrological inputs and wastewater loading. The Urdaibai estuary is a shallow unique habitat for many species. It receives urban inputs, especially from the village of Gernika, which cause its environmental degradation (Cortazar et al., 2008). The estuarine ecosystem is also highly variable. During spring and summer, intense phytoplankton blooms usually develop in the upper reaches of the estuary (Madariaga, 2002). Several previous studies have been carried out in order to characterize the short-term variability of phytoplankton in the Urdaibai estuary from field observations. However, scarce modelling was used due to the complex geometry of the estuary and the high variability of river flow and nutrient fluxes. For this study, we set up and implemented a two-dimensional water quality model for the estuary. Here, we present the calibration and validation of the model based on data measured at different meteorological and hydrodynamic conditions. Once calibrated, the model was used to study the effect of variations in hydrodynamic conditions on estuarine water quality. The model was also applied to analyse future scenarios of nutrient input reduction by the Gernika WWTP. The results and conclusions of these studies are described.

2. Study area The Urdaibai Estuary is a meso-macrotidal temperate estuary located in the Basque Country (43◦ 22 N, 2◦ 40 W), north of Spain (Fig. 1). An important geomorphological feature of the estuary is its shallowness (on average 2.6 m deep in the main channel). It has a drainage area of 149 km2 , being about 12.5 km long, from its mouth to the city of Gernika, and with a maximum width of 1.2 km. The estuary receives inputs from several rivers and it is wellmixed to partially mixed depending on river discharge and tidal range. The strongest salinity gradient occurs generally within less than 4 km of the upper estuary, where salinity oscillates between near zero and more than 30 psu. The estuary passes from a wellmixed state during periods of low river discharge and spring tides, to partially mixed during neap tides or freshets. Thus, the estuary is river dominated during freshets but tidally dominated during summer, when river discharges tend to be low (Ruiz et al., 1994). Water residence time is less than 1 day on average in the lower estuary due to the large tidal amplitude. Conversely, in the upper estuary it experiences drastic changes depending on river flow, oscillating from 1 day to more than 3 months during prolonged dry periods (Revilla et al., 2002). The Oka River basin has a drainage area of 67 km2 , with the main stream covering a distance of 25 km. The average rainfall is over 1400 mm per year and the overall runoff coefficient is 0.64. Other rivers entering the estuary are the Mape and the Golako. The drainage area of the latter is 34 km2 , while the basin drainage of the Mape River covers an area of 20 km2 . The estuary is surrounded by relatively broad tidal flats along its lower reaches and by salt-marshes in its middle part. These geomorphological features determine the estuarine ecosystem, which can vary considerably in volume and flushing rates, strongly affecting its biological and chemical properties (Madariaga, 2002). The catchment area is essentially rural, and industrial activities are mainly concentrated in the town of Gernika, which with a population of 16255 inhabitants, is located in the upper estuary. As a consequence, estuarine waters receive the discharge of a domestic sewage treatment plant effluent near Gernika that introduces

1195

a considerable amount of nutrients (mainly ammonia and phosphate) to the estuary, causing eutrophication (Fraile et al., 1992). 3. Model description The model reproduces estuarine water movements during a tidal cycle using short time steps from hydrodynamic model results (velocity components, water surface level) obtained by previously running a two-dimensional hydrodynamic model. This working scheme allows a larger time step to be selected in order to model eutrophication than that used for the velocity field calculation. But in the case of shallow estuaries, this methodology has the disadvantage that during hydrodynamic data interpolation inside the water quality model there is a loss of information about the real instant in which wetting and drying of the tidal flat occur. To avoid this problem, caused by using hydrodynamic results stored separately in time, a subroutine has been included in the model that allows water surface level to be computed in the entire estuarine domain from two consecutive water level recordings provided by the hydrodynamic model (García et al., 2002). In this way the exact instant in which shallow areas get wet or dry is accurately computed. 3.1. The hydrodynamic model The hydrodynamic and transport models used in this work solve the two-dimensional vertically integrated hydrodynamic and transport equations. The numerical computation is carried out on a spatial domain that represents the entire estuary through a finitedifference grid. The system of equations is expressed in Cartesian coordinates (x increasing eastward and y increasing northward) and following the hydrostatic assumption and Boussinesq approximation. The hydrodynamic model uses the so-called alternating direction implicit technique (ADI), to integrate the depth-averaged mass and momentum equations in the space–time domain, which are expressed as: ∂H ∂(UH) ∂(VH) + + =0 ∂t ∂x ∂y

(1)

∂(UH) ∂(U 2 H) ∂(UVH) ∂ gH 2 ∂0 + + = fVH − gH − 20 ∂x ∂t ∂x ∂y ∂x +

 1  xz() − xz(−h) + He 0

∂e +H ∂y



∂U ∂V + ∂y ∂x



∂2 U ∂2 U + 2 ∂x ∂y2



+ 2H

∂e ∂U ∂x ∂x



(2)

∂(VH) ∂(UVH) ∂(V 2 H) ∂ gH 2 ∂0 + + = −fUH − gH − 20 ∂y ∂t ∂x ∂y ∂y +

 1   − yz(−h) + He 0 yz()

∂e +H ∂x



∂U ∂V + ∂y ∂x



∂2 V ∂2 V + 2 ∂x ∂y2



+ 2H

∂e ∂V ∂y ∂y



(3)

where  is the water level (m), g is the gravitational acceleration (m s−2 ), H = h +  is the total water depth (m), h is the undisturbed water depth (m), U and V (m s−1 ) are the vertical integrated velocities, 0 is the average water density, e is the horizontal eddy viscosity,  iz() and  iz(−h) are the wind and bed stresses and f is the Coriolis parameter. The bottom shear stress is represented as a

1196

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Fig. 1. Study area and location of hydrodynamic and water quality sampling stations.

quadratic function of velocity parameterized in terms of the Chezy friction coefficient (C): Ui 1 =g  H iz(−h)



 HDx

∂T ∂x

 +

∂ ∂y

 HDy

∂T ∂y

 (7)

U2 + V 2

(4)

C2H

The Chezy friction term can be expressed as a variable dependent on depth, using the following formula: C = Mh1/6

(5)

where M is the inverse of the Manning friction coefficient. Salinity and temperature variations along the estuary are solved from the depth-averaged transport equation for these variables: ∂(UHS) ∂(VHS) ∂ ∂(HS) + + = ∂t ∂x ∂y ∂x

∂(UHT ) ∂(VHT ) ∂ ∂(HT ) + + = ∂t ∂x ∂y ∂x



∂S HDx ∂x



∂ + ∂y



∂S HDy ∂y



(6)

where S is the salinity (psu), T is the water temperature (◦ C) and Dx and Dy are the diffusivities in x and y. Density is calculated from spatially variable salinity and temperature through the expression suggested by UNESCO (1981). We chose a time interval of 5 s as time step for the hydrodynamic model in order to assure stable conditions. At open boundaries we prescribed the water level and freshwater flow, while at the open sea boundary, we prescribed tidal elevation from the Andersen–Grenoble version 95.1 model (Le Provost, 2001).

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

1197

The penetration of incoming solar radiation is described by the Lambert–Beer equation: Iz = I0 exp(−Ke z)

(11) (ly day−1 )

where Iz is the light intensity at depth z calculated from light surface intensity (I0 ) and the extinction coefficient (Ke ). Daily I0 , issued by the National Meteorological Institute from the meteorological station G057 located at Mungia, was used in the model. The total light attenuation coefficient is determined by the effect of water and chlorophyll a, and can be expressed for the Urdaibai estuary as a function of suspended organic matter and phytoplankton biomass as follows: Ke = 0.36948 + 0.99577

C

6

2.7

+ C4



(12)

The light-limitation function of Steele and Baird (Thomann and Mueller, 1987) is used in the model. Vertically averaged G(I) over a given water depth is integrated as: G(I) =

2.718 [exp(−˛1 ) − exp(−˛0 )] Ke H

(13)

where ˛1 =

I0 exp(−Ke H) Is

(14)

˛0 =

I0 Is

(15)

Fig. 2. Variables of the water quality model.

3.2. The water quality model The water quality model solves the system of differential equations that describe the main chemical and biological interactions. The variables of the model are the following (see Fig. 2): ammonia (C1 ), nitrate (C2 ), phosphate (C3 ), phytoplankton biomass (C4 ), dissolved BOD (C5 ), suspended BOD (C6 ), sediment BOD (C7 ) and dissolved oxygen (C8 ). The spatial and temporal evolution of these variables is influenced by external factors, such as incident solar radiation, or urban and freshwater discharges. The water quality model is coupled to the hydrodynamic model through the transport equation, which integrates the advection and the diffusion properties of the flow, as well as the basic processes occurring in the estuary water column: ∂(UHCi ) ∂(HCi ) ∂(VHCi ) ∂ + + = ∂t ∂x ∂y ∂x +

∂ ∂y

 HDy

∂Ci ∂y

 HDx

∂Ci ∂x



(8)

where Ci is the depth-averaged concentration for water quality (i), and Ri describes the chemical reaction terms, corresponding to the interaction equations for state variables. Algal growth is described by a first-order kinetic expression where the first-order growth rate is defined as the difference between the growth (Gp ) and the death (Dp ) rates. Phytoplankton growth (Gp ) is a function of light G(I), temperature G(T) and nutrients G(N), assuming the classical approach that these effects are multiplicative (Chapra, 1997): Gp = G(T )G(I)G(N)

(9)

The effect of temperature is expressed as (Thomann and Mueller, 1987): G(T ) = Gmax GT −20 max

G(N) = min



(C1 + C2 ) C3 , KmN + (C1 + C2 ) KmP + C3



(16)

The preferential uptake of ammonia as compared to nitrate for phytoplankton growth kinetics is described as: PNH3 = C1

(KmN

C2 KmN + + C1 )(KmN + C2 ) (C1 + C2 )(KmN + C2 )



(17)

Phytoplankton mortality is described as the sum of phytoplankton endogenous respiration and zooplankton grazing. The endogenous respiration rate is expressed as proportional to algal biomass and the carbon/chlorophyll (CCHL) ratio:

 + Ri H

being Is the saturating light intensity calculated as the weighted average of light intensity for the previous 3 days (Chau and Jin, 1998). All of the individual nutrient limitations are computed by a Michaelis–Menten-type expression and the minimum value is chosen for G(N) (Chao et al., 2007):

(10)

were Gmax is the maximum daily growth rate of phytoplankton at 20 ◦ C (day−1 ) under optimal light and nutrient conditions and Gmax is the temperature coefficient.

kr = 37.85 CCHL C4

(18)

The CCHL ratio (mg C/mg Chl-a) is considered a variable which is affected by light, temperature and algae (Cloern et al., 1995): CCHL = 0.003 + 0.0154exp(0.05 T ) exp(−0.059Iav )G(N)

(19)

where Iav is the depth-averaged daily irradiance. The effect of temperature over endogenous respiration is given by: kr (T ) = kr (20 ◦ C)rT −20

(20)

where kr is the endogenous respiration rate and  r is a temperature coefficient. In the water column, biological oxidation of NH3 to NO3 and denitrification are all considered a function of temperature and dissolved oxygen (Chau and Jin, 1998). The dissolved oxygen concentration at saturation is considered dependent on both temperature and salinity. To compute the oxygen saturation value we used the following equation (Lopes et al., 2008):





Cs = 14.652 − 0.0841S + T 0.0026S − 0.41022 + T ˇ(S, T )

(21)

1198

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Table 1 Water quality model interactions. Variable

Interaction equation

Ammonia-N

dC1 dt

Nitrate-N

dC2 dt

=

(T −20) kn n

Phosphate-P

dC3 dt

=

(T −20) Yd2 kd d C5

Phytoplankton-C

dC4 dt

= G(T )G(I)G(N) − (kr r

Dissolved BOD

dC5 dt

= −kd d

Suspended BOD

dC6 dt

=

(T −20)

= Yd kd d



(T −20)

C5 + Ys ks s C8 (C8 +knit )



(T −20)

Settled BOD

=

Dissolved oxygen

dC8 dt

= ka a

Vs H

C6 −

Vr H

(T −20)

(T −20)

Yb kb b

(T −20) Ys2 ks s C6

+

+ kdm ) −

Vph



H

+

(T −20) kdn dn

(T −20) apc (kr r



kNO (C8 +kNO )





C8 (C8 +knit )



C1

C2

+ kdm )(1 − fop )C4 − G(T )G(I)G(N)apc C4

C4

C7 −

+ kdm )(anc fon + apc fop )C4 −

(T −20) kb b C7

+

Vph H

Vs H

(T −20)

C6 − ks s

C4

(Cs − C8 ) + G(T )G(I)G(N)

C7 ) −

(T −20) Yb2 kb b C7

(T −20)

+ kdm )anc (1 − fon )C4 − G(T )G(I)G(N)PNH3 anc C4 − kn n

C5

C7 + (kr r

dC7 dt

(T −20)

C7 + (kr r

C1 − G(T )G(I)G(N)(1 − PNH3 )anc C4 −

(T −20)

(T −20)

Vr H

+



(T −20)

C6 + Yb kb b

 32 12

+

48 a (1 14 nc

− PNH3 ) C4 −

SOD (T −20) SOD H

with ˇ(S, T ) = 0.007991 − 0.0000374S − 0.000077774T



C6

(22)

The re-aeration kinetic constant ka is influenced by water temperature, flow characteristics (flow velocity and water depth) and meteorological conditions (wind). The developed model allows Eq. (8) to be solved by using an explicit finite-difference scheme based on the split operator approach, in which advection and diffusion processes are computed independently for each time step. This approach facilitates the use of different numerical methods for solving each process (Komatsu et al., 1997). Hence, advective transport is computed by an upwind scheme while diffusion is described by a centred scheme. In a further step, changes in water quality concentration caused by reaction processes are computed. The time step for the transport model was 30 s, in order to assure stable solutions for transport equations. The whole system of equations which describe the interactions of 8 state variables is presented in Table 1. The symbols along with the units used for the model are listed in Table 2. The discharge of the Gernika WWTP is treated as a point-source input of organic matter and nutrients into the water quality model. At the discharge cell, a mass balance assuming complete mixing is made, due to its small water depth. 4. Model calibration 4.1. Field data A model calibration may be defined as an operation by which specific values, distributions, or a range of variations are given to the floating free model parameters, so that the model results fit optimally to a set of field observations (Lopes et al., 2008). To calibrate the model, both hydrodynamic and water quality data were needed. The hydrodynamic model was calibrated and validated with data measured between January 1998 and February 1998 at six stations distributed along the estuary. The collected data included surface elevation, velocities and freshwater inputs from the Oka and Mape rivers. The location of the sampling stations is shown in Fig. 1. Two tide gauge stations were located at the middle part of the estuary (Murueta shipyards) and at the inner part (near the Gernika WWTP), respectively, during a time period of 15 days. A tidal distortion is observed as the tide moves landwards. This distortion is more pronounced during spring tides at low tide. Velocity measurements were carried out at four different locations (stations H1, H2, H5 and H7) distributed along the estuary during neap and spring tides. Maximum velocities were registered

(T −20) 32 k  C4 12 r r

−

(T −20) 64 k  14 n n



C8 (C8 +kn )



(T −20)

C1 − (Yd kd d

(T −20)

C5 + Ys ks s

C6 +

during spring tides near the mouth of the estuary with values over 1 m s−1 . In the upper estuary maximum velocities were about 0.8 m s−1 . The velocity field during spring tides was more than twice the magnitude found during neap tides. Time series of water quality variables were measured for a 4-day period in May and July 1999, at four locations in the estuary (stations WQ1, WQ3, WQ5 and WQ7), which are indicated in Fig. 1. Vertical profiles of salinity, temperature, dissolved oxygen and chlorophyll a were taken. The sampling carried out in May was intended to evaluate the estuary recovery after a period of heavy rain. That of July, in contrast, was conducted during a relatively dry period in order to see the daily changes in water quality associated with changes in radiation. In May, salinity remained low and constant at the uppermost site during the study period. In contrast, station WQ5 experienced strong changes in water salinity, which varied in the surface from 1.1 psu, on the first day, to 12.5 psu on the last one. Temperature experienced drastic changes along the estuary, mainly at station WQ5 which ranged from 16 to 21 ◦ C. Oxygen concentrations remained near saturation level during the sampling period. Temperature experienced a significant increase at the inner stations, while it remained constant in the sea area. Chlorophyll a values were low, both in the river and at the upper estuary. The average Oka River discharge at the Muxika gauging station in this sampling period was less than 0.5 m3 s−1 , and the tidal range varied between 2.0 and 2.7 m. During the summer the main physical feature was the decrease in the oxygen concentration from the mouth towards the upper estuary. Salinity and temperature remained rather constant. In the upper estuary, the water column stratification caused by low river discharges enhanced phytoplankton growth and plankton metabolism, leading to anoxic conditions for most of the water column. The concentration of chlorophyll a was typical of the summer, with peaks at station WQ8, and the minimum values in the sea area. At station WQ1 values lower than 2 ␮g l−1 were observed while at station WQ8, the daily changes in the concentration of chlorophyll, which ranged between 14.8 and 50.2 ␮g l−1 , were due to the effect of solar radiation. The average Oka River discharge during this period was about 0.13 m3 s−1 and the tidal range varied between 1.7 and 2.0 m. The spatio-temporal variability in the estuary hydrographic properties during the two surveys is summarized in Table 3. Table 4 shows the nutrient concentrations measured at the mouth of the main rivers flowing into the estuary. Benthic fluxes of dissolved oxygen were derived from data collected in five sampling stations along the estuary by Ortega et al. (2005). The sediment oxygen demand (SOD) varied from 1.5 gO2 m−2 day−1 at the mouth of the estuary to 8.4 gO2 m−2 day−1

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Fig. 3. Sea surface elevation from model results (continuous line) and from the tide gauge located at the Murueta shipyards and at the Gernika WWTP.

Fig. 4. Depth-averaged velocities measured (black dots) and modelled (continuous line) during neap tides.

1199

1200

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Table 2 Coefficients and constants used in the model. Parameter

Description

Value

Units

Fr d kd ks kb d s b Vr Vs Yd Ys Yb Yd2 Ys2 Yb2 anc apc kn n kdn  dn knit kNO fon fop  SOD Gmax Gmax kmN kmP kr r kdm Vph a

Fraction of dissolved BOD BOD dissolved degradation rate BOD suspended degradation rate BOD settled degradation rate Temperature coefficient for dissolved BOD degradation Temperature coefficient for suspended BOD degradation Temperature coefficient for settled BOD degradation Resuspension rate of particulate BOD Settling rate of particulate BOD Yield factor for release of ammonia by ammonification of dissolved BOD Yield factor for release of ammonia by ammonification of suspended BOD Yield factor for release of ammonia by ammonification of settled BOD Yield factor for release of phosphate from dissolved BOD Yield factor for release of phosphate from suspended BOD Yield factor for release of phosphate from settled BOD Nitrogen/carbon ratio Phosphorus/carbon ratio Nitrification rate Temperature coefficient for nitrification Denitrification rate Temperature coefficient for denitrification Half saturation constant of oxygen for nitrification Half saturation constant for oxygen limitation denitrification Fraction of dead algae recycled to organic nitrogen Fraction of dead algae recycled to organic phosphorus Temperature coefficient for SOD Maximum specific growth rate of phytoplankton Temperature coefficient for growth rate of phytoplankton Half saturation constant for uptake of inorganic nitrogen Half saturation constant for uptake of inorganic phosphorus Endogenous respiration rate of phytoplankton Temperature coefficient for endogenous respiration Maximum grazing rate by zooplankton Settling speed of phytoplankton Temperature coefficient for re-aeration

0.5 0.4 0.4 0.2 1.047 1.047 1.047 0.1 0.2 0.1 0.1 0.1 0.01 0.01 0.01 0.15 0.009 0.05 1.088 0.1 1.05 2.0 0.1 0.5 0.5 1.09 2.0 1.067 0.025 0.001 0.1 1.045 0.1 0.1 1.024

– day−1 day−1 day−1 – – – m day−1 m day−1 g NH3 -N/g BOD g NH3 -N/g BOD g NH3 -N/g BOD g PO4 -P/g BOD g PO4 -P/g BOD g PO4 -P/g BOD mg N/mg C mg P/mg C day−1 – day−1 – gO2 m−3 gO2 m−3 – – – day−1 – mg N l−1 mg N l−1 day−1 – day−1 m day−1 –

Fig. 5. Depth-averaged velocities measured (black dots) and modelled (continuous line) during spring tides.

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Fig. 6. Time series of depth-averaged salinity (䊉, measured; —, modelled) and temperature (, measured; - - -, modelled) during spring conditions.

Fig. 7. Time series of depth-averaged salinity (䊉, measured; —, modelled) and temperature (, measured; - - -, modelled) during summer conditions.

1201

1202

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Fig. 8. Depth-averaged dissolved oxygen (䊉, measured; — modelled) and phytoplankton-C (, measured; - - -, modelled) during summer conditions.

Fig. 9. Depth-averaged ammonia (䊉, measured; —, modelled) and phosphate (, measured; - - -, modelled) during summer conditions.

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

1203

Fig. 10. Depth-averaged dissolved oxygen (䊉, measured; —, modelled) and phytoplankton-C (, measured; - - -, modelled) during spring conditions.

at a location near the Gernika WWTP. Different SODs were considered for the whole estuary taking into account the bed characteristics and the distance to measuring points. The discharge flow of the Gernika WWTP ranged between 1560 and 6990 m3 day−1 with an average value of 3970 m3 day−1 , a mean 5-day BOD of 60.6 mgO2 l−1 , and an average concentration of ammonium higher than 52 mg NH4 + l−1 .

ment with measured data but small deviations are observed during spring tides, mainly at the shipyards of Murueta. However, this discrepancy which was due to the complexity inherent in the geometric representation of the estuary and the processes of tidal wave deformation occurring along it, was not significant. Modelled time series of velocities occurring along the estuary for four different locations were compared with measurements taken during neap tide (Fig. 4) and spring tide (Fig. 5). The computed velocities were generally quite close to the measured values. The largest differences were observed at stations 1 and 2, which were located at the upper site of the estuary. Therefore, it can be concluded that the hydraulic model reproduced water movement satisfactorily throughout the estuary. The calibration procedure for the dispersion coefficient was performed against the observed salinity for the four stations of the estuary. Figs. 6 and 7 show the comparison between modelled and observed salinity and temperature. Model salinity results compared well with field measurements. A constant dispersion coefficient of 0.5 m2 s−1 was obtained.

4.2. Hydrodynamic model calibration The model calibration was performed adjusting the bottom friction coefficient for the entire estuary. In this study the best adjustment between model results and field observations was achieved through bottom roughness parameterized from a variable Chezy coefficient dependent on water depth. A friction coefficient M of 45 and an eddy viscosity of 1 m2 s−1 were obtained. Fig. 3 shows the comparison between computed and observed water surface elevation time series for the two stations located in the estuary. The model results are in general in good agree-

Table 3 Mean values (±STD) of hydrodynamic and water quality parameters at each sampling station. Survey

Station

Variable Salinity (psu)

Temperature (◦ C)

Dissolved oxygen (mg l−1 )

Chl-a (␮g l−1 )

PO4 3− (mg l−1 )

NH4 + (mg l−1 )

Spring

WQ1 WQ3 WQ5 WQ8

34.2 25.3 11.2 0.1

± ± ± ±

0.7 3.6 5.0 0.0

17.3 18.8 19.0 18.6

± ± ± ±

0.5 1.0 1.8 1.6

8.1 7.4 7.0 9.1

± ± ± ±

0.3 0.3 0.5 1.4

1.60 2.44 3.10 0.90

± ± ± ±

0.77 1.10 2.61 0.30

0.31 0.32 0.24 0.37

± ± ± ±

0.15 0.11 0.06 0.32

0.06 0.14 0.46 1.36

± ± ± ±

0.02 0.09 0.14 0.79

Summer

WQ1 WQ3 WQ5 WQ8

34.7 31.4 24.5 14.9

± ± ± ±

0.2 0.6 1.8 3.3

22.0 22.8 23.3 23.1

± ± ± ±

0.4 0.4 0.5 0.4

7.2 5.5 3.4 2.6

± ± ± ±

0.3 0.2 0.6 2.0

1.13 2.57 2.98 29.93

± ± ± ±

0.33 1.17 0.56 17.50

0.43 0.60 1.20 2.66

± ± ± ±

0.22 0.12 0.08 0.50

0.05 0.30 1.50 4.92

± ± ± ±

0.01 0.09 0.19 2.17

1204

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Table 4 Mean values (±STD) of water quality parameters in the main rivers. River

Survey

NO3 − (mg l−1 )

PO4 3− (mg l−1 )

NH4 + (mg l−1 )

Chl-a (␮g l−1 )

Oka

Spring Summer

4.00 ± 0.61 4.86 ± 0.52

– 0.16 ± 0.02

– 0.49 ± 0.36

0.84 ± 0.17 0.23 ± 0.10

Golako

Spring Summer

3.26 ± 0.55 2.16 ± 0.25

– 0.12 ± 0.05

– 0.57 ± 0.24

– –

Mape

Spring Summer

1.57 ± 0.65 1.64 ± 0.37

– 0.12 ± 0.05

– 0.50 ± 0.18

– –

4.3. Water quality model calibration and validation The calibration procedure for the water quality model was performed against data measured along the 4-day survey in July 1999. Water quality comparisons between field data and model results are presented for dissolved oxygen and phytoplankton in Fig. 8. Fig. 9 shows the comparison of nutrient concentrations measured and modelled at stations WQ1 to WQ8. The modelled curves show an important amplitude variation with tidal oscillations that could not be measured during the field survey since only one daily sample was taken at each station. This type of oscillations has been found in other mesotidal estuaries, such as the Ria de Aveiro lagoon (Lopes et al., 2008). These fluctuations of the amplitude of the concentration made it more difficult to compare model to data. As is shown in Figs. 8 and 9, the model reproduces the general trend of data. Concerning dissolved oxygen, the mean relative error is lower than 10% in stations WQ1 to WQ5. The model overestimates phytoplanktonC concentrations for stations WQ3 and WQ5 and underestimates the phytoplankton-C concentrations for stations WQ1 and WQ8. The mean relative error is of the order of 20%. In the case of nutri-

ents, the model reproduces the order of magnitude of data, even though the estimated relative error is of the order or greater than 20%. The calibrated parameter values are presented in Table 2. Thus, the model could be used as a predictive tool for computing algal growth dynamics and other water quality processes affecting dissolved oxygen concentrations in shallow water bodies. Nevertheless, we developed a validation process using the parameters obtained considering the data acquired during May 1999. This process focused on verifying the model applicability to different environmental conditions. The comparison between computed and observed dissolved oxygen and phytoplankton levels is presented in Fig. 10. Fig. 11 shows the same comparison but for nutrient levels. In these figures it can be observed that the model reproduces the general trend of observed values and the order of magnitude of data. The mean relative error is lower than 10% for dissolved oxygen and phosphate concentrations. Concerning phytoplankton-C and ammonia concentrations, the mean relative error is of the order or greater than 20%. The model underestimates the phytoplankton-C concentration for station WQ1 and overestimates the phytoplankton levels in the innermost stations.

Fig. 11. Depth-averaged ammonia (䊉, measured; —, modelled) and phosphate (, measured; - - -, modelled) during spring conditions.

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Fig. 12. Model results at the Gernika WWTP discharge.

Fig. 13. Time series showing the effect of the hydrodynamic input variability on depth-averaged phytoplankton levels in spring time.

1205

1206

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

Fig. 14. Time series showing the effect of the hydrodynamic input variability on depth-averaged phytoplankton levels in summer time.

The model results for both periods at the Gernika WWTP discharge are presented in Fig. 12. High nutrient levels as well as low dissolved oxygen concentrations in summer conditions were obtained as a consequence of urban waste water inputs.

5. Modelling of the estuary response to environmental changes Previous studies carried out in the Urdaibai estuarine ecosystem (Madariaga and Orive, 1989) showed that phytoplankton biomass and production in the upper estuary was low during seasonal periods of high river flow and short residence time, but highly increased

during the summer season, when lower river flow and higher residence time occurred. As has been mentioned before, two field surveys were carried out during the year 1999 considering different environmental conditions. The first one was conducted at the end of May, trying to represent spring conditions after an increase in freshwater discharges. The second one was carried out at the end of July, under summer conditions. The comparison between the two periods shows that the inner zone of the estuary had a higher concentration of oxygen when salinity was very low, whereas when salinity increased the oxygen concentration decreased. In the same way, station 5, also in the Gernika channel, only went into hypoxic conditions when salinity increased. The measured data reflected

Fig. 15. Comparison of mean phytoplankton levels at stations WQ3, WQ5 and WQ8 for several pollutant loads of the Gernika WWTP discharge.

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

1207

Fig. 16. Comparison between the phytoplankton biomass concentration levels computed in the current situation and in the hypothetical scenario of eliminating the discharge of the Gernika WWTP during spring conditions (—, without WWTP discharge; - - -, current situation).

Fig. 17. Comparison between the phytoplankton biomass concentration levels computed in the current situation and in the hypothetical scenario of eliminating the discharge of the Gernika WWTP during summer conditions (—, without WWTP discharge; - - -, current situation).

1208

A. García et al. / Ecological Modelling 221 (2010) 1194–1208

a significant decrease in the phytoplankton concentration as the input of fresh water from the river Oka increased, reaching the minimum level after a maximum peak in river discharge. To analyse the influence of the hydrodynamic input variability on changes in the phytoplankton level, water quality simulations were extended to a 15-day period. In this way the effect of tidal variations on phytoplankton dynamics was taken into account. Figs. 13 and 14 show the phytoplankton concentrations obtained at stations WQ1, WQ5 and WQ8 for spring and summer conditions, respectively. These results highlight the influence that stems from the periodic oscillation of the astronomical tide, i.e., the concentration levels of phytoplankton are higher during periods coinciding with neap tides, with a gradual reduction as the tide rises to periods of spring tides. The relationship between oxygen concentration and salinity reflects the impact of water residence time on productivity and water quality in the estuary. Productivity can be estimated from phytoplankton biomass, which in previous studies has been shown to be highly correlated with phytoplankton and bacterial production and respiration rates in the community. On the other hand, oxygen concentration is a good indicator of water quality. The inverse relationship between productivity and water quality in the inner estuary makes it necessary to apply mathematical models to assist in making decisions about river flow management and reduction of waste water loads. Several simulations were run considering different lower pollutant load levels for the Gernika WWTP discharge with respect to the current situation. Also, a scenario of total suppression of the waste water discharge was analyzed. In such situations, water quality indicator variables, especially nutrients and phytoplankton biomass, were simulated. Fig. 15 presents the model results for the different input levels studied. The curves shown in this figure indicate the mean 15-day phytoplankton-C concentrations obtained at stations WQ3, WQ5 and WQ8 as a function of the WWTP discharge level. It can be observed that the reduction in the nutrient input caused a lower concentration of phytoplankton within the system. Obviously, the lowest phytoplankton levels were associated to the zero discharge. This phenomenon was more pronounced for spring conditions, that corresponded to higher river discharges (see also Fig. 16), with a decrease of over 50% in the phytoplankton concentration in stations WQ3 and WQ8. In July, this reduction was smaller (Fig. 17), due to nutrient concentration decrease was partly neutralized by the improvement in the environmental conditions for the development of phytoplankton, with respect to those of May. A phytoplankton concentration decrease of around 20% was obtained for water quality stations WQ3 and WQ8. 6. Conclusions A two-dimensional mathematical model for computing eutrophication changes in the shallow macrotidal Urdaibai estuary (Basque Country), considering the water domain variability related to wetting and drying phenomena of the shallowest areas, was developed and applied. The model was calibrated and verified both with hydrodynamic and water quality data measured in several field surveys. Water quality data collected in May characterized estuarine conditions after a period of strong rainfall, while data measured in July reflected dry weather conditions.

The reduction in the magnitude of the currents within the estuary due to a decrease in river discharges and tidal range led to an increase in available nutrients and therefore a significant increase in phytoplankton biomass, especially in the innermost stations. This could increase their potential risk of eutrophication. The results of modelling considering several reduced input levels from the Gernika WWTP, concluded that there was a decrease in the phytoplankton concentration in the estuary. The maximum decrease was obtained with the “elimination of the WWTP discharge” scenario. This reduction was more pronounced in spring, when the weather and the increased river discharge were more favourable to this decline. In such conditions a decrease of about 50% in phytoplankton concentrations was predicted. This phenomenon was not so clear under summer conditions. References Chan, T.U., Hamilton, D.P., Robson, B.J., Hodges, B.R., Dallimore, C., 2002. Impacts of hydrological changes on phytoplankton succession in the Swan River, Western Australia. Estuaries 25, 1406–1415. Chao, X., Jia, Y., Shields Jr., F.D., Wang, S.Y., Cooper, C.M., 2007. Numerical modelling of water quality and sediment related processes. Ecol. Model. 201, 385– 397. Chapra, S.C., 1997. Surface Water-Quality Modeling. McGraw-Hill, Singapore, 844 pp. Chau, K.W., Jin, H., 1998. Eutrophication model for a coastal bay in Hong Kong. J. Environ. Eng. 124, 628–638. Cloern, J.E., Grenz, C., Vidergar-Lucas, L., 1995. An empirical model of the phytoplankton chlorophyll:carbon ratio—the conversion factor between productivity and growth rate. Limnol. Oceanogr. 40, 1313–1321. Cortazar, E., Bartolomé, L., Arrasate, S., Usobiaga, A., Raposo, J.C., Zuloaga, O., Etxebarria, N., 2008. Distribution and bioaccumulation of PAHs in the UNESCO protected natural reserve of Urdaibai, Bay of Biscay. Chemosphere 72, 1467– 1474. Dyer, K.R., 1997. Estuaries: A Physical Introduction. Wiley, Chichester, 195 pp. Fraile, H., Franco, J., Ruiz, A., Orive, E., 1992. Algunos datos sobre variables indicadoras del estado trófico del estuario de la Reserva de Urdaibai. Kobie, pp. 33– 37. Le Provost, C., 2001. Ocean tides. In: Fu, L., Cazenave, A. (Eds.), Satellite Altimetry and Earth Sciences: A Handbook of Techniques and Applications. Academic Press, 463 pp. García, A., Revilla, J.A., Medina, R., Álvarez, C., Juanes, J.A., 2002. A model for predicting the temporal evolution of dissolved oxygen concentration in shallow estuaries. Hydrobiologia 475–476, 205–211. Komatsu, T., Ohgushi, K., Asai, K., 1997. Refined numerical scheme for advective transport in diffusion simulation. J. Hydraul. Eng. 123 (1), 41–50. Lopes, J.F., Silva, C.I., Cardoso, A.C., 2008. Validation of a water quality model for the Ria de Aveiro lagoon, Portugal. Environ. Model. Softw. 23, 479–494. Madariaga, I., Orive, E., 1989. Spatio-temporal variation of size-fractioned primary production in the Gernika estuary. J. Exp. Mar. Biol. Ecol. 127, 273–288. Madariaga, I., 2002. Short-term variations in the physiological state of phytoplankton in a shallow temperate estuary. Hydrobiologia 475–476, 345–358. Na, E.H., Park, S.S., 2006. A hydrodynamic and water quality modelling study of spatial and temporal patterns of phytoplankton growth in a stratified lake with buoyant incoming flow. Ecol. Model. 199, 298–314. Nixon, S.W., 1995. Coastal marine eutrophication: a definition, social causes, and future concerns. Ophelia 41, 119–129. Ortega, T., Ponce, R., Forja, J., Gómez-Parra, A., 2005. Fluxes of disolved inorganic carbon in three estuarine Systems of the Cantabrian Sea (north of Spain). J. Mar. Syst. 53, 125–142. Revilla, M., Ansotegui, A., Iriarte, S., Madariaga, I., Orive, E., Sarobe, A., Trigueros, J.M., 2002. Microplankton metabolism along a trophic gradient in a shallowtemperate estuary. Estuaries 25, 6–18. Ruiz, A., Franco, J., Orive, E., 1994. Suspended particulate matter dynamics in the shallow mesotidal Urdaibai Estuary (Bay of Biscay, Spain). Neth. J. Aquat. Ecol. 28, 309–316. Thomann, R.V., Mueller, J.A., 1987. Principles of Surface Water Quality Modelling and Control. Harper Collins Publishers, New York, 644 pp. UNESCO, 1981. Background papers and supporting data on the International Equation of State of Seawater. Unesco Technical Papers in Marine Science 38, France, 193 pp.