Hydrodynamic controls on oxygen dynamics in a riverine salt wedge estuary, the Yarra River estuary, Australia

Oxygen depletion in coastal and estuarine waters has been increasing rapidly around the globe over the past several decades, leading to decline in water quality and ecological health. In this study we apply a numerical model to understand how salt wedge dynamics, changes in river flow and temperature together control oxygen depletion in a micro-tidal riverine estuary, the Yarra River estuary, Australia. Coupled physical–biogeochemical models have been previously applied to study how hydrodynamics impact upon seasonal hypoxia; however, their application to relatively shallow, narrow riverine estuaries with highly transient patterns of river inputs and sporadic periods of oxygen depletion has remained challenging, largely due to difficulty in accurately simulating salt wedge dynamics in morphologically complex areas. In this study we overcome this issue through application of a flexible mesh 3-D hydrodynamic– biogeochemical model in order to predict the extent of salt wedge intrusion and consequent patterns of oxygen depletion. The extent of the salt wedge responded quickly to the sporadic riverine flows, with the strength of stratification and vertical density gradients heavily influenced by morphological features corresponding to shallow points in regions of tight curvature (“horseshoe” bends). The spatiotemporal patterns of stratification led to the emergence of two “hot spots” of anoxia, the first downstream of a shallow region of tight curvature and the second downstream of a sill. Whilst these areas corresponded to regions of intense stratification, it was found that antecedent conditions related to the placement of the salt wedge played a major role in the recovery of anoxic regions following episodic high flow events. Furthermore, whilst a threshold salt wedge intrusion was a requirement for oxygen depletion, analysis of the results allowed us to quantify the effect of temperature in determining the overall severity and extent of hypoxia and anoxia. Climate warming scenarios highlighted that oxygen depletion is likely to be exacerbated through changes in flow regimes and warming temperatures; however, the increasing risk of hypoxia and anoxia can be mitigated through management of minimum flow allocations and targeted reductions in organic matter loading. A simple statistical model ( R2 > 0.65) is suggested to relate riverine flow and temperature to the extent of estuary-wide anoxia.


