One- and three-dimensional biogeochemical simulations of two differing reservoirs

Ecological Modelling 174 (2004) 143–160

One- and three-dimensional biogeochemical simulations of two differing reservoirs J.R. Romero∗ , J.P. Antenucci, J. Imberger Centre for Water Research, University of Western Australia, Crawley, WA 6009, Australia

Abstract A one-dimensional hydrodynamic model was coupled to an aquatic ecological model and applied to two differing reservoirs with one set of parameters. Simulations over 2 years of small (area = 5 km2 ) and shallow (average depth = 9 m) Prospect Reservoir (1988–1990), and large (volume = 2 km3 , area = 82 km2 , length = 50 km) and deep (maximum depth = 90 m) Lake Burragorang (1992–1994) were validated with field data. The dominant fate processes of the nitrogen and phosphorus cycles in the water column of these reservoirs were identified with model output. CAEDYM was then coupled to a three-dimensional hydrodynamic model (ELCOM) to simulate a flood underflow through Lake Burragorang. Advective transport was the dominant mechanism that caused biogeochemical variations during the flood, though settling was important for riverine-derived particulates. Validation of ecological models with both one- and three-dimensional hydrodynamics drivers and across multiple aquatic settings is advocated as an additional approach to limit ranges of biogeochemical parameters and assess ecological representations of lacustrine systems. © 2004 Elsevier B.V. All rights reserved. Keywords: Reservoirs; Biogeochemical modeling; Nitrogen; Phosphorus

1. Introduction Application of numerical water quality models to evaluate aquatic management strategies is widely used for lakes, reservoirs, rivers, estuaries, and coastal zones. A range of modeling approaches is available to evaluate aquatic ecosystems as summarized in Jørgensen and Bendoricchio (2001). Here, we consider the ‘reductionist/analytical’ approach with a general plankton model. In contrast to other modeling approaches that capture the dominant processes with a few parameters, the aim of these general plankton models is to model the major biogeochemical processes. For such general plankton models, a goal is to ∗ Corresponding author. Tel.: +61-8-9380-1685; fax: +61-8-9380-7115. E-mail address: [email protected] (J.R. Romero).

establish parameter values that are independent of the aquatic setting, and to identify those parameters that are site-specific. Such approaches have been considered in the past in lakes and reservoirs with one-dimensional (1D) frameworks such as CE-QUAL-R1 (USCE, 1995), MINILAKE (Riley and Stefan, 1988), and DYRESMWQ (Hamilton and Schladow, 1997), and with multidimensional models such as WASP (Ambrose et al., 1993). These 1D water quality models (CE-QUALR1, MINILAKE, DYRESM-WQ) were generally integrated into the architecture of an existing vertical mixing model. In the case of the multi-dimensional WASP modeling system, output from a prior hydrodynamic simulation is needed to provide transport for biogeochemical model runs. Desktop computational power is now sufficient to evaluate general plankton models in a ‘multi-system’ (several lakes and/or

0304-3800/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2004.01.005

144

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

reservoirs) validation framework to bound parameter values. Further, coupling of biogeochemical models with both 1D and 3D hydrodynamic drivers allows for evaluation of time and space scale dependencies of parameters. In this study a general plankton model was coupled to a 1D vertical mixing model and applied to two differing reservoirs with a single parameter set. This 1D configuration was used to evaluate the seasonal dynamics and to quantify the primary fluxes of nitrogen and phosphorus in two differing reservoirs. The general plankton model was also coupled to a three-dimensional (3D) hydrodynamic model to evaluate spatial and temporal biogeochemical distributions during a flood in one of the reservoirs. The aim of this study was to evaluate whether a simple ecological configuration and one set of biogeochemical parameters could be coupled to both 1D and 3D hydrodynamic drivers to adequately simulate the nutrient cycles in two differing reservoirs.

2. Study sites Both reservoirs are located in (Prospect Reservoir) or near (Lake Burragorang) Sydney, Australia (Fig. 1A). Prospect Reservoir is a moderately sized (5.25 km2 ) reservoir with two primary inflows (Pipeline and Upper Canal) along the western margin and four outlets around the perimeter (Fig. 1B). Prior to 1996 both of the inflows discharged directly into Prospect Reservoir. Selective withdrawal from Lake Burragorang is the source of the Pipeline waters. Those of the Upper Canal derive from reservoirs to the south of Sydney. The average residence time was low at ca. 1 month over the 2 years of this study (February 1988 to February 1990). Past modeling studies of Prospect Reservoir have been conducted to predict the effects of diversion of inflows from the reservoir and the cessation of withdrawals (Hamilton, 1999; Hamilton et al., 1995; Schladow and Hamilton, 1995). Lake Burragorang is a large reservoir (volume of 2 km3 , length of 60 km) (Fig. 1C) with a water balance that is dominated by occasional floods and consistent withdrawals for Sydney supply. The seven primary inflows into Lake Burragorang are gauged with most of the river confluences further than 30 km from the dam. The long, narrow and deep morphology of the reser-

voir results in the development of horizontal gradients of dissolved oxygen because of the greater sediment area to hypolimnetic volume ratio in the upper reaches (Romero and Imberger, 1999).

