Controls on oxygen dynamics in a riverine salt-wedge estuary – a three-dimensional model of the Yarra River estuary , Australia

Introduction Conclusions References


Introduction
Estuaries provide an important role in the filtering and transformation of nutrients transported from catchment areas into marine environments.Global trends including climate change, increased population, industrialization and agriculture have led to a decline in estuarine health across the world (Kennish, 2001;Vitousek et al., 1997;D'Avanzo and Kremer, 1994).Of particular concern is widespread depletion of oxygen in bottom waters affecting benthic and pelagic habitats, and leading to a decline in biodiversity (Kemp et al., 2009).Reduced oxygen leads to altered nutrient budgets, including increased sediment nutrient release often further contributing to eutrophication and deteriorating water quality (Kennish, 2002;Middelburg and Levin, 2009;Webster and Harris, 2004).Quantifying patterns in estuarine oxygen depletion is a crucial step towards understanding the extent to which changing global trends will affect the essential role of estuaries in nutrient transformation (Pena et al., 2010).
In salt wedge dominated riverine estuaries, patterns of observed anoxia and hypoxia are highly dynamic varying over both diurnal (tidally driven) and seasonal (dependent on river flow and temperature) time scales (Cooper and Brush, 1991;Paerl and Pinckney, 1998;Roberts et al., 2012).Balancing inputs of tidally driven salt wedge intrusion and river inflow sets up stratification patterns that prevent mixing between oxygen depleted bottom waters and the atmosphere.Oxygen depletion depends on the biogeochemical processes of organic matter mineralization in both the water column and sediments and is dependent on temperature.In this way the extent of oxygen depletion depends both on physical circulation and seasonal climatic patterns (Kemp et al., 2009).For oceanic climates typical of estuaries in the south east of Australia, precipitation patterns are sporadic and occur throughout the year leading to episodic Figures

Back Close
Full high flow events both in the warmer summer months and colder winter months.For example, the Yarra River estuary, in the south east of Australia (Fig. 1), is a narrow urban river estuary with tight curvature that is dominated by a 22 km salt-wedge (Beckett et al., 1982).Monthly measurements of oxygen have recorded regions of anoxia propagating the length of the estuary during low flow periods a complete absence of anoxia during peaks in fresh water inflow (Roberts et al., 2012).Continuous data loggers have recorded variations of 80-240 mmol m −3 over a diurnal time scale (Roberts et al., 2012).The studies highlight that oxygen dynamics in riverine estuaries can be high dynamic over a wide range of timescales and differ in their response, relative to larger estuarine systems (Kim et al., 2010;Wang and Justić, 2009).
Of concern is that prolonged periods of oxygen depletion in the Yarra River estuary alter nitrogen cycling pathways and can lead to increased nitrogen load input into Port Philip Bay (Grace et al., 2003).It is therefore vital to unravel the dynamic and complex behaviors of the estuary that lead to patterns of observed oxygen depletion and develop modeling tools able to predict how changing flow regimes, brought about by climate change, land-use change or environmental flow management, impact on patterns of oxygen depletion (Kennish, 2002).
Estuarine models able to support our understanding of essential physical and biogeochemical processes have developed rapidly over the past two decades (e.g.Capet et al., 2013;Cerco and Cole, 1993;Chan et al., 2002;Wang et al., 2011).In particular, numerical models are frequently used to predict the physical processes that affect the development of patterns of hypoxia and anoxia (Chen et al., 2003;Huang et al., 2011;Wang and Justić, 2009;Warner et al., 2005).These models have found that riverine flow (Ralston et al., 2010), tidal forces (Scully et al., 2009) and storm driven mixing (Wang and Justić, 2009) all play a significant role in determining the extent of salt wedge intrusion.Estuarine topography is also a major factor able to mediate the dominance of each driving force (Savenije, 1993).Circulation models have been used to support observations that deep estuaries with moderate river flow and tidal forcing are dominated by baroclinic pressure gradients and in contrast the circulation of shal-Introduction

Conclusions References
Tables Figures

Back Close
Full low estuaries tend to be dominated by riverine and tidal flow (Ralston et al., 2010).More recently, hydrodynamic models have been coupled to biogeochemical models to directly predict oxygen concentrations in response to patterns of stratification and circulation (Borsuk et al., 2001;Capet et al., 2013;Scully, 2013;Sohma et al., 2008).The advantage of coupled models is the ability to differentiate between biotic and abiotic processes to determine the key driving factors behind observed patterns of oxygen depletion.
Whist there has been significant progress made in modeling morphologically complex coastal environments (Chen et al., 2003(Chen et al., , 2007;;Liu et al., 2007;Luyten et al., 2003), simulating the finely resolved pycnocline of a salt wedge riverine estuary with tight curvature has remained a challenge (Kurup et al., 2000;Li et al., 2005;Oey et al., 1985;Warner et al., 2005).Progress has been made with respect to horizontal grid coordinates able to accurately represent more complex estuarine bathymetries, including a move from Cartesian (Simanjuntak et al., 2011) to curvilinear (Burchard et al., 2004) and finite element models (Wang and Justić, 2009).The recent development of the finite volume method that combines finite element grid structures with finite-difference methods has led to improved computational efficiency in estuarine modeling studies (Chen et al., 2003;Wang and Justić, 2009).Most models use sigma-coordinates in the vertical which may be less suited to resolve sharp pycnoclines in a riverine estuary due to problems with numerical diffusion.Depending on the bathymetric complexity of a particular system, the choice of physical model is crucial to the correct simulation of the salt-wedge dynamics and therefore critical in determining the extent of oxygen depletion.
In this study the finite volume 3-D TUFLOW-FV model platform, which uses a hybrid sigma-z co-ordinate scheme in the vertical, is coupled with a biogeochemical model, FABM, to simulate oxygen dynamics within the Yarra River estuary.Specifically, it is the aim of the study to: (a) evaluate model performance in capturing the extent and strength of salt-wedge intrusion and oxygen depletion using a number of model fit metrics, (b) define conditions of river flow, tidal range and topography that govern the Introduction

Conclusions References
Tables Figures

Back Close
Full strength and extent of salt wedge propagation and stratification, (c) compare patterns of stratification against bottom oxygen concentrations to determine the conditions underlying these patterns, and (d) establish the relationship between temperature, flow and tidal range that govern the magnitude and extent of oxygen depletion.The model is assessed against data from a recent field study that measured both salinity and oxygen profiles along a transect of the estuary over a ten month period (Roberts et al., 2012).The knowledge gained from this application is discussed in terms of supporting management of the estuary, such as environmental flow allocation and the impacts of climate change, and provides the basis for simulating the processes of nutrient cycling and assimilation.

Study site
The Yarra River estuary lies within the city of Melbourne and extends from the edge of the harbor at Bolte Bridge ∼ 18.5 km upstream to Dight's Falls (Fig. 1).Dight's Falls acts as a natural upstream boundary with a seasonally variant, perennial, fresh, oxygenated riverine inflow.The estuary is narrow (width from 30-120 m generally wider, downstream narrowing upstream) and shallow (maximum depth 8 m, minimum depth 1 m and average depth 4 m).For the study period 3 September 2009 to 3 September 2010, upstream flow (entering the model domain at Dight's Falls) ranged from minimum 2 m 3 s −1 (January-February) to maximum of 60 m 3 s −1 (September-November) with an average flow of 8 m 3 s −1 (Fig. 2a).This period represented a dry year, with higher than average late spring flows.Residence time during this period, calculated as total volume divided by river flow, ranges from 2.5 h to 22 days dependent on river discharge rate with an average of 6.5 days.A second inflow, Gardiner's Creek enters the Yarra River estuary approximately half way down the model domain (Fig. 1).Flow from Gardiner's Creek is sporadic, highly correlated to local rainfall and characterized Introduction

Conclusions References
Tables Figures

Back Close
Full by long periods of zero flow and short peaks of up to 50 m 3 s −1 (Fig. 2a).The estuary experiences a semi-diurnal tide, average tidal range ∼ 1.4 m (Fig. 2b).The climate for Melbourne is temperate with warm summers (December-February) and cooler winters (June-August), with temperature ranging from 4-47 • C (Fig. 2c).Rain falls throughout the year with greatest fall from August to September (Fig. 2d).
The estuary exhibits a classic salt-wedge circulation (Beckett et al., 1982), with the extent of salt wedge intrusion highly variable and responding on two time scales of hours-days (tide and river flow) and weeks-months (river flow).The semi-diurnal tide drives the salt wedge up and downstream twice per day and variation in upstream river discharge leads to changes in upstream extent of the salt wedge on a seasonal scale.
Monthly measurements of dissolved oxygen spanning the years 2009 to 2011 have shown hypoxia to be a widespread and common phenomenon during periods of low river flow (Roberts et al., 2012).These measurements found evidence for the presence of anoxic bottom waters extending the full extent of the salt wedge up to ∼ 22 km from Port Philip Bay during a period of extreme low flow.

Hydrodynamic platform
TUFLOW-FV, a finite volume 3-D hydrodynamic model (www.TUFLOW.com;BMT WBM, 2013) is used to simulate the mixing and transport of the Yarra River estuary.TUFLOW-FV adopts a flexible mesh (in plan view), consisting of triangular and quadrilateral elements of different size that are suited to simulating areas of complex estuarine morphometries.The vertical mesh discretization has options for sigma-coordinates, or z coordinates, with the latter also allowing multiple surface Lagrangian layers.
The finite volume numerical scheme solves the conservative integral form of the non-linear shallow water equations (NLSWE).The scheme also simulates the advection and transport of scalar constituents such as salinity and temperature as well as the state variables from the coupled biogeochemical model.The equations are solved Introduction

Conclusions References
Tables Figures

Back Close
Full in 3-D with baroclinic coupling from both salinity and temperature using the UNESCO equation of state (Fofonoff and Millard, 1983).Both 1st and 2nd order spatial integration schemes are available.The temporal integration scheme is explicit and uses both mode splitting and dynamically varying timesteps to maximize computational efficiency subject to Courant and Peclet stability constraints (Fig. 3).
Surface momentum exchange and heat dynamics are solved where appropriate meteorological boundary conditions are supplied.Calculation of surface heat exchange accounts for shortwave, long wave, latent and sensible heat.In the current application, turbulent mixing of momentum and scalars has been calculated using the Smagorinsky scheme in a horizontal plane and with the General Ocean Turbulence Model (GOTM) (Umlauf and Burchard, 2003;Umlauf et al., 2005) for vertical mixing.

Biogeochemical model
TUFLOW-FV was dynamically coupled to the Framework of Aquatic Biogeochemical Models (FABM; Bruggeman et al., 2011).The Aquatic Ecosystem Dynamics (AED) biogeochemical model library included within FABM was compiled to simulate the oxygen dynamics and configured to include state variables for dissolved oxygen, particulate and dissolved organic matter.The processes of mineralisation and sedimentation of particulate organic matter, oxygen exchange through the surface water and oxygen flux across the sediment/water interface are represented by the model (Table 1).

Model setup
The model domain comprised a finite volume grid of 397 trapezoidal cells ranging from 213 to 7585 m 2 in area.A single cell in the cross section along the thalweg of the estuary was chosen to reduce computational time and avoid pressure gradient errors introduced due to steep sloping river bank cells.Higher resolution was applied to regions of tight curvature and lower in straight stretches of the river.In the vertical, 4 Lagrangian layers were applied in the surface layer (subject to elevation change) above Introduction

Conclusions References
Tables Figures

Back Close
Full The methods selected for the vertical mixing model included: a second order turbulence mixing model; a dynamic k-epsilon style method for total kinetic energy; a dynamic dissipation rate equation for length scale method; and a constant stability function.Parameter values are included in Table 2.
Meteorological data were supplied by the Bureau of Meteorology and included air temperature, wind speed and direction, dew point temperature, precipitation and total cloud cover (Fig. 2).Hourly measurements of dew point temperature, air temperature and precipitation were taken from the Melbourne Regional Office Station.Relative humidity was calculated from wet bulb and air temperature using (Bolton, 1980).Hourly measurements of wind speed and direction were taken from Melbourne Airport station and converted to wind speed in the x and y direction.Finally cloud cover measured every 3 h from Essendon Airport was input into TUFLOW-FV to calculate net long wave radiation.
Gauged river discharge supplied by Melbourne Water at Dight's Falls and Gardiner's creek every 5 min were used as input flow rates while 5 min inflow concentration of salinity and oxygen were linearly interpolated from weekly sampling.Inflow temperature was determined by taking a 24-point moving average of air temperature, linearly interpolated to 5 min intervals.This procedure was checked against weekly sampled temperatures and found to have strong correlation.Surface water elevations supplied Introduction

Conclusions References
Tables Figures

Back Close
Full Corresponding salinity, temperature and oxygen values were linearly interpolated from monthly sampling data supplied by Monash University (Roberts et al., 2012).The bottom values from the closest profile sampled to Spencer Street Bridge were used in the interpolation.
Parameters for the oxygen model (refer to Table 2) were estimated, taken from literature or based on laboratory studies undertaken on sediment and water column samples along the model domain (Roberts et al., 2012).For determination of sediment oxygen demand parameters, a MATLAB curve fit using a first order Michaelis Menten limitation function to measured sediment flux data from sediment cores was applied (Fig. 4).
Temperature multipliers in the range 0.8-1.2 were used to correct the measured oxygen flux data and the multiplier with the greatest fit (minimum RMSE = 14.1) was chosen (Table 2).

Model evaluation
Water elevation measurements from 4 stations in the domain, (Spencer Street, Burnley, Hawthorn and Johnson Street refer to Fig. 1) were compared to model output.For salinity and oxygen, 9 along channel profile transects of surveys conducted between September 2009 and June 2009 (refer to Roberts et al., 2012) were used.For these data comparisons were made for time series data, field profiles, degree of stratification as well as upstream extent of salt wedge, hypoxia and anoxia.Time series comparisons of surface and bottom salinity and oxygen were made at three representative stations: Morell Bridge (in the down stream end of the estuary); Scotch College (middle of the estuary) and Bridge Road (upstream end of the estuary).
In this study we evaluate five alternative measures of model performance.The purpose of calculating alternative model fit parameters are two fold: to enable comparison to similar modeling studies with no consistency in method of performance evaluation, and secondly because different methods of model evaluation tell us different things Introduction

Conclusions References
Tables Figures

Conclusions References
Tables Figures

Back Close
Full where N is the number of observations, O i and P i , the "i -th" observed and model predicted data and Ō and P the mean observed and model predicted data respectively.NMAE is a measure of the absolute deviation of simulated values from observations, normalized to the mean; a value of zero indicated perfect agreement and greater than zero an average fraction of the discrepancy normalized to the mean.Similarly RMSE is a measure of the average square error with values near zero indicating a close match.MEF is a measure of the square of the deviation of simulated values from observations, normalized to the standard deviation of the observed data.For both MEF and MSS, a maximum value of one indicates perfect fit; zero indicates that the model provides equal predictive skill as assuming mean observed data.Ralston et al. (2010) showed MEF to be an unreliable measure of model fit, however we have included this metric as it is commonly used allowing us to compare model fit with similar estuary modeling studies.The correlation coefficient gives an indication of the linear relationship between observed and predicted data and is the most common measure for assessing aquatic models (Ahronditsis and Brett, 2004).

Model evaluation
The model performed well in measure of fit for surface water elevation in all four gauged stations (Table 3).All correlation coefficients were > 0.95, the lowest measure of fit Introduction

Conclusions References
Tables Figures

Back Close
Full Simulated time series of the surface and bottom concentrations of salinity and oxygen followed both the seasonal and episodic patterns observed in the field measurements (Fig. 5).Comparisons of along channel transects showed good agreement between model and data including the extent of salt-wedge propagation and associated regions of oxygen depletion (Fig. 6).Measurements of model fit for salinity and oxygen varied significantly from position in the water column and region (Table 4).Whilst measures of MSS, MEF and r 2 were greater overall for salinity (0.42/0.85/0.72)than oxygen (0.29/0.80/0.65)measures of NMAE were lower for oxygen (0.28) than salinity (0.47).
For salinity, measures of model fit indicated that the model predicted bottom salinities better than surface and surface oxygen better than bottom (Table 4).The model simulated salinity of the lower reaches of the estuary better than the upstream end, with negative values of MSS indicating that the model predicted worse than using an observational average for the regions > 10 km upstream of Spencer Street Bridge.For simulated oxygen the patterns of model fit were similar although the model performed better than salinity in the upper reaches of the estuary.In fact, for r 2 the highest value was in the region > 15 km from Spencer Street Bridge (r 2 = 0.74).

Comparison of seasonal patterns of stratification and associated hypoxia/anoxia
Seasonal patterns of simulated output describing the extent of the salt wedge and associated anoxia illustrate strong response to both period of continuous flow and episodic high flow events (Fig. 7).The spatial extent of the salinity stratification showed a general response to seasonal patterns of river flow with the greatest extent during the summer months of low flow and lowest during winter months of higher flow.Exceptions to the general pattern occurred in response to changes in river flow.In the summer 9811 Introduction

Conclusions References
Tables Figures

Back Close
Full Although influenced by patterns of river flow and stratification, the response of the spatial extent of low oxygen to episodic high flow events was less pronounced.In addition the spatial/seasonal patterns (Fig. 7) highlight two regions of low oxygen or hypoxic "hot spots".These correspond to the regions from 5-10 km upstream and 12-17 km upstream of Spencer Street Bridge.

Physical controls on patterns of salt wedge intrusion
Two length scales were computed from both model results and field data to represent the extent of salt-wedge intrusion.The first, L 2 , is the distance from the downstream model boundary (at Spencer Street Bridge) to the point at which the salinity of the bottom waters is less than 2 psu.This length scale, often referred to as the salt intrusion length (Fisher, 1974), represents the extremity of the influence of the tidal pulse on the estuary.The second, L 15 , is the distance from the downstream model boundary to the point at which the depth averaged salinity equals 15 psu.This length scale is representative of the point at which the ocean and river flows meet (i.e. they are in equal proportion).A further non-dimensional number used to determine estuary type is a measure of the ratio of strength of river to tidal flow (R/V ), where R is the volume of river flow and V the volume of tidal flow over a single flood tidal period.
To elucidate the various effects of physical drivers on the salt wedge dynamics of the Yarra River estuary, a comparison of L 2 against R (Fig. 8) and L 2 , L 15 against R/V (Fig. 9) was made.There was a significant amount of scatter in the comparison of L 2 to R but a general inverse logarithmic pattern was apparent.The comparison of L 2 , L 15 against R/V highlighted three inflection points corresponding to areas of tight curvature (horseshoe bends) occurring at ∼ 6500, ∼ 9700 and ∼ 15 800 m respectively.This model of salt-wedge intrusion divided the simulated results into 4 zones corresponding to 4 alternative river flow regimes.To obtain a measure of L 2 and L 15 from the Introduction

Conclusions References
Tables Figures

Back Close
Full field data, salinity measurements were interpolated from monthly profiles at approximately 4-6 sites.These crude interpolations of the field data are not accurate to a fine scale resolution but they do fit within the point of inflection salt wedge model (Fig. 9), providing a further system-level validation of the model.For zone 1, during periods of high river flow L 2 extended beyond the first point of inflection, L 15 remained near to zero.For zone 2, as river flow decreased, L 2 extended towards the second point of inflection and L 15 increased towards and slightly beyond the first point of inflection.For zone 3, L 2 steadily increases towards the third inflection point and L 15 remains steady at a point approximately half way between the first and second inflection point whilst the flow ratio (R/V ) stays within 2-5.For zone 4, under low river flow conditions L 2 extends to the upper domain of the estuary at Dights Falls and L 15 rapidly increased beyond the second inflection point towards the third.Deviations from the inflection point model occurred under conditions of high fluvial flow rates (R/V > 10) in zone 3.Under these conditions whilst L 2 extended beyond the second point of inflection L 15 remained less than the first point.

Physical and biogeochemical controls on patterns of hypoxia/anoxia
To differentiate between driving forces of stratification and sediment oxygen demand on patterns of oxygen depletion in the bottom waters, comparisons of L 2 (extent of salt-wedge intrusion) and bottom water temperatures were made against both percent anoxia (DO < 2 g L −1 ; 62.5 mmol m −3 ) and hypoxia (DO < 4 g L −1 ; 125.0 mmol m −3 ) (Fig. 9a and b respectively).Both percent estuarine anoxia and hypoxia were strongly positively correlated with the temperature of the bottom layer whilst conditions of depleted oxygen followed similar patterns of salt wedge intrusion thresholds.
Model results indicated that bottom water hypoxia is almost non-existent when salt intrusion falls below the distance to the first point of inflection described in Sect. the dominant limiting factor with temperatures below 14 • C having less than ∼ 20 % of maximum hypoxia, 14-16 • C from ∼ 20-50 % maximum hypoxia and for temperatures of greater than 22 • C, the maximum extent of hypoxia occurred in the estuary.
For patterns of anoxia, the model results indicate that anoxia in the bottom waters was avoided when the length of the salt intrusion (L 2 ) is below ∼ 8000 m.This point corresponds to the presence of a sill beyond which dense water is trapped resulting in high residence times and leading to anoxic conditions.Again, once the salt intrusion passes the threshold length scale, the extent of anoxia was highly correlated with bottom water temperatures.For temperatures less than ∼ 15 • C anoxia rarely occurred, for temperatures from ∼ 15-17 • C up to ∼ 20 % anoxia only occurred for L 2 >∼ 14 000 m (i.e. when the salt wedge covered a significant portion of the estuary).Similarly, the maximum extent of anoxia occurred at temperatures > 22 • C with an almost linear increase from L 2 = 8000 : 16 500 m.
To test the response of oxygen depletion in the bottom waters to temperature, a comparison was made between model calculated percentage hypoxia and experimentally derived rates of sediment oxygen demand against temperature (Fig. 10).Included is an Arrhenius temperature function to mimic the temperature limitation function included in the model applied to both sediment oxygen demand and oxygen consumption due to the mineralization of organic matter, it takes the form: where a is a constant representing optimal growth and θ is a temperature multiplier.Whilst the scatter in the simulated output of percent hypoxia against temperature can be accounted for by limitations due to antecedent bottom oxygen concentrations, the maximum oxygen depletion for a given temperature clearly follows the Arrhenius temperature function accounted for in the model parameterization of sediment oxygen demand.Since a similar pattern was observed in the experimental estimations of oxygen sediment flux, this validates the inclusion of the Arrhenius temperature function in the simulated equation for sediment oxygen demand.Introduction

Conclusions References
Tables Figures

Back Close
Full mixing is under-predicted.Calibration of mixing parameters and a higher resolution bathymetry may lead to improved (less stratified) predictions of salinity distribution.Measures of model fit for oxygen were generally lower than for salinity and surface water elevations.This difference is consistent with a general trend in modeling of aquatic ecosystems (Arhonditsis and Brett, 2004;Burchard et al., 2006) where model fit measures decrease from the physics, through nutrients to higher order plankton (Allen et al., 2007a).Whilst salinity is dependent only on the hydrodynamic model equations, oxygen concentrations may additionally include error related to biotic factors represented in the model, which are often less controlled and not as easily parameterised.It is interesting to note that for NMAE the model performed better for oxygen than salinity.
Since NMAE is a measure of the absolute error, normalized to the mean the difference in model performance is attributed to the greater relative spread of data for oxygen (standard deviation = 94 % of mean) compared to salinity (standard deviation = 49 % of mean).These results emphasise the need to calculate a number of measures of model to fit to evaluate model performance against variables with a range of variance.
For salinity, higher values of model fit for the bottom layer compared to the surface layer are indicative of a better representation of the salt-wedge propagation and less on mixing in the surface layers.On the other hand for oxygen, higher values of model fit for the surface layer compared to the bottom reflect the dominance of air-water surface exchange at the top, and highlight a need to focus calibration methods on improved parameterisation of sediment oxygen demand.Simulated outputs indicated significant diurnal variations in the position and extent of both the salt wedge and oxygen depleted waters similar to those observed by data loggers employed in 2008 outside the simulated period (Roberts et al., 2012).Whilst a simplified grid of just a few hundred cells gave acceptable levels of model fit related to the propagation and strength of the salt wedge intrusion, a grid with increased resolution in the cross section would allow improved predictability of the fine resolution of the pcynocline and lateral mixing.
For prediction of oxygen concentrations, closeness of model fit were limited by the error associated with parameter estimations and the assumption of homogeneous sed-Introduction

Conclusions References
Tables Figures

Back Close
Full iment oxygen demand.Whilst a best fit Arrhenius oxygen limitation was used to parameterize the oxygen sediment flux parameters using measured data from laboratory studies there was a substantial amount of scatter in the data.This scatter indicates that other limitation factors were at play including sediment organic matter and that the sediment flux is also dependent on sediment type and therefore location in the estuary.Errors associated with the experimental method also contribute to uncertainty in estimation of model parameters and predicted concentrations.
Each model skill equation used to measure the performance of the model results gave slightly different indications of model fit which is why a number of measures is preferable as they tell us different things about model performance (Bennett et al., 2013;Stow et al., 2009).Traditional methods of fit such as r 2 , NMAE and RMSE were useful in determining absolute deviations from observed and calculated for each spatial domain in the vertical and horizontal helped direct focus on further model refinement.Measures of model fit that relates error to variability in observational error, such as MEF and MSS, on the other hand are useful measures of model efficiency that can be compared to similar modeling studies (Allen et al., 2007b;Stow et al., 2009).For this study we found that in regions with low variance in observed data such as the upper reaches of the estuary values of MEF and MSS gave little benefit in evaluation of model performance.On the other hand in the middle regions with more dynamic patterns of salinity and oxygen depletion, MEF and MSS indicated a strong model performance.
The different measures of model performance demonstrated that whilst the general dynamics of the shifting salt-wedge and associated anoxia were well reproduced by the model, refinement of boundary conditions and grid morphometry is required to improve model performance near the domain boundaries.

Physical controls on patterns of salt wedge intrusion
The seasonal variation of the extent of salt-wedge intrusion in response to the hydrological flow regime is typical of Australian estuaries with shallow geomorphology (Eyre, 1998).For these estuaries wet season flows push the salt wedge toward the Introduction

Conclusions References
Tables Figures

Back Close
Full estuarine mouth leaving the bulk of the estuary in a well-mixed regime.This response differs from general theory based on observations made on estuaries in the Northern Hemisphere where maximum stratification occurs during periods of highest flow (Geyer, 2010).Simple models relating the extent of salt water intrusion to rates of river flow, useful to estuarine managers, have long been studied and were reviewed by Savenije (1993).Whilst simple empirical relationships, based on laboratory studies are often limited to homogeneous cross sections, the power-law (L 2 ∼ Q α ) relating salt intrusion length to river discharge (Uncles and Stephens, 1993) can be universally applied and used to compare estuaries.For the Yarra River estuary our estimate of α = −0.35 is greater than the values of −0.2 for the Hudson River estuary (Oey, 1984) and Sumjin River estuary (Shaha and Cho, 2009) and −0.14 for San Francisco Bay (Monismith et al., 2002) but less than the values of −0.42 to −1.0 estimated for the shallow, convergent Guadalquivir estuary (Díez-Minguito et al., 2013).It is interesting to note that Uncles and Stephens (1993) derived a similar value of −0.27 for the Tamar estuary, a shallow, narrow riverine estuary with tight curvature.The stronger dependence of salt-wedge intrusion on flow could be attributable to the riverine geometry (shallow/narrow) of the Yarra River estuary as opposed to the classical bay topography and morphometry (deep/wide).Whist 3-D hydrodynamic models are often computationally intensive and require large amounts of field data, they are useful as a virtual laboratory to explore the dominating factors that force salt wedge intrusion in estuaries (Parsa and Etemad-Shahidi, 2011).In this study the model results, verified by field data, illustrated that the effect of curvature dominated patterns of salt-wedge intrusion.Whist bottom water salinities can push upstream of regions of tight curvature (points of inflection) greater threshold river flows are required to push the salt wedge beyond these points.Under conditions of high river flow, when the salt-wedge is pushed downstream into the first zone the model predicted that the salinity of the bottoms waters may occasionally remain in the second and third zones.This would be expected to occur during the summer months of high stratification when the warmer fresh riverine flow pushes over the more saline Introduction

Conclusions References
Tables Figures

Back Close
Full bottom waters.Whilst comparison of salt intrusion length (L 2 ) and salt wedge intrusion (L 15 ) are rarely reported for riverine estuaries a comparison of the saline intrusion length to river flow in the Guadalquivir estuary showed three distinct zones corresponding to the different processes that dominate each section (Díez-Minguito et al., 2013), demonstrating the utility of this measure.

Physical and biogeochemical controls on patterns of hypoxia/anoxia
Episodic high flow events that occurred following periods of low flow led to a complete breakdown of stratification and reduction in hypoxia.This observation differs from that of larger estuaries such Chesapeake Bay and Narragansett Bay where response time to river discharge for the hypoxic volume is much greater than for stratification with no direct response on a synoptic time scale (Codiga, 2012;Scully, 2013).The presence of two "hot spots" of reduced oxygen concentration in the Yarra River estuary domain highlights the importance of topography in the determination of patterns of oxygen depletion in narrow estuaries.The regions of low oxygen occurred in deeper waters on either side of a region of shallow water surrounding a sill.In this case the deoxygenated water is trapped in the "deep holes" requiring a threshold of flow induced mixing to replenish bottom water oxygen, similar to that observed in the Tone River estuary during periods of low flow (Ishikawa et al., 2004).It is often difficult to differentiate between driving factors that determine patterns of oxygen depletion in estuaries, formed by the combination of physically driven stratification and biological oxygen demand of the water column and sediments, based on often limited spatial observation (Wang and Justić, 2009).In this study we compared percent oxygen depletion against key estuarine system parameters (such as flow, temperature and extent of salt wedge intrusion) to determine dominant forces in the formation of anoxia and hypoxia.Whilst we found no clear correlation between the length of saline intrusion and distribution of anoxia and hypoxia, a parabolic relationship was evident between L 2 and the maximum extent of oxygen depletion.As discussed above, the length of saline intrusion was closely aligned with river flow, however, whilst 9819 Introduction

Conclusions References
Tables Figures

Back Close
Full saline intrusion responds instantly to changes in mixing dynamics, oxygen depletion is a function of both the antecedent oxygen concentrations (i.e. is delayed by the time spent constrained over sediments under stratified conditions allowing oxygen to deplete) and temperature (follows an Arrhenius relationship) so the correlation is not as clear (Borsuk et al., 2001).To take into account antecedent conditions we plotted percent anoxia against temperature separately for each of the four stratification zones and found a strong linear relationship (Fig. 12).The model predicts an almost seven fold increase in gradient from when the antecedent position of the salt-wedge is in zone 1 to zone 4.This roughly equates to an estuarine average 0.9 % increase in anoxia to a 6.7 % increase in anoxia per 1 • C temperature increase dependent on the antecedent salt wedge conditions.There is concern for estuaries worldwide, that global warming may exacerbate regions of oxygen depletion through reduced oxygen solubility, enhanced stratification and reduced winter mixing (Conley et al., 2009;Pena et al., 2010).The close link between patterns of anoxia, temperature and river flow in the Yarra River estuary highlights the need to increase our understanding of how the estuarine health of shallow estuaries may respond to expected changes in climate resulting from global warming.A review of recent climate patterns in southeastern Australia demonstrated a reduction in wet season rainfall and associated inflows as well as increased summer temperatures (Murphy and Timbal, 2008).A study of rainfall trends in the Yarra River catchment found generally decreasing trends with the months of January and June as exception (Barua et al., 2013).Of particular concern was the significant decrease in rainfall observed for the end of Autumn (Barua et al., 2013).In the present study we have found that whilst increased temperatures in the summer months will lead to increased anoxia and hypoxia, a reduction in winter flows will have less effect and a shift to increased spring and summer flows may have the opposite effect in suppressing oxygen depletion.Using the power law to determine position of salt water intrusion (L 2 ) in the four stratification zones and the linear relationships between percent anoxia and temperature for each zone we predicted the expected probability of anoxia under Introduction

Conclusions References
Tables Figures

Back Close
Full a 2 • C temperature increase.We then ran the full hydrodynamic-biogeochemical model increasing air temperature, inflow temperatures and boundary condition temperatures by 2 • C and re-calculated the expected probabilities of anoxia distribution in the estuary.
A comparison of the two methods against the current climate scenario showed that the simple regression analysis method overestimates the amount of anoxia under a 2 • C temperature rise (Fig. 13).The difference can be explained by the subtleties of dependence of antecedent conditions not accounted for in the simple model.This exercise demonstrates how high frequency model output can be used to develop simple empirical models to relate parameters of estuarine health such as anoxia to environmental forcing parameters such as river flows and temperature.These simple empirical models can be useful tools to enable policy makers to make quick management decisions such as determining environmental flow allocation.

Conclusions
Numerical models that simulate both the hydrodynamic and biogeochemical processes in aquatic ecosystems are useful tools in the determination of driving forces behind patterns of oxygen depletion of salt wedge estuaries.In this study we have been able to explain through analysis of model output, patterns of dependence not easily predicted using field studies, which is a key benefit of modeling applications (Arhonditsis et al., 2006).Alternative measures of model fit were assessed and used to resolve future direction of parameterization, calibration and model refinement and development.
The salt-wedge intrusion of the Yarra River estuary was dominated by riverine flow rates, with greatest propagation during extended periods of low flow.Differences in bottom water salt intrusion and water column salt wedge intrusion are separated into four zones where points of inflection correspond to regions of tight curvature and shallow water.Model results indicated that whilst a threshold salt-wedge intrusion is a requirement for oxygen depletion, temperature was the dominating factor in determining the extent and persistence of hypoxia and anoxia in the estuary.This co-dependence

Conclusions References
Tables Figures

Back Close
Full has major implications for expected temperature rises and timing of high flow events.A simple empirical model was developed relating river flow and temperature to extent of anoxia with in the estuary.This model will allow decision makers to determine the minimum environmental flow required to maintain acceptable levels of oxygen given seasonal variation in temperature.Introduction

Conclusions References
Tables Figures

Back Close
Full      Full  Full     • C climate change using regression analysis (refer Fig. 12); red -simulated by 3-D model under +2 • C climate change.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |−1.0 m AHD, and below this datum up to 37 0.2 m z coordinate vertical layers were used.Depths and position data from two hydrographic surveys of the Yarra River estuary supplied by Parks Victoria and Melbourne Water were used to determine the grid bathymetry.The first survey conducted in 2004 was carried out from a boat logging depth and position at a rate of three soundings per metre run and included depths close to but not including the river banks.A second survey in 2009 measured 160 cross sections within the model domain.Data from both surveys were used and interpolated to obtain a single depth for each cell.A bed elevation smoothing function (40 point moving average) was then applied in order to minimise bathymetrically induced numerical diffusion effects on simulated vertical mixing.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | by Melbourne Water, taken every 5 min at the open boundary at Spencer Street Bridge were used to force the tidally driven flow for the open boundary.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5. Correlation coefficient (r) Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | occurred at Burnley Station located just downstream of the first horseshoe bend.Model skill scores were all reasonably high and NMAE reasonably low with the exception of Burnley Station.Patterns of model fit for water level elevation were consistent for all four measures of skill score.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |months the stratification was episodically pushed back to the downstream end of the domain during high flow events.Similarly during the winter months extended periods of low flow result in stratification extending towards the upstream end of the domain.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3.3 (∼ 6500 m).Once L 2 has passed beyond the first point of inflection maximum extent of hypoxia increases on a parabolic arc from ∼ 45 % to almost 100 % when L 2 extends to Dights Falls.At a given salt intrusion length, model results indicate that temperature is 9813 Discussion Paper | Discussion Paper | Discussion Paper | Whilst rigorous analysis of the accuracy or "model fit" of numerical models is considered crucial, particularly when a model is used for predictive management decisions(Stow et al., 2009), quantitative measurements of model fit are rarely presented(Arhonditsis and Brett, 2004;Mooij et al., 2010).A further challenge in the comparison of model performance is a lack of consistent model fit metrics reported.In this study we included five alternative measures of model fit in order to compare against similar estuarine modelling studies.The measures for surface water elevations were generally high and in close agreement with those measured by Ralston et al. (2010 MSS = 0.78-0.92r = 0.99) andWarner et al. (2005 MSS = 0.85-0.95).The greatest deviation from observed data occurred at Burnley (5.7 km) and can be attributed to the over simplification of bathymetry around the island in the grid cells just downstream of the sampling station.In calculating the measures of model fit for surface elevations, the number of sampling points in the time series record (N = 17 524) was relatively high which could account for the observed consistency in model fit measures as a large sample space minimises bias and sensitivity to phase errors.In general, the high measures of model fit for surface elevation were an indication that the grid bathymetry, inflows and tidal boundary forcing were well represented in the model.Much progress has been made in the accurate prediction of stratification in the lower reaches of estuaries(Ralston et al., 2010) with less success in the upper reaches of narrow, shallow estuaries(Sharples et al., 1994).In general for this application, measures of model fit for salinity decreased with upstream distance.Values of MSS and r 2 compared well to Ralston et al. (2010) for bottom salinity (0.70-0.94/0.88-0.99)and surface salinity (0.29-0.86/0.77-0.94)and Warner et al. (2005) for salinity (MSS = 0.85).The slope and offset calculations suggest that the model is under-predicting lower salinities and over predicting higher salinities indicating the degree of vertical Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Nomenclaturet
= time (seconds) h = layer height (m) Note that for atmospheric oxygen exchange this is the height of the surface layer and for all sediment fluxes this is the height of the bottom layer.T = temperature ( • C) O 2 = concentration of dissolved oxygen (mmol O m −3 ) POC = concentration of particulate organic carbon (mmol C m −3 ) k = oxygen transfer coefficient (m s flux of oxygen across the sediment water interface (mmol O m −2 d −1 ) K O 2 sed = half saturation constant for oxygen dependence of sediment oxygen flux (mmol O m −3 ) θ O 2 sed = temperature multiplier for temperature dependence of sediment oxygen flux (dimensionless) ω POC = settling rate of particulate organic matter (m s −1 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Study site (see inset for location), computational grid and bottom depth relative to Australian height datum (m).

Table 1 .
Equations used to parameterize the processes represented in the biogeochemical model.

Table 2 .
Parameters used in the coupled model simulations.