Introduction
Estuaries provide an important role in the filtering and transformation of carbon and nutrients transported from catchments into marine environments.Global trends, including climate change, increased urbanization and agriculture, have led to a decline in estuarine health across the world (Kennish, 2002;Vitousek et al., 1997;D'Avanzo and Kremer, 1994).Of particular concern is widespread depletion of oxygen in the bottom waters of estuaries, affecting benthic and pelagic habitats and leading to a decline in biodiversity (Kemp et al., 2009;Vaquer-Sunyer and Duarte, 2008).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).

Published by Copernicus Publications on behalf of the European Geosciences Union.
There is further concern for estuaries worldwide and 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), though it remains unclear how this may manifest within different estuarine types.Quantifying patterns in estuarine oxygen depletion is a crucial step towards understanding the extent to which changing global trends will affect the essential role that estuaries play in biogeochemical cycling and habitat provision.
This study focuses on salt wedge-dominated riverine estuaries, where patterns of observed anoxia and hypoxia are highly dynamic, with large variations over both diurnal and seasonal timescales (Cooper and Brush, 1991;Paerl et al., 1998;Roberts et al., 2012).Oxygen depletion is driven by the biogeochemical processes associated with organic matter mineralization in both the water column and sediments, and is strongly dependent on temperature.However, the ultimate extent depends not only the rate of organic matter loading to the estuary, but also the physical circulation and seasonal climatic patterns, primarily since the balance of input from tidally driven oceanic water and river inflows sets up stratification patterns that prevent mixing between oxygen depleted bottom waters and the atmosphere (Kemp et al., 2009).For oceanic climates, typical of estuaries in the southeast of Australia, precipitation patterns are sporadic and occur throughout the year, leading to episodic high flow events both in the warmer summer months and colder winter months that occur on top of seasonal winter rainfalls and associated higher flow rates.
In our study area in south-eastern Australia, the Yarra River estuary, the salt wedge propagates far inland up the river channel during periods of low flow, leading to regions of anoxia spanning the length of the estuary (Beckett et al., 1982;Butler and Smith, 1985).In contrast, a complete absence of anoxia has been observed in response to high flow events (Roberts et al., 2012).Continuous data loggers have also recorded variations of 80-240 mmol m −3 over a diurnal timescale (Roberts et al., 2012).These observations highlight that oxygen depletion in low volume estuaries with sporadic patterns of flow can be highly dynamic over a wide range of timescales, and thereby differ in their response relative to larger estuarine systems that are more widely reported (Kim et al., 2010;Wang and Justić, 2009).
Of concern for the Yarra River estuary is that (i) prolonged periods of low oxygen can reduce the overall rate of denitrification and can lead to increased nitrogen loads exported into Port Philip Bay (Grace et al., 2003); (ii) low oxygen has led to a loss of integrity of benthic habitat and fish stocks (Coleman et al., 2011), and (iii) recent climate patterns in the region have recorded a reduction in wet season rainfall and associated inflows, along with increased summer temperatures (Murphy and Timbal, 2008;Barua et al., 2013).Provision of environmental flows to the estuary has recently been a management approach to improve estuarine health (Coleman et al., 2011), in tandem with long-term efforts to reduce organic matter loading to the system.The underlying controls on oxygen depletion, however, remain unclear, making it difficult to devise scientifically-based management strategies.This study is therefore motivated by the need to unravel the dynamic and complex behaviours of the estuary that lead to patterns of observed oxygen depletion in order to understand how changing flow regimes, brought about by climate change, land-use change or environmental flow management, impact on patterns of oxygen depletion.
Estuarine models able to simulate the essential physical and biogeochemical processes that drive patterns of oxygen depletion 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;Xu et al., 2012).Models have been used to demonstrate that estuarine topography is a major factor in mitigating the dominant force driving patterns of salt wedge intrusion (Savenije, 2012(Savenije, , 1993)), including riverine flow (Ralston et al., 2010;Xu et al., 2012), tidal forces (Scully et al., 2009) and storm-driven mixing (Wang and Justić, 2009).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 shallow estuaries tend to be dominated by riverine and tidal flow (Ralston et al., 2010).Hydrodynamic models have also been coupled to biogeochemical models to predict oxygen concentrations in response to patterns of stratification and circulation (Borsuk et al., 2001;Scully, 2013;Sohma et al., 2008) and differentiate between biotic and abiotic processes to determine the key driving factors behind observed patterns of oxygen depletion (Bruce et al., 2011;Capet et al., 2013).
In this study, our aim was to simulate oxygen depletion in the Yarra River estuary on a finely resolved spatial and temporal scale required to capture the essential dynamics associated with varied topography and tight curvature.For this purpose, a 3-D flexible mesh model able to resolve the salt wedge dynamics was coupled with a simple oxygen dynamics model.The specific aims were therefore to (a) validate the model approach within the context of simulating oxygen dynamics within a riverine estuary; (b) define conditions of river flow, tidal range and bathymetry that govern the strength and extent of salt wedge propagation and stratification; (c) establish the relationship between temperature, flow and tidal range that govern the magnitude and extent of oxygen depletion; and (d) explore management implications in relation to minimum environmental flow allocation or changes in sediment oxygen demand that would be required to mitigate increased anoxia under predicted climate change scenarios.Based on this analysis, a simple statistical model is presented that may be used to support management efforts to reduce the severity of oxygen depletion.A1.