3. Description of models The Dynamic Reservoir Simulation Model (DYRESM) is a one-dimensional hydrodynamic model that simulates the temperature, salinity, and density in lakes and reservoirs. The model is based on a Lagrangian architecture that models the lake as horizontal layers of uniform properties (i.e. temperature and salinity) (Fig. 2). The layers move vertically, expanding and contracting in response to mass fluxes and mixing processes. The model explicitly simulates fundamental mixing mechanisms in stratified lakes and reservoirs (Imberger and Patterson, 1981, 1990). A new version of DYRESM has recently been completed with updates to the mixing parameterisations and a upgrade to the model architecture (Gal et al., 2003) to allow coupling to an aquatic ecological model (see the subsequent text). The updated version of DYRESM has been successfully applied to simulate the seasonal thermal evolution of lake Kinneret over 3 years (Gal et al., 2003). The Estuary and Lake Computer Model (ELCOM) is a 3D hydrodynamic model (Fig. 2) that has been applied to investigate internal waves in lakes (Hodges et al., 2000) and floods in reservoirs (Romero and Imberger, 2003). ELCOM is based on the unsteady Reynolds-averaged, hydrostatic, Boussinesq, Navier–Stokes, and scalar transport equations with an eddy-viscosity approximation for the horizontal turbulence correlations (Hodges et al., 2000). The model solves the unsteady Reynolds-averaged equations using a semi-implicit method with quadratic Euler–Lagrange discretization of momentum advection (Cassulli and Cheng, 1992) and a conservative ULTIMATE QUICKEST approach for scalar transport (Leonard, 1991). A one-dimensional mixed-layer model (Imberger and Patterson, 1990; Spigel et al., 1986) is included for turbulence closure of vertical Reynolds stress terms, and thus the turbulent fluxes (Hodges et al., 2000; Laval et al., 2003). Either of these two hydrodynamic models (DYRESM, ELCOM) can be readily coupled to the Com-

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

145

Fig. 1. (A) Location of Sydney region in Australia. (B) Bathymetry and locations of inflows, outflows, the meteorological station, and the sampling station in Prospect Reservoir. (C) Shoreline, river confluences, and sampling stations in Lake Burragorang. (D) Straightened bathymetry of Lake Burragorang for 3D simulations and corresponding sampling station locations.

146

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

Fig. 2. Schematic representation of the modeled nitrogen, phosphorus, and dissolved oxygen processes in CAEDYM. The (∗ ) symbol represents dissolved oxygen consumption from nitrification. Fluxes between the water column, sediment, and atmosphere are depicted as well as nutrient and dissolved oxygen fluxes from inflows and outflows. Instantaneous atmospheric losses of nitrogen from denitrification in the water column were modeled. The difference in discretization between the 1D and 3D models is highlighted in the lower portion of the diagram.

putational Aquatic Ecosystems Dynamic Model (CAEDYM). CAEDYM contains process descriptions for primary production, secondary production, nutrient cycling, and oxygen dynamics (Griffin et al., 2001; Romero and Imberger, 2003). State variables (Table 1) and ecological parameterizations (Table 2) of CAEDYM are similar to those of DYRESM-WQ (Hamilton and Schladow, 1997). CAEDYM was configured to simulate the dynamics of phosphorus, nitrogen, dissolved oxygen, and algae (Fig. 2) with one set

of parameter values (Table 3) for the three cases of this study. CAEDYM has been applied to the winter flood in Lake Burragorang of this study previously (Romero and Imberger, 2003), and zooplankton–phytoplankton interactions in the Swan River Estuary of Western Australia (Griffin et al., 2001). During a simulation, ecological processes (Fig. 2, Table 2) are updated by CAEDYM after each ELCOM time step (i.e. transport, mixing, meteorology, inflows, and outflows). Filterable reactive phospho-

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

147

Table 1 State variables simulated by CAEDYM Variable

Common name

Process description

DO

Dissolved oxygen

FRP POP

Filterable reactive phosphorus Particulate phosphorus = TP − IP − FRP where TP is total P and IP is algal internal P Ammonium Nitrate Particulate nitrogen = TN − IN − NO3 − NH4 where TN is total N and IN is algal internal N Diatoms Greens Blue-greens Areal mass of diatoms on the sediments Areal mass of green algae on the sediments Areal mass of blue-green algae on the sediments

Algal production/respiration, organic decomposition, nitrification, surface exchange, sediment oxygen demand Algal uptake, organic mineralization, sediment flux Mineralization, settling, algal mortality/excretion

NH4 + NO3 − PON AD AG AB ADmass AGmass ABmass

Algal uptake, nitrification, organic mineralization, sediment flux Algal uptake, nitrification, denitrification, sediment flux Mineralization, settling, algal mortality/excretion Growth, respiration, settling, resuspension Growth, respiration, settling, resuspension Growth, respiration, settling, resuspension Settling, resuspension Settling, resuspension Settling, resuspension

Water temperatures (T) and currents modeled by the hydrodynamic models DYRESM and ELCOM. In the table above a list of modeled mechanisms for each state variable is provided. All state variables also have inflow and withdrawal fluxes. The units of the state variables are: DO, mg DO l−1 ; FRP/POP, mg P l−1 ; NH4 + /NO3 − /PON, mg N l−1 ; Ai , mg chla l−1 ; Aimass , g chla m−2 .

rus (FRP) processes include algal uptake, sediment fluxes, and organic matter mineralization. Phytoplankton uptake of dissolved inorganic nitrogen was analogous to FRP except for a preference factor (PN ) for NH4 + over NO3 − , and the additional processes of nitrification and denitrification. Sources of particulate organic nitrogen (PON) and phosphorus (POP) included inflows, algal mortality and excretion, and settling from higher layers, while mineralization and settling to lower layers were losses. Because shear stresses have not been validated with DYRESM, no resuspension of detritus was modeled. However, algal resuspension was simulated through application of a small critical shear stress (Table 3) to ensure that phytoplankton at the sediment–water interface in shallow regions remained in suspension. The bed shear stress was computed with currents from the hydrodynamic drivers, and CAEDYM estimates of the orbital velocities from surface waves (Romero et al., 2003). The oscillatory currents generated by surface waves only influence the upper several meters of the water column. Hence, any phytoplankton that settled onto the sediments of the upper several meters of the water column were simulated to resuspend because of surface waves during moderate winds. Simulated dissolved oxygen (DO) processes included exchange across the air–water