Study site
The Yarra River estuary is a highly modified urban estuary that lies within the city of Melbourne, Australia.Over the passage of European settlement, the river has been widened, deepened and the course modified to allow for greater shipping accessibility and to protect the city from flooding.Primarily to alleviate flooding downstream, the Upper Yarra Reservoir was constructed in 1957, which reduced flow to ∼ 50 % of it natural levels.Early settlements used the river for sewage and industries have a long history of using the river and tributaries to dump industrial pollutants.The current flows through the system are highly regulated with water extracted to meet urban and agricultural water demands.
Prior to European settlement the estuary would have been classified as a mature riverine estuary (Roy et al., 2001) with the river flow being the dominant energy source and the coastal environment classified as a wave-dominated delta (Gillanders et al., 2011).With the removal of falls ∼ 7 km upstream of the river mouth and reduced river flow, increased inland propagation of the salt wedge has expanded the estuary domain.According to the estuarine classification of Savenije (2012), it is considered a "riverine" estuary, with "delta" shape, an "ideal" tidal influence, "fixed bed" geology and "positive" salinity.The estuarine Richardson number is generally high and sufficient to maintain a salt wedge for most of the year.
The estuary mouth empties into the northern end of Port Philip Bay and extends ∼ 22 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 to 120 m, generally wider downstream, narrowing upstream and shallow, maximum depth 8 m, minimum depth 1 m and average depth 4 m (Fig. 2).A shallow region or sill of ∼ 1.6 m depth was created by rubble left behind from the construction of a bridge, located ∼ 15.4 km upstream of the estuary mouth (Figs. 1 and 2).
For the study period from 3 September 2009 to 3 September 2010, upstream flow (entering the model domain at Dight's Falls) ranged from a minimum 2 m 3 s −1 (January-February) to a maximum of 60 m 3 s −1 (September-November) with an average flow of 8 m 3 s −1 (Fig. 3a).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 estuary approximately half-way down the model domain (Fig. 1).Flow from Gardiner's Creek is sporadic with short peaks of up to 80 m 3 s −1 , which are highly correlated to local rainfall, followed by long periods of zero flow (Fig. 3a).
The micro-tidal estuary experiences a semi-diurnal tide with an average tidal range of 0.65 m.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 timescales of hours-days (tide and river flow) and weeks-months (seasonal changes in river flow).The climate for Melbourne is temperate with warm summers (December-February) and cooler winters (June-August), with temperatures ranging from 4 to 43 • C (Fig. 3c).Rain falls throughout the year with greatest rainfall from August to September (Fig. 3d).
Historical measurements of dissolved oxygen have shown oxygen-depleted bottom waters to be a widespread and common phenomenon during periods of low river flow (Grace et al., 2003;Roberts et al., 2012).The factors that determine oxygen depletion in the Yarra River estuary are driven by a complex interaction of hydrodynamics, temperature, organic matter loads from rivers, creeks and drains, water column biological oxygen demand and chemical oxygen demands as well as sediment oxygen demand (SOD).Whilst no comprehensive study of the predominant causes of oxygen depletion in the Yarra River estuary has been undertaken, recent studies suggest that SOD from the long-term accumulation of organic-rich sediment exacerbated under stratified conditions is a dominant factor (Grace et al., 2003;Roberts et al., 2012).Measurements taken near Dight's Falls showed that following stratification, concentrations of bottom-water dissolved oxygen declined rapidly, indicating that sediment oxygen demand was a main cause of oxygen depletion in this part of the estuary (Coleman et al., 2011).

Hydrodynamic platform
TUFLOW-FV, a finite volume 3-D hydrodynamic model (BMT WBM, 2013), was developed for complex coastal environments and has been 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, and 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 in 3-D with baroclinic coupling from both salinity and temperature using the UNESCO equation of state (Fofonoff and Millard, 1983).Both 1-and 2-order spatial integration schemes are available.The temporal integration scheme is explicit and uses both mode splitting and dynamically varying time steps to maximize computational efficiency subject to Courant and Peclet stability constraints (Fig. 4).
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 Table 1.Equations used to parameterize the processes represented in the biogeochemical model.

Ordinary Differential Equations dO
/h sediment oxygen demand (4) 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.
sed = temperature multiplier for temperature dependence of sediment oxygen flux (dimensionless) ω POC = settling rate of particulate organic matter (m s −1 ) Rmax POC miner = maximum rate of mineralisation of particulate organic carbon (/d) K miner = half saturation constant for oxygen dependence of mineralisation rate (mmol O m −3 ) θ miner = temperature multiplier for temperature dependence of mineralisation rate (dimensionless) 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 mineralization 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 starts with the downstream boundary condition at Spencer Street Bridge (chainage = 7.4 km) and ends at the upstream boundary condition at Dight's Falls (chainage = 21.9 km).The model domain comprised a flexible mesh 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 riverbank cells.Higher resolution was applied to regions of tight curvature and lower resolution in straight stretches of the river.In the vertical, 4 Lagrangian layers were applied in the surface layer (subject to elevation change) above −1.0m AHD, and below this datum up to 37 z coordinate vertical layers of 0.2 m were used.Australian Height Datum (AHD) is defined as the mean sea level for thirty tide gauges around the coast of the Australian continent.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 mesh 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 riverbanks.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 (40point moving average) was then applied in order to minimize bathymetrically induced numerical diffusion effects on simulated vertical mixing.The methods selected for the vertical mixing model included a second order turbulence mixing model with a dynamic k epsilon method for total kinetic energy, a dynamic dissipation rate equation for length scale method and a constant stability function (refer to details 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. 3).Hourly measurements of dew point temperature, air temperature and precipitation were taken from the Melbourne Regional Office (MRO) station (Fig. 1).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 (MA) station and converted to wind speed in the x and y direction.Finally, cloud cover measured every 3 h from Essendon Airport (EA) station was input into TUFLOW-FV to calculate net long wave radiation.Using meteorological data measured at some distance from the study site (MRO ∼ 1-5 km, MA ∼ 20 km & EA ∼ 10 km), particularly wind data taken ∼ 10 km from the Yarra River estuary, has induced predictive error associated with mixing and stratification.It was the only available data, however, and general diurnal and seasonal trends were captured in the meteorological boundary conditions.
Gauged river discharge supplied by Melbourne Water at Dight's Falls and Gardiner's Creek (Fig. 1) 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 (r 2 = 0.83).This process of interpolation enabled general seasonal variation patterns to be reflected in the inflow boundary conditions but will introduce some error.Surface water elevations supplied by Melbourne Water, taken every 5 min at Spencer Street Bridge, were used to force the tidally driven flow for the downstream open boundary.The water elevation over the study period for Spencer Street Bridge ranged from 0.24 to 1.12 m with an average of 0.59 m (Fig. 3b).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.Similarly, while this method could account for general seasonal variation, it introduced significant error in the downstream open boundary condition.
Parameters for the oxygen model (refer to Table 2) were estimated from the 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 locally taken sediment cores was applied (Fig. 5).Arrhenius temperature multipliers in the range 0.8-1.2 were used to standardize the measured oxygen flux data and the multiplier with the greatest fit (minimum RMSE = 14.1 mmol m −3 ) was chosen (Table 2).

Model evaluation
Since the model is a simplified representation of the estuary, its utility to test our hypotheses is subject to a suitable level of validation.A comprehensive evaluation process was conducted, comparing model simulated output to measured water elevation and salinity and dissolved oxygen profiles.Several methods for model evaluation were conducted based on Arhonditsis and Brett (2004), who highlighted the need for aquatic ecosystem models to be more rigorously assessed.For brevity in the main paper, a full description of the methods, results and discussion of the model evaluation is provided in the Supplement.

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 periods of continuous flow and episodic high flow events (Fig. 6).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 rapid fluctuations in river flow.In the summer 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 resulted in stratification extending towards the upstream end of the domain.
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. 6) highlight two regions of low oxygen or hypoxic "hot spots".These correspond to the regions from 12-16 km upstream and 18-20 km upstream of the river mouth.

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 river mouth (at Port Philip Bay; Fig. 1) 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 upstream extent or toe of the salt wedge.The second, L 15 , is the distance from the ocean to the point at which the depth-averaged salinity equals 15 psu.Assuming a salinity difference of ∼ 34 psu from the ocean water to the fresh water inflow at Dight's Falls, this point represents the cell where approximately half the water in the column has been derived from the ocean and half from fresh water inflows.A further non-dimensional number used to determine estuary type is a measure of the ratio of strength of river to tidal flow, referred to as the estuary number: where Q f is the river inflow rate, t f the tidal period and V t the volume of tidal flow over a single flood tidal period (i.e. the tidal prism).
To elucidate the various effects of physical drivers on the salt wedge dynamics of the Yarra River estuary, comparisons of L 2 against Q f (Fig. 7) and L 2 , L 15 against N (Fig. 8a) were made.There was a significant amount of scatter in the comparison of L 2 to Q f but a general inverse logarithmic pattern was apparent, fitting the curve L 2 = 3.3e 4 Q −0.36 (R 2 = 0.63).The comparison of L 2 , L 15 against N highlighted three inflection points corresponding to shallow regions in areas of tight curvature (horseshoe bends), occurring at ∼ 13.3, ∼ 15.7 and ∼ 20.7 km respectively.These points of inflection have been included in both Fig. 1 (a plan view of the estuary showing curvature) and Fig. 2.These model results of the 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 field data for each of the monthly sampling dates during the study period, salinity measurements were interpolated along the thalweg between profiles taken at approximately 4-6 sites.These crude interpolations of the field data are not accurate to a fine resolution, but nonetheless they fit within the simulated relationship model (Fig. 8a), 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 increased towards the third inflection point and L 15 remained steady at a point approximately half way between the first and second inflection point whilst the estuary number, N, stayed within range 2-5.For zone 4, under low river flow conditions, L 2 extended to the upper domain of the estuary at Dight's Falls and L 15 rapidly increased beyond the second inflection point towards the third.Deviations from the inflection point model occurred under conditions of very high inflow rates (N > 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.
To elucidate which topographical features (i.e.sill or channel curvature) were responsible for the points of inflection and patterns observed in Fig. 8a, two hypothetical scenarios with alternative numerical domains were run: 1.A straightened domain (all curvature removed) with original bathymetry (Fig. 8b) and; Hydrol 2. The original computational domain with the sill removed.This was achieved by increasing the sill depth to 2.8 m (the original depth of the river at that point) and smoothing the surrounding bathymetry (Fig. 8c).
When the curvature was removed from the mesh structure (Fig. 8b), the first inflection point was no longer apparent, while the second and third inflection point remained clear.
On the other hand, when the sill was removed (Fig. 8c) all three inflection points remained, although the second and third were less pronounced.

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 system-wide 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. 9).Both percent estuarine anoxia and hypoxia were strongly positively correlated with the temperature of the bottom layer whilst conditions of depleted oxygen responded to minimum levels of salt wedge intrusion.Model results indicated that bottom water hypoxia is almost non-existent when salt intrusion falls below the distance to the first point of inflection, as described in Sect.3.2 (∼ 13.3 km).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 Dight's Falls.At a given salt intrusion length, model results indicate that temperature is 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 did not occur when L 2 was below ∼ 15.4 km.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 > ∼ 21 km (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 = 15.4 : 23.9 km.
To assess the relationship between oxygen depletion in the bottom waters and temperature, a comparison was made between model-calculated percentage hypoxia and experimentally derived rates of sediment oxygen demand against temperature (Fig. 10).Whilst the scatter in the simulated output of percent hypoxia against temperature can be accounted for by dependence on 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 temperature scaling approach adopted when simulating sediment oxygen demand.

Physical controls on patterns of salt wedge intrusion
The almost immediate response of stratification extent to river flow is typical of Australian estuaries, with shallow geomorphology and hydrological flow regimes that are characterized by low flows interspersed with episodic high flow or flood events (Eyre, 1998).For these estuaries, increased flows push the salt wedge toward the estuarine  mouth, leaving the bulk of the estuary in a well-mixed regime.This response differs from larger estuaries with strong seasonality where maximum stratification occurs during periods of highest flow and mixing during low flow (Geyer, 2010;Scully, 2013).These hydrologically dynamic conditions provide the relatively unique context for which this study was undertaken.
Simple models relating the extent of saltwater 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 or estuaries with ideal shape, 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 (Table 3).For the Yarra River estuary, our estimate of α = −0.35 is similar to the Tamar Estuary, a shallow, narrow riverine estuary with tight curvature and the Swan River estuary in Western Australia that has similar morphometry and hydrological regime.As the variance in power-law exponents can be attributed to morphometric differences between estuaries (Parsa and Etemad-Shahidi, 2011), we have attributed the stronger dependence of salt wedge intrusion on flow to the riverine geometry (shallow/narrow) of the Yarra River estuary as opposed to the classical bay topography and morphometry (deep/wide).
Whilst empirical models are useful tools to determine basic relationships between river flow and salt wedge intrusion, the 3-D model used in this study enabled us to study specific patterns of stratification in the estuary.In particular, we found that the dynamics of the salt wedge were dominated by both bottom topography and the effect of channel  −20) .curvature (Fig. 8).Whilst bottom water salinities can push upstream of shallow regions or 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 main salt wedge is pushed downstream into the first zone, the model predicted that the salinity of the bottom waters would remain in the second and third zones.This occurs during the summer months when strong stratification exists.Exceptions to these observations occurred under conditions of high river flow that we attribute in part to the peak flows from Gardiner's Creek, which, although short in duration, at times exceeded that of the upstream Yarra River inflow rates.Gardiner's Creek is located between the first and second points of inflection where the pattern was disrupted.This supports the importance of storm events causing peak flows in smaller estuary tributaries, which can lead to the disruption of salinity stratification and potential breakdown of anoxic bottom waters in the regions surrounding the point of inflow.
By plotting a direct comparison of salt intrusion length (L 2 ) against the 15 psu halocline (L 15 ) (Fig. 8), we were able to understand the reasons for the high scatter observed in the plot of flow rate against L 2 (Fig. 7).As well as the effect of Gardiner's Creek during peak flow events, antecedent conditions determine how the strength of salt intrusion responds to flow rates.Díez-Minguito et al. (2013) showed a non-linear relationship in their time series plot of L 2 , L 15 and river flow for the Guadalquivir estuary and similarly determined three distinct estuarine zones corresponding to the different processes that dominate each.The results of simulations using alternative mesh structures provided insights as to the cause and effect of physical boundaries to patterns of flow and mixing (Fig. 8).From these simulations it was apparent that curvature had a greater influence on patterns of mixing in the deeper, ocean end of the estuary and shallow regions or sills in the upper reaches of the estuary.The importance of dividing estuaries into zones was a useful exercise in understanding how topographical features dictate how hydrodynamic controls affect position and extent of salt wedge propagation as well as controlling patterns of oxygen depletion, which are described below.