interface, photosynthetic production, phytoplankton respiration, organic matter decomposition, and nitrification. CAEDYM was configured with three phytoplankton groups (diatom, chlorophyte, cyanophyte) representative of Prospect Reservoir and Lake Burragorang (Fig. 2, Table 1). CAEDYM tracks nitrogen and phosphorus in algae (IP, IN), organic matter (POP, PON), and dissolved inorganic forms available for algal uptake (FRP, NH4 + , NO3 − ) to ensure mass balance of nutrients. In this application, constant internal nutrient to algal biomass ratios (YP:chla , YN:chla ) were modeled, though CAEDYM can also be configured to simulate dynamic nutrient to chla ratios. All of the nutrient state variables are explicitly tracked to ensure mass balance, regardless of whether dynamic or static internal nutrients are modeled. PON and POP were modeled as one detrital particle with the same parameters for decomposition and settling (Table 3). In short, the simplest nutrient cycle (nutrient → algae → detritus → nutrient in Fig. 2) is examined as an ecological representation of the meso-oligotrophic reservoirs considered in this study. Many of the model parameters in this application were ‘fixed’ from previous studies though several had to be calibrated over the three model settings considered here (Table 3).

148

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

Table 2 A summary of the differential equations for the nitrogen, phosphorus, dissolved oxygen, and algal cycles in CAEDYM Nitrogen

    D,G,B
∂NO3 − DO Kd =− µi (1 − PNi )YN:chla Ai + kn f(T) NH4 + − kd f(T) NO3 − ∂t Kn + DO Kd + DO i          nitrification

algal NO3 uptake

∂NH4 + =− ∂t

D,G,B




i

denitrification







 f(T)SNH4 KDOS /(KDOS + DO) DO kAN KON + kON DO µi PNi YN:chla Ai − kn f(T) NH4 + + f(T)PON + Kn + DO KON + DO h            nitrification

algal NH4 uptake





organic matter mineralization

dissolved sediment flux



D,G,B
vN ∂PON kAN KON + kON DO = PON − kr f(T)YN:chla Ai f(T)PON + ∂t z KON + DO        i   settling

 PNi =

organic matter mineralization

NH4 + NO3 − + (NH4 + KNi )(NO3 −

 + KNi )

mortality and excretion

NH4 + KNi + (NH4 + NO3 − )(NO3 − + KNi )



Phosphorus

  D,G,B
∂FRP kAP KOP + kOP DO f(T)SP KDOS /(KDOS + DO) µi YP:chla Ai + =− f(T)POP + ∂t KOP + DO h   i        dissolved sediment flux organic matter mineralization

algal FRP uptake

  D,G,B
∂POP vP kAP KOP + kOP DO = POP − kr f(T)YP:chla Ai f(T)POP + ∂t z KOP + DO        i   settling

organic matter mineralization

mortality and excretion

Dissolved oxygen D,G,B
∂DO DO = ka (Oa − DOw ) + {µi (1 − kp ) − kri f(T)}YO:C YC:chla Ai − f(CDO ) − kn f(T) NH4 + YO:N       ∂t Kn + DO i    air –water flux mineral    algal oxygen consumption/production

−

nitrification

f(T)FSOD DO/(DO + KSOD ) h    sediment oxygen demand

Phytoplankton group Ai where i = D for diatoms, i = G for green algae and i = B for blue-green algae          ∂Ai vi  αA (τ − τcA )/τr Aimass /(Aimass + Kimass ) = µi − kri f(T) + Ai +      ∂t z  h         growth  respiration

settling

resuspension

fi (T) = θ T −20 + θ k(T −a) + b   NH4 + + NO3 − I f(N) = f(I) = 1 − exp − NH4 + + NO3 − + KNi IK

µi = fi (T)µmaxi min[f(I), f(N), f(P)] f(P) =

FRP FRP + KPi

f(T) = θ T −20

The rates of change of these state variables from transport and mixing are simulated by the hydrodynamic drivers (1D DYRESM, 3D CAEDYM), and are not listed here. Algal groups include diatoms, chlorophytes (green algae), and cyanophytes (blue-green algae), which were configured with constant internal nutrient to chla ratios. The variable h is the layer (1D) or grid cell (3D) thickness above the sediments and z is the layer or grid cell thickness in the water column. DOw is the dissolved oxygen in the surface layer. The term f(CDO ) represents dissolved oxygen consumption from organic matter decomposition (DO sub-model). The organic carbon sub-model is not described here because it is beyond the scope of this study but is detailed in Romero et al. (2003). PNi is the NH4 preference factor of algal group i over NO3 . τ r is a reference shear stress set to 1 N m−2 .

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

149

Table 3 Nitrogen and phosphorus parameters used in the simulations with all rates at 20 ◦ C and temperature multipliers (θ) set to 1.08 (Ambrose et al., 1988), except for algal temperature functions Parameter

Value

Description

References/remarks

DO Oa ka KDOS

Equation Equation 0.5

Equivalent DO at air–water interface (mg DO l−1 ) Transfer coefficient from equation DO 1/2 saturation constant for nutrient sediment fluxes (mg DO l−1 ) Static DO consumption by sediments (g DO m2 per day) DO 1/2 saturation constant for sediment oxygen demand (mg DO l−1 )

Riley and Skirrow (1975) Wanninkhof (1992) Pickering (1994)

FSOD

0.3

KSOD

5.0

Pickering (1994) Calibrated

FRP, POP YP:chla kOP kAP vP a SP KOP

0.3 0.010 0.003 −0.1 0.0005 5.0

Fixed algal P to chla ratio (mg P (mg chla)−1 ) Aerobic POP mineralization rate (per day) Anaerobic POP mineralization rate (per day) Settling velocity of POP (m per day) Maximum FRP sediment flux (g P m−2 per day) DO 1/2 saturation constant for POP mineralization (mg DO l−1 )

Stoichiometry relation Schladow and Hamilton (1997) 30% of aerobic rate, as kAN Calibrated Pickering (1994) Calibrated

NO3 , NH4 , PON YN:chla kON kAN vN b SNH4 YO:N

9.0 0.010 0.003 −0.1 0.01 3.43