Hydrodynamic controls on patterns of hypoxia and 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 as Chesapeake Bay and Narragansett Bay in the United States, where response time to river discharge for the hypoxic volume is much greater than for stratification with no direct response on a synoptic timescale (Codiga, 2012;Scully, 2013).The presence of two "hot spots" of reduced oxygen concentration in the Yarra River estuary domain highlights the importance of local 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 similar to observations on the Strymon River in the Mediterranean climate of northern Greece (Haralambidou et al., 2010).In these cases the deoxygenated water is trapped in "deep holes", requiring a threshold of flow-induced mixing to replenish bottom water oxygen, similar also to that observed in the Tone River estuary in Japan 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 at the system scale.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.When L 2 was in the fourth zone, anoxia occurred at temperatures as low as 13 • C, consistent with measurements of dissolved oxygen in the upper reaches of the Yarra River (Coleman et al., 2011).While the length of saline intrusion is closely aligned with river flow (the salt wedge responds rapidly to freshwater input), 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 where oxygen is allowed to deplete) and temperature (follows an Arrhenius scaling relationship), so that the correlation is less clear.To our knowledge there is no analytical solution for the extent of anoxia in estuaries in response to changes in flow or other boundary conditions, and so we attempted to develop a multi-variate statistical relationship from the numerical model.Most attempts led to a poor relationship, thus highlighting the multiple interacting factors that ultimately combine to manifest in low oxygen conditions.To demonstrate the effect of antecedent conditions on the response of anoxia to temperature, we plotted percent anoxia against temperature separately for each of the four stratification zones (Fig. 11).For zone 1 there was an almost complete absence of anoxia, and for zones 2-4 we fitted a power relationship with strongest fit when the salt wedge toe was sitting in zones 3 and 4. In practical terms, these results reinforce that, during the summer months when temperatures are higher and antecedent conditions tend to favour extended salt wedge intrusion, a higher flow threshold must be met in order to alleviate anoxic conditions.General management practices tend to increase environmental allocations of freshwater in the winter and reduce flows in the summer to meet increased urban and agricultural demands.Altering these trends to increase storage in winter and send higher flushing pulses of river flow to the estuary in the summer months to drive the salt wedge beyond the critical zones where anoxia risk is highest would help to alleviate oxygen depletion and restore ecological health.