Fixed algal N to chla ratio (mg N (mg chla)−1 ) Aerobic PON mineralization rate (d−1 ) Anaerobic PON mineralization rate (per day) Settling velocity of PON (m per day) Maximum NH4 sediment flux (g N m−2 per day) Nitrification stoichiometry ratio of DO to N (mg DO (mg N)−1 ) Nitrification rate (per day) DO 1/2 saturation constant for nitrification (mg DO l−1 ) Denitrification rate (per day) DO 1/2 saturation constant for denitrification (mg DO l−1 ) DO 1/2 saturation constant for PON mineralization (mg DO l−1 )

Stoichiometry relation As kOP As kAP As vP Pickering (1994) Stoichiometry relation

kn Kn kd Kd

0.02 2.0 0.01 0.5

KON

5.0

Diatoms (D), green algae (G), blue-green algae (B) YO:C 2.67 Photosynthetic stoichiometry ratio of DO to C (mg DO (mg C)−1 ) 40 Ratio of C to chla (mg C (mg chla)−1 ) YC:chla 0.014 Fraction of algal DO lost to photosynthetic respiration kp 1.3, 0.8, 1.1 Maximum growth rates of algae (per day) µmaxD,G,B krD,G,B 0.14, 0.09, 0.07 Algal respiration, mortality, and excretion (per day) 0.007, 0.005, 0.008 P 1/2 saturation constant for algal uptake (mg P l−1 ) KPD,G,B KND,G,B IKD,G,B

0.06, 0.06, 0.03 60c , 100, 130

vD,G,B TSTD,G,B TOTD,G,B TMTD,G,B kD,G,B aD,G,B

−0.05, −0.02, −0.01 16, 20, 20 20, 28, 30 29, 35, 39 2.4, 3.2, 2.2 24.7, 30.2, 30.1

N 1/2 saturation constant for algal uptake (mg N l−1 ) Light 1/2 saturation constant for algal limitation (␮E m−2 s−1 ) Algal settling velocities (m per day) Standard temperature for algal growth (◦ C) Optimum temperature for algal growth (◦ C) Maximum temperature for algal growth (◦ C) Exponent parameter 1 for algal temperature function Exponent parameter 2 for algal temperature function

Hamilton and Schladow (1997) Calibrated Calibrated As SP As KOP

Stoichiometry relation Griffin et al. (2001) Ambrose et al. (1993) USCE (1995) Schladow and Hamilton (1997) 0.002–0.05 (USCE, 1995), calibrated Schladow and Hamilton (1997) Hamilton and Schladow (1997) Calibrated Griffin et al. (2001) Griffin et al. (2001) Griffin et al. (2001) Computed Computed

150

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

Table 3 (Continued ) Parameter

Value

Description

References/remarks Computed Low critical shear stress ensures algal resuspension in the upper water column High so rate not limiting Low, so all available algae resuspended

bD,G,B τ CD,G,B

0.20, 0.08, 0.18 0.001

Offset parameter for algal temperature function Critical shear stress for algal resuspension (N m−2 )

αD,G,B KmassD,G,B

0.000008 0.00001

Resuspension rate of algae (g chla m−2 s−1 ) 1/2 saturation constant of available phytoplankton mass on sediments for algal resuspension (g chla m−2 )

The parameters (a, b, and k) for the algal growth temperature function were determined so that maximum productivity occurred at TO (optimum temperature), zero productivity was simulated at the lethal temperature of TM (maximum temperature), and the power relation proportional to θ T −20 was modeled below TS (the standard temperature). Higher PON and POP settling rates (vP and vN ) for the 3D simulations are footnoted. ‘Calibrated’ parameters in italics were varied during simulations to provide good validation results for both reservoirs. All other parameters were assigned values based on previous studies and a literature review. a PON settling velocity set to −2 m per day for the 3D flood simulation. b POP settling velocity set to −5 m per day for the 3D flood simulation. c Ritchey and Romero (unpublished data).

4. Simulation inputs Meteorology from both reservoirs was not available over the simulation period, so measurements of short-

wave radiation, net longwave radiation, air temperature, and vapour pressure from on-lake stations during 2002–2003 were used. Daily wind run and rainfall were available from the shoreline of Prospect Reser-

Fig. 3. (A) Gauged daily discharge from Warragamba Pipeline and Upper Canal, and total withdrawals of Prospect Reservoir from 1988 to 1990. Measured daily to weekly (B) temperature, (C) FRP, (D) TP, (E) NH4 + , (F) NO3 − , and (G) TN.

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

voir during 1988–1990 and nearby weather stations for Lake Burragorang during 1992–1994, and served as inputs. In Prospect Reservoir from February 1988 to February 1990, both the Upper Canal (200 ML per day) and the Pipeline (1000 ML per day) discharged directly into the reservoir (Fig. 3A). Both of these inflows were sampled frequently (daily to weekly) for temperature, dissolved oxygen, and nutrients (Fig. 3). Temperature (Fig. 3C) and phosphorus levels (Fig. 3C and D) in the Pipeline were higher in the second year because withdrawals were extracted from the metalimnion (1989–1990) rather than the hypolimnion (1988–1989) of Lake Burragorang. The high NH4 of the Upper Canal was from chloramination prior to discharge into Prospect Reservoir (Fig. 3E). From July 1992 to July 1994 drought conditions in Lake Burragorang led to a low average inflow of 500 ML per day (Fig. 4A) relative to the long-term 3-year average of ca. 4000 ML per day from 1960 to

151

2002. The residence time over the 2 years of the simulation period was 8 years, substantially longer than 1 month of Prospect Reservoir. Withdrawals of 1000 ML per day at 30–50 m below the surface corresponded to the upper hypolimnion during seasonal thermal stratification (Fig. 4B). Inflow data for the rivers near the confluence with Lake Burragorang is temporally sparse because of the remote and rugged nature of the lower catchments. Inflow biogeochemistry was estimated from correlations with discharge and water temperatures for the seven rivers (Fig. 4C–G). No substantial floods occurred over these two drought years, so details of these inflow estimates are not provided here. A well-monitored flood event from late June to mid-July 1997 in Lake Burragorang served as the third case of this study for a 3D simulation with ELCOM–CAEDYM. A detailed description of the field study and validation of ELCOM–CAEDYM for this flood is given elsewhere (Romero and Imberger,

Fig. 4. Daily discharge from July 1992 to July 1994 of (A) two of the primary inflows (Wollondilly and Kowmung Rivers) into Lake Burragorang and (B) withdrawals at the dam. (C) Measured water temperatures at the Wollondilly gauging station, and representative estimates of nutrient inputs for the Wollondilly and Kowmung Rivers of (D) FRP, (E) TP, (F) NH4 + , (G) NO3 − , and (H) TN.

152

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

Fig. 5. (A) Gauged discharge of the major rivers during the 1997 flood event and examples of the measurements from the Wollondilly + River station including (B) river temperature, (C) FRP, and (D) PON(= TN − NO− 3 − NH4 ). Lines on FRP and PON plots are based on regression relations as a function of discharge.

2003). Approximately 90% of the discharge during this flood event was from the Wollondilly River, which had a gauging station and automatic water sampler. Hence the total fluxes of water (Fig. 5A), heat (Fig. 5B), and nutrients (Fig. 5C and D) into the reservoir were well measured for this flood event.

5. Results and discussion 5.1. 1D simulation of Prospect Reservoir The stratification and biogeochemical dynamics of Prospect Reservoir were largely controlled by external forcing from the inflows (Warragamba Pipeline and Upper Canal with 80 and 20% of the discharge, respectively) and withdrawals (supply to the city of Sydney) that led to a low residence time (ca. 30 days). DYRESM–CAEDYM reproduced thermal stratification well over the 2 years (Fig. 6A) with hypolimnetic temperatures maintained at ca. 15 ◦ C by the Warragamba Pipeline inflows (Fig. 3B). Constant Pipeline temperatures occurred over the 2 years because the source of these inflows was withdrawals from the upper hypolimnion of Lake Burragorang (consistent tem-

perature) with subsequent transport via the Pipeline (prevention of atmospheric heating). In contrast, the Upper Canal inflows had strong temperature seasonality (Fig. 3B), therefore these inputs inserted above the hypolimnion during thermal stratification (Fig. 6A). High DO in both inflows (not shown) maintained oxic conditions in both the surface and bottom waters in correspondence with field profiles (Fig. 6B). The simulated total phosphorus (TP = FRP + POP + IP, where IP = YP:chla chla) and FRP reproduced the field data, except the modeled FRP was too low in the summer mixed layer (Fig. 6C and D). Variations of FRP and TP in the reservoir tracked external phosphorus inputs of the two major inflows (Fig. 3C and D), which again highlighted the dominance of external forcing on Prospect Reservoir. NH4 + , NO3 − , and total nitrogen (TN = NH4 + + NO3 − + PON + IN, where IN = YN:chla chla) also followed the observed seasonal patterns (Fig. 6E–G). Though reservoir NO3 − levels (Fig. 6E) tracked the concentrations of the two inflows (Fig. 3F), the high NH4 + of the Upper Canal (Fig. 3E) was substantially diluted upon discharge into the reservoir (Fig. 6F). The phytoplankton biomass (as indexed by chla) was in the range of observations, but was overesti-

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

153

Fig. 6. Comparison between field data (symbols) and 1D simulation (lines) of Prospect Reservoir at 2 m (∗, thin line) and 17 m (䊊, bold line) from 1988 to 1990 for (A) temperature, (B) dissolved oxygen, (C) FRP, (D) TP, (E) NO3 − , (F) NH4 + , (G) TN, and (H) chla.

mated during the spring bloom of the second year (algal growth too high) and underestimated in the mixed layer during summer stratification (modeled FRP too low) (Fig. 6H). As in other temperate reservoirs, diatoms are a dominant group in the winter and spring phytoplankton assemblage of Prospect Reservoir (Romero, unpublished data). Often a spring diatom bloom is silica (Si) limited, an element not considered in this modeling study. However, reactive Si levels in the surface waters were in excess of 3 mg SiO2 l−1 from 1988 to 1990 (Romero, unpublished data), which precludes Si limitation of the spring diatom bloom. In summary, the algal, nitrogen, and phosphorus state variables tracked the observed seasonal variations with the 1D configuration (DYRESM–CAEDYM). Next the same parameter set was applied to a simulation of two drought years for the much larger Lake Burragorang system.

5.2. 1D simulation of Lake Burragorang The simulation of the 1992–1994 drought reproduced the observed seasonal patterns of temperature, dissolved oxygen, and nutrients at the surface and bottom of Lake Burragorang (Fig. 7). Thermal stratification was reproduced well except in the surface waters during autumnal cooling (Fig. 7A), perhaps because of the application of the 2002–2003 meteorological forcing (except for wind speed and rainfall) to the 1992–1994 period. Simulated dissolved oxygen observations in the bottom and surface waters were reproduced well (Fig. 7B) with the sediment parameters from Prospect Reservoir. No observations were available to evaluate FRP over the simulation period, which generally remained low (ca. 0.002 mg P l−1 ) in the surface waters except during winter mixing (Fig. 7C). Simulated TP was greater than observations

154

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

Fig. 7. Comparison between field data (symbols) and 1D simulation (lines) of Lake Burragorang at 3 m (∗, thin line) and 75 m (䊊, bold line) from July 1992 to July 1994 for (A) temperature, (B) dissolved oxygen, (C) FRP (no field data), (D) TP, (E) NO3 − , (F) NH4 + , (G) TN, and (H) chla.