Consequences of increased temperature
As a simple scenario to predict the effect of increased temperatures on patterns of anoxia in the Yarra River estuary, we reran the 3-D model, increasing the temperature of all boundary conditions (air temperature, inflow temperatures and temperature of downstream water level boundary) by 2 • C. Whilst a true climate change scenario would also account for changes in sea level and the timing and intensity of rainfall, run-off and associated inflow rates, this exercise focused on isolating the effect of increased temperature to analyse the sensitivity of oxygen dynamics to increases in SOD with discussion included on how these effects are mediated by discharge.A time series of %anoxia for this "warming" scenario plotted against the "base" scenario showed that during build-up of anoxia, increased temperatures resulted in an increase spatial extent of anoxia but that sporadic river flows in the summer months were sufficient to effectively reset the anoxia independent of the temperature increase (Fig. 12).This exercise demonstrated the usefulness of the model in predicting probability of indicators of estuarine health under increased temperatures but did not take into account shifts in seasonal flow patterns.Due to the high correlation between flow rate and levels of anoxia, predictions of reduced winter flows and a shift to increased spring and summer flows would aid in suppressing oxygen depletion during periods of increased temperature.Future applications of the model will account for changes in flow volumes and seasonal variation in order to determine the net affect of climate change scenarios.
The statistical model can also be a simple tool to understand how estuarine anoxia may respond to flow management and climate change scenarios.To compare predictions of anoxic probability distribution using empirical methods against those made using the 3-D model, we (i) determined the position of saltwater intrusion (L 2 ) based on river flow using the power-law equation derived in Sect.3.2; (ii) divided estimates of L 2 into the four stratification zones and used the empirical relationships between percent anoxia and temperature for each zone to predict %anoxia for each tidal period in the simulated time period; and (iii) calculated the expected probability of anoxia.
The procedure was repeated for the 2 • C temperature increase scenario.The results from both scenarios were compared to the 3-D model results (Fig. 13).This comparison of the two methods showed that the simple regression analysis method overestimates the amount of anoxia both for the current climate and 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, and therefore demonstrates the difficulty in defining simple relationships between flow and oxygen in such dynamics environments.Nonetheless, the simple empirical models can be used by policy makers to make decisions related to estuary management and environmental flow allocation.Any increases in anoxia due to climate change is of concern to the Yarra River estuary for two main reasons: (1) the estuary provides ecological habitat that is reduced by the presence of extensive low oxygen zones; and (2) any flowinduced alteration of hypoxia in the river shifts the rate of system-scale nitrogen assimilation through de-nitrification.A comprehensive management plan for the control of anoxia in the Yarra River estuary would need to include both a reduction strategy for sediment oxygen demand and allocation of minimum environmental flows.As we found that temperature was a strong driving factor for oxygen depletion (for example, by rough estimation if the annual average bottomwater temperature were to increase by 2 • C (17.2 to 19.2 • C), the maximum sediment oxygen demand would have to decrease by ∼ 26 % (∼ 85 mmol m −3 to ∼ 63 mmol m −3 ) to maintain the current levels of hypoxia (Fig. 10)), timing of environmental flows is critical for effective water quality improvement.Whilst this study is specific to the Yarra River estuary, the approach used here highlights the importance of flow management in controlling the extent and duration of anoxia and hypoxia, which can be applied to other estuaries with similar flow regimes and morphometry.