(Fig. 7D), possibly from overestimation of external nutrient loading during several small flood events. Hypolimnetic NO3 − was simulated well with a decrease each year during holomixis and relatively constant levels through seasonal stratification (Fig. 7E). Surface NO3 − decreased during seasonal stratification from algal uptake similar to field data. Ammonium was reproduced well at the surface though concentrations were overestimated in the bottom waters (Fig. 7F). Denitrification in the bottom waters of a small eutrophic lake was hypothesized to maintain constant hypolimnetic levels of dissolved inorganic nitrogen even with sources from the sediments and organic matter decomposition during seasonal thermal stratification (Krivtsov et al., 2001). Sediment fluxes of NH4 + were low in both of the reservoirs of this study (see Fig. 8E), so modeled organic de-

composition rates (the major source of NH4 + , see Fig. 8E) may have been too high in the bottom waters of Lake Burragorang. The inter-annual reduction of TN from 1992/1993 to 1993/1994 was captured by the simulation (Fig. 7G), but the modeled levels were lower than observations during the second year. Uncertainties in the external TN loading from the inflows during small flood events may account for these discrepancies between modeled and observed levels. Simulated chla at the surface (3 m) was greater than the field data, but had similar seasonal patterns (Fig. 7H). Diatoms dominate the spring phytoplankton bloom in Lake Burragorang (Romero, unpublished data), so that inclusion of Si dynamics into CAEDYM may reduce the algal biomass peak via limitation by this third macronutrient. Diatom blooms not only lead to a depletion of Si in the surface waters, but also

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

155

Fig. 8. (A) Cumulative change of FRP at 2 m during the 1D simulation of Prospect Reservoir. The total change (mg l−1 per simulation) over the model runs of 730 days for Prospect Reservoir at 2 and 17 m, and Lake Burragorang at 3 and 75 m are illustrated for each process for the state variables of (B) FRP, (C) POP, (D) NO3 − , (E) NH4 + , and (F) PON. Processes depicted include algal uptake (uptak), organic matter mineralization (miner), dissolved sediment flux (sedms), algal excretion and mortality into the organic matter pool (ex/mo), nitrification (nitri), denitrification (denit), and vertical mixing (hydro).

of phosphorus levels through sedimentation, thereby decreasing phosphorus availability and phytoplankton concentrations in the summer (Krivtsov et al., 2000, 2001). Next, biogeochemical rate output is used to identify the dominant mechanisms of the nitrogen and phosphorus cycles. 5.3. Seasonal biogeochemical fluxes CAEDYM allows one layer (i.e. depth) to be nominated during 1D simulations to output the incremental change to state variables from each fate mechanism during every time step. This process-based output at 3 and 17 m for Prospect Reservoir and 3 and 75 m for Lake Burragorang was used to identify the dominant biogeochemical flux paths in the surface and bottom waters of these two systems. As the same parameter set was used in both simulations of these reservoirs, a comparison of nutrient fluxes may provide insight into reservoir biogeochem-

ical cycling. For example, FRP at 2 m below the surface in Prospect Reservoir was simulated to have algal uptake as the only loss balanced by mineralization of organic matter and hydrodynamics (Fig. 8A). Hydrodynamics includes fluxes from inflows and outflows in addition to transport and mixing. If these results are expressed as the cumulative change in nutrient concentration over the duration of the simulation (2 years, 730 days), then comparisons between reservoirs (Prospect and Burragorang) and depths (surface and bottom waters) can identify the dominant fate mechanisms. Algal uptake was the only loss of FRP in the surface waters of both reservoirs, which was balanced by replenishment through transport (i.e. vertical mixing) and in situ POP mineralization (Fig. 8B). In the bottom waters POP mineralization was the dominant source balanced primarily by dilution during winter holomixis. In the surface waters, algal mortality and excretion was the dominant source of POP and mineralization was a major loss mechanism in both reservoirs (Fig. 8C).

156

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

However, settling losses of POP from the surface waters of Lake Burragorang were large though substantial replenishment occurred during winter holomixis (vertical transport to the surface). In contrast, POP in Prospect Reservoir was much more influenced by the hydrodynamics (inflow and outflows in particular). Nitrogen cycling in these two systems was similar to phosphorus, though the additional mechanisms of nitrification and denitrification also contributed substantial fluxes. NO3 − at the surface of both reservoirs decreased primarily from algal uptake and to a lesser degree denitrification, whereas nitrification was the primary source (Fig. 8D). PON mineralization was the primary source of NH4 + in the surface waters, whereas nitrification and algal uptake were the primary losses (Fig. 8E). The dominant flux paths for PON cycling (Fig. 8F) were as POP (Fig. 8C). Because oxygen levels did not approach anoxia in the bottom waters, dissolved NH4 + (Fig. 8E) and FRP (Fig. 8B) fluxes were negligible from the sediments. In general, the process-based rates for nitrogen cycling in Prospect Reservoir were similar to Lake Burragorang with a few exceptions. Nitrification was substantially greater in Prospect Reservoir because of the high NH4 + loading from the Upper Canal. Settling losses from the surface waters and subsequent replenishment during winter holomixis was greater in Lake Burragorang because of greater dependence on internal cycling during the drought conditions (low inflow, long residence time). Lastly, inflows (i.e. hydrodynamics) increased PON and POP in Prospect Reservoir, whereas in Lake Burragorang consistent withdrawals decreased particulate nutrient levels. 5.4. 3D simulation of a winter flood through Lake Burragorang As in many water supply reservoirs, flood events provide most of the water into Lake Burragorang, which can cause short-term (weeks) degradation in water quality. Because physical processes such as transport, mixing, and settling largely control spatial and temporal distributions of nutrients during floods; accurate modeling of such events requires 3D simulation. Next a 3D simulation of a moderate-sized flood with the same parameter set as the prior 1D simulations is described.

The June 27–28, 1997 flood resulted in a cool nutrient laden underflow that traversed the reservoir in ca. 7 days (Fig. 9, left panels). The underflow clearly dominated the biogeochemistry during and several weeks after the event. The cool well-oxygenated underflow displaced the pre-flood hypolimnion upwards, which resulted in a mid-depth anoxic region. High levels of nutrients from the underflow remained as a ‘pool’ ‘trapped’ in the profundal regions of the reservoir near the dam. Particle settling was evident from the large decrease of PON from the gauging station (Fig. 5C) to the mid-reservoir site with a subsequent modest decrease to the dam (Fig. 9). Such winter underflow events, if cool and turbid, may cause persistent stratification through the winter months (Ferris and Tyler, 1992). Because of the complex bathymetry of Lake Burragorang (Fig. 1C), model idealizations (Fig. 1D) of the actual reservoir morphometry were used to improve simulations and reduce run times. The idealization involves ‘straightening’ the basin morphology so that the Cartesian grid is aligned with the streamwise and cross-stream axes (Romero and Imberger, 2003). A comparison with field profiles at two stations (DWA2 and DWA27 in Fig. 1C and D) served to validate the accuracy of the modeled underflow movement and biogeochemical concentrations. The only modifications to the biogeochemical parameters from the previous 1D simulations were increases to the settling velocities of POP and PON (Table 2). Floods convey particles with different characteristics into Lake Burragorang than internal particle sources such as algal mortality. The 1D drought simulation compared well with field data because the river discharge was low over the 2 years from 1992 to 1994, and thus did not input large loads of riverine PON. In order to model flood events, at least two organic particle size classes are needed to account for differential settling between external (riverine) and internal (algal detritus) sources of particulate nutrients. Comparisons of isopleths of T, DO, and PON between field and simulation profiles at mid-reservoir (DWA27 at 20 km from the dam) and near dam (DWA2 at 0.5 km from the dam) stations illustrate that ELCOM–CAEDYM reproduced the event well (Fig. 9, right panels). In particular, the upward displacement of the anoxic hypolimnion and the temporal and spatial evolution of PON were captured. Next,

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

157

Fig. 9. Field data (left panels) compared to model output (right panels) at two stations (DWA27 and DWA2, Fig. 1) for temperature, DO, and PON during the 1997 flood event.

158

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

Fig. 10. The two dominant mechanisms for changes to the state variables of (A) NO3 (hydrodynamics and nitrification) (B) FRP (hydrodynamics and mineralization), and (C) POP (hydrodynamics and settling) during the 1997 flood simulation at 10 m above the sediment–water interface at DWA27 are represented as the cumulative change over the model run.

quantification of the most important processes is addressed with the biogeochemical rate output from the 3D simulation. Similar to the 1D simulations, a time series grid cell can be nominated through the CAEDYM user interface to output the rate of change of state variables from fate mechanisms during 3D simulations. The time series grid cell was configured at 10 m above the bottom at the mid-reservoir station (ca. DWA27). The hydrodynamics (in particular advective transport) was the primary mechanism that caused variations of NO3 − , FRP, and POP during the flood (Fig. 10). The next dominant processes for NO3 − and FRP was nitrification and mineralization, respectively, though substantially less important than advective transport. For POP, settling contributed to lower concentrations within the head of the underflow as it passed the time series grid cell. In summary, during physical perturbations such as floods, the hydrodynamics cause most of the spatial and temporal variations of state variables. Hence accurate hydrodynamic numerical modeling is a precursor to confident simulation of the nutrient cycling during and after such events.

6. Conclusions There are advantages and disadvantages to using complex mechanistic biogeochemical models such as CAEDYM relative to other modeling approaches. By incorporating most of the dominant mechanisms into a general plankton model, biogeochemical fluxes can be quantified. By extension, this may provide insight

into the development of suitable amelioration strategies for poor water quality settings. A major disadvantage to these complex models is that a sensitivity analysis over the entire parameter space is impossible. Establishment of generic parameter values for general plankton models has remained elusive. The application of CAEDYM with one parameter set coupled to the 1D hydrodynamics model, DYRESM, to these two differing reservoirs provided the following conclusions. (1) A simple ecological configuration (nutrients, particulate organic matter, and phytoplankton) with one parameter set is sufficient to capture most of the seasonal and interannual vertical variability of nitrogen and phosphorus levels in these meso-oligotrophic reservoirs. However, improved parameter estimation for algal growth, particle settling, particle resuspension, and dissolved sediment fluxes is needed. Incorporation of the Si biogeochemical cycle is also needed to capture the spring diatom bloom and sedimentation dynamics that causes a depression in summer phosphorus concentrations (Krivtsov et al., 2000, 2001). (2) ‘Particle closure’, whereby the same parameter values for settling and decomposition of particulate organic nutrients are prescribed, yields agreement with field observations for these two reservoirs under 1D conditions. (3) The settling velocities of the particulate organic matter (i.e. detritus) (Table 3) to simulate the simple nutrient cycle considered here (Fig. 1) were much lower than field measurements (USCE, 1995). This indicates that an additional state vari-

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

able for ‘dissolved organic matter’ is needed in CAEDYM. Application of ELCOM–CAEDYM with the same parameter set to a flood in Lake Burragorang highlighted the importance of accurate hydrodynamic modeling to reproduce temporal and spatial distributions of biogeochemistry during such perturbations. This application also identified a shortcoming of CAEDYM where only one particle size class for particulate organic matter is currently configurable. In the case of floods, the characteristics of particulate organic matter (i.e. density, diameter, decomposition rates) from riverine sources are substantially different from internal sources such as phytoplankton mortality. The use of a ‘multi-system’ and ‘multi-hydrodynamic driver’ validation approach provided insight into appropriate ecological representation of lacustrine systems and parameter values of biogeochemical processes. Improved ecological representation of CAEDYM to allow for several types of particulate and dissolved organic matter and improved shear stress estimates from DYRESM and ELCOM is underway. Future cross-system validation studies will be extended to other lakes and reservoirs to further bound suitable biogeochemical parameter values over a range of lacustrine settings. Acknowledgements Amir Deen (Sydney Catchment Authority, SCA) and Robert Craig (SCA) provided data and coordinated financial support. Assistance in model applications was provided by Ben Hodges (University of Texas, Austin) and Chris Dallimore (CWR) with ELCOM; David Hamilton (CWR) and Matt Hipsey (CWR) with CAEDYM; and Alan Imerito (CWR) with DYRESM. This research has been funded through a Sydney Catchment Authority project. Five anonymous reviewers and the editor of this special issue (Vladimir Krivtsov) are gratefully acknowledged for valuable suggestions and comments to improve the manuscript. References Ambrose, R.B., Wool, T., Connolly, J.P., Schanz, R.W., 1988. WASP4, A Hydrodynamic and Water Quality Model:

159

Model Theory, User’s Manual, and Programmer’s Guide. Environmental Research Laboratory, US EPA, Athens, GA. Ambrose Jr., R.B., Wool, T.A., Martin, J.L., 1993. The Water Quality Analysis Simulation Program, WASP5, Part A: Model Documentation. Environmental Research Laboratory, US EPA, Athens, GA. Cassulli, V., Cheng, R.T., 1992. Semi-implicit finite difference methods for three-dimensional shallow water flow. Int. J. Numer. Meth. Fluids 25, 629–648. Ferris, J.M., Tyler, P.A., 1992. The effects of inflow and outflow on the seasonal behaviour of a stratified reservoir in temperate Australia—a 20 year analysis. Arch. Hydrobiol. 126, 129–162. Gal, G., Imberger, J., Zohary, T., Antenucci, J.P., Anis, A., Rosenberg, T., 2003. Simulating the thermal dynamics of Lake Kinneret. Ecol. Model. 162, 69–86. Griffin, S.L., Herzfeld, M., Hamilton, D.P., 2001. Modelling the impact of zooplankton grazing on phytoplankton biomass during a dinoflagellate bloom in the Swan River Estuary, Western Australia. Ecol. Eng. 16, 373–394. Hamilton, D.P., 1999. Numerical modelling and reservoir management: applications of the DYRESM model. In: Tundisi, J.G., Straskraba, M. (Eds.), Theoretical Reservoir Ecology and its Applications. Backhuys Publishers, pp. 153–173. Hamilton, D.P., Schladow, S.G., 1997. Water quality in lakes and reservoirs. Part I—model description. Ecol. Model. 96, 91–110. Hamilton, D.P., Schladow, S.G., Fisher, I.H., 1995. Controlling the indirect effects of diversions on water quality in an Australian reservoir. Environ. Int. 5, 583–590. Hodges, B.R., Imberger, J., Saggio, A., Winters, K.B., 2000. Modeling basin scale waves in a stratified lake. Limnol. Oceanogr. 45, 1603–1620. Imberger, J., Patterson, J.C., 1981. A dynamic reservoir simulation model, DYRESM: 5. In: Fischer, H.B. (Ed.), Transport Models for Inland and Coastal Waters. Academic Press, New York, pp. 310–361. Imberger, J., Patterson, J.C., 1990. Physical limnology. Adv. Appl. Mech. 27, 303–475. Jørgensen, S.E., Bendoricchio, G., 2001. Fundamentals of Ecological Modelling, third ed. Elsevier, Amsterdam, 530 pp. Krivtsov, V., Bellinger, E., Sigee, D., Corliss, J., 2000. Interrelations between Si and P biogeochemical cycles—a new approach to the solution of the eutrophication problem. Hydrol. Process. 14, 283–295. Krivtsov, V., Sigee, D., Bellinger, E., 2001. A one-year study of the Rostherne Mere ecosystem: seasonal dynamics of water chemistry, plankton, internal nutrient release, and implications for long-term trophic status and overall functioning of the lake. Hydrol. Process. 15, 1489–1506. Laval, B., Hodges, B.R., Imberger, J., 2003. Reducing numerical diffusion effects with a pycnocline filter. J. Hydraul. Eng. (ASCE) 129, 215–224. Leonard, B.P., 1991. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection. Appl. Mech. Eng. 88, 17–74. Pickering, S., 1994. Prospect Reservoir sediment metal and nutrient release studies. Australian Water and Technology Report DSQ11WS, Sydney, Australia.

160

J.R. Romero et al. / Ecological Modelling 174 (2004) 143–160

Riley, J.P., Skirrow, G., 1975. Chemical Oceanography, second ed. Academic Press, London. Riley, M.J., Stefan, H.G., 1988. MINILAKE: A dynamic lake water quality simulation model. Ecol. Model. 43, 155–182. Romero, J.R., Imberger, J., 2003. Effect of a flood underflow on reservoir water quality—data and 3D modeling. Arch. Hydrobiol. 162, 69–86. Romero, J.R., Imberger, J., 1999. Seasonal horizontal gradients of dissolved oxygen in a temperate Australian reservoir. In: Tundisi, J.G., Straskraba, M. (Eds.), Theoretical Reservoir Ecology and its Applications. Backhuys Publishers, pp. 211–226. Romero, J.R., Hipsey, M., Antenucci, J.P., Imberger, J., 2003. CAEDYM Science Guide, Version 2.0. Centre for Water Research, University of Western Australia.

Schladow, S.G., Hamilton, D.P., 1997. Water quality in lakes and reservoirs. Part II—model calibration, sensitivity analysis and application. Ecol. Model. 96, 111–123. Schladow, S.G., Hamilton, D.P., 1995. Effect of major flow diversion on sediment nutrient release in a stratified reservoir. Mar. Freshwater Res. 46, 189–195. Spigel, R.H., Imberger, J., Rayner, K.N., 1986. Modeling the diurnal mixed layer. Limnol. Oceanogr. 31, 533–556. USCE, 1995. CE-QUAL-R1: A Numerical One-Dimensional Model of Reservoir Water Quality; User’s Manual. Instruction Report E-82-1 (revised edition). Department of the Army, US Corps Engineers, Washington, DC, 427 pp. Wanninkhof, R., 1992. Relationship between gas exchange and wind speed over the ocean. J. Geophys. Res. 97, 7373– 7381.