Conclusions
Stratification patterns in the Yarra River estuary were dominated by riverine flow rates, bottom topography and river morphometry.A plot of salt intrusion length (L 2 ) and the 15 psu halocline (L 15 ) suggested that patterns of salt wedge intrusion were separated into four zones, divided by three inflection points corresponding to shallow regions in areas of tight curvature (horseshoe bends) occurring at ∼ 13.3, ∼ 15.7 and ∼ 20.7 km.Model results indicated that whilst a threshold salt wedge intrusion of ∼ 15 km is a requirement for the occurrence of anoxic conditions, temperature was the dominating factor in determining the extent and persistence of hypoxia and anoxia in the estuary.The model indicated that if the antecedent conditions of salt wedge intrusion were in the third or fourth zone, anoxia would occur at temperatures as low as 13 • C and that temperatures greater than 22 • C could result in almost compete coverage of anoxia in bottom waters.This co-dependence 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 the extent of anoxia within the estuary.This model can support decision makers in determining minimum environmental flow requirements to maintain acceptable levels of oxygen under increased temperatures.
The flexible mesh model approach adopted here demonstrated that it was able to simulate the coupled hydrodynamic-biogeochemical processes and understand the driving forces behind patterns of oxygen depletion of salt wedge estuaries.In particular, 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 modelling applications, as defined by Arhonditsis et al. (2006).

Fig. 1 .
Fig. 1.Study site (see inset for location), computational mesh and bottom depth elevation relative to Australian Height Datum (m).Sampling stations described in TableA1.

Fig. 2 .
Fig. 2. Depth and cross sectional area of Yarra River estuary as a function of distance from Port Philip Bay.Vertical black lines represent the 3 points of inflection described in the results.

Fig. 7 .
Fig. 7.Length of saline intrusion (L 2 = length of bottom salinity > 2 psu) against river discharge; * represent field data interpolated along the thalweg from monthly samples taken at 5-6 discrete stations.

Fig. 8 .
Fig. 8. Length of estuary where the depth-averaged salinity = 15 psu (L 15 ) against limit of salinity intrusion (L 2 ) over the range of river to tidal ratios; * represent interpolated field data.Vertical black lines refer to the position of the three points of inflection.(a) Using the original mesh; (b) using a straightened mesh, sill included; and (c) using the original mesh with sill removed.

Fig. 10 .
Fig. 10.Comparison of experimental rates of sediment oxygen demand (*) and area percentage hypoxia (DO < 4 mg L −1 ) (•) against temperature.The solid line represents an Arrhenius temperature function of the form: F max θ (T −20) .

Fig. 11 .
Fig. 11.% Anoxia (%A) as a function of bottom water temperature (T ) for the four zones in the Yarra River estuary.

Fig. 13 .
Fig. 13.Probability distribution of percent anoxia in the Yarra River estuary: blue -simulated by 3-D model; red -estimated by empirical model for (a) under 2009-2010 climate and (b) under +2 • C climate change.

Table 2 .
Parameters used in the coupled model simulations.

Table 3 .
Comparison of power-law exponents for different estuaries.