A Comparison of Bulk Estuarine Turnover Timescales to Particle Tracking Timescales Using a Model of the Yaquina Bay Estuary

The ability to determine a bulk estuarine turnover timescale that is well defined under realistic conditions is in high demand for estuarine research and management. We compare how turnover timescales vary with tidal and river forcing from idealized forcing scenarios using a three-dimensional circulation model of the Yaquina Bay estuary in order to understand the limitations and benefits of different timescale methods for future application. Using model results, we compare bulk formula approaches—the tidal prism method, freshwater fraction method, and a relatively new estuarine timescale calculation method based on the total exchange flow (TEF)—to directly calculated timescales from particle tracking in order to assess the utility of the bulk formula timescales. All of the timescales calculated had similar magnitudes during high river discharge but varied significantly at low discharge and had different dependences on tidal amplitude. Even in the application of a single estuary-averaged timescale, we did not find that any of the bulk timescales described the estuary over a realistic range of tidal and river discharge forcing. During high discharge, the Yaquina Bay timescale is on the order of 2–5 tidal cycles based on the particle tracking analysis, but during low discharge, the turnover time varies across methods and spatial considerations appear to be more important.


Introduction
We compare five turnover timescales from idealized forcing scenarios using a three-dimensional circulation model of the Yaquina Bay estuary. Bulk formula timescales-that is, the freshwater fraction method and a new method for determining estuarine turnover time using the total exchange flow (TEF) at the estuary mouth (MacCready 2011)-are compared to particle tracking methods that estimate transit time and flushing time of the estuary. A third bulk method, the tidal prism method, is briefly discussed, but due to its limited applicability to the Yaquina Bay estuary, only a brief analysis is presented. We examine how these timescales vary with river and tidal forcing as well as how sensitive they are to chosen boundaries that define the estuary volume in order to understand the limitations and benefits of these timescale methods for future application. Under steady state forcing with complete mixing, these timescales would all be equivalent. However, these conditions are not typically met in estuaries and differences between the timescales occur across the range of forcing conditions. The questions that we ask are the following: how do the timescales differ from each other and, specifically, how do the timescales from bulk formulas compare to those directly calculated using particle tracking methods when applied to a numerical model of the Yaquina Bay?
In practice, the ability to calculate a bulk estuarine timescale that is constrained and determinable under realistic conditions is a tool in high demand for estuarine research and management. Turnover time incorporates hydrological and geomorphological information into a single parameter that can readily be related to estuarine ecology by comparison to biogeochemical rates and is therefore proposed as an estuarine classification scheme (Jay et al. 2000). It has been widely applied in water quality, ecosystem management, and budget estimation studies and found to be correlated with distributions and transformations of nitrogen (Wulff et al. 1990;Nixon et al. 1996;Dettmann 2001;Galloway et al. 2003) and other nutrients (Wulff et al. 1990; Lee and Wong 1997;Abril et al. 2002) as well as biomass accumulation (Delesalle and Sournia 1992;Ferguson et al. 2004) across a wide range of estuarine classifications. Across these types of analyses, a range of methods for computing estuarine timescales have been applied, including variations of the tidal prism model (Delesalle and Sournia 1992), the freshwater fraction method (Dettmann 2001;Ferguson et al. 2004), results from hydrodynamics models (Lee and Wong 1997;Dettmann 2001), and the ratio of estuary volume to freshwater inflow (Nixon et al. 1996). Each of these techniques relies on underlying assumptions and simplifications that are important considerations when selecting that timescale for a given application.
The turnover time ratio is the fundamental mathematical basis for the most commonly applied methods for calculating estuarine timescales. The general expression of turnover time of a reservoir is, for scalars, the ratio of the mass in the reservoir to the renewal rate, or in terms of volumetric parcels, the ratio of the reservoir volume to volumetric flow from the sources or sinks of the reservoir (Bolin and Rodhe 1973;Fischer et al. 1979;Monsen et al. 2002). In this analysis, this definition is applied to estuaries, which requires making the assumptions that the volume of an estuary is well defined, that the estuary is in steady state, and that the estuary is well mixed. The sources and sinks for estuaries are generally the rate of freshwater inflow (Q R ) and the volumetric exchange flows at the mouth (Q in,out ), which can be spatially separated from each other (Fig. 1). We discuss three methods that are based on this fundamental theory. The first is a new method proposed by MacCready (2011). The other two are conventional methods in estuarine research: the tidal prism method and the freshwater fraction method.
The definitions and application of estuarine turnover timescales have varied across the literature. Residence time is perhaps the most widely and inconsistently applied of these terms. Bolin and Rodhe (1973) propose that the term residence time can be used to describe the average time for newly incorporated particles (or water masses) to leave an estuary. This umbrella definition applies to many of the uses of the term, such as to refer to localized timescales of transport (Barcena et al. 2012), or as the time for a particle or water parcel to leave an estuary for the first time in contrast to accounting for subsequent reentries due to the tidal nature of estuaries (de Brauwere et al. 2011), but it is not specific to any one of these interpretations. Similar terms such as turnover time and flushing time are often employed, sometimes as synonyms of residence time. Here, for clarity, we specifically define each of the three bulk timescales and refer to them according to the method by which they are calculated. Transit time and flushing time are specifically defined by the particle tracking technique employed. Similar to the definitions used by Vallino and Hopkinson (1998), transit time refers to the time it takes for a passive particle to transverse the estuary from an upriver location, and flushing time refers to the typical time for particles released uniformly throughout the estuary volume to leave the estuary.
In order to avoid ambiguity in the application of any of these timescales, at least the boundaries of the estuary or the local region being flushed, the initial position or concentration of the water (or tracer, or particles) being flushed, and the initial time should be explicitly defined (Monsen et al. 2002;de Brauwere et al. 2011). In this analysis, these parameters as well as the methodology and specific assumptions underlying each timescale are specified.
The new method under investigation is calculated using the total exchange flow (TEF) volume flux quantities, Q in,out , and the corresponding flux-weighted salinities S in,out averaged over a tidal cycle (MacCready 2011). The details of the isohaline analysis used to determine the estuarine exchange flow through a cross section are described by MacCready (2011), who proposes that these results can be applied to calculate an estuarine timescale, referred to here as τ TEF : where the brackets indicate averaging over a tidal cycle (MacCready 2011;Sutherland et al. 2011). The TEF timescale is a bulk timescale based on the general definition of turnover time in a reservoir, which can be thought of either as the ratio of the mass of salt in the estuary to the salt flux into the estuary or as the ratio of salinity-weighted estuary volume (V TEF ) to the volumetric inflow of water from the ocean (Q in ). Defining the estuary in terms of salt provides an unambiguous upriver boundary that accounts for changes in estuary length and salinity structure in a way that a physical location-based boundary cannot. The only ambiguity is in a reasonable selection of a cross section to represent the oceanic boundary. The application of TEF analysis has been shown to lead to a calculation of exchange flow and the transformation of estuarine water that is applicable across estuarine classification regimes (MacCready 2011;Chen et al. 2012). The resulting timescale accounts for all physical tidal and subtidal processes that influence exchange, such as dispersion, that occur within the estuary. The tidal prism method similarly provides a bulk estimation of turnover time within an estuary, τ TP . τ TP is calculated by the ratio of the estuarine volume to the tidal exchange: where P is the tidal prism volume and T T is the tidal period (Zimmerman 1988). The definition of V varies across the literature and can be the estuary volume at low water, at high water, or the mean of the low water and high water volumes. The differences between these definitions are of the order (P). Under the assumption that P<<V, the differences between definitions is small and there is no clear argument in favor of a particular definition of V in the application of this method.
Here, V is defined as the tidally averaged volume between the mouth of the estuary, a (Fig. 2), and an upriver boundary based on the upriver extent of salinity intrusion, b ( Fig. 2; discussed in the "Methods" section): (Fischer et al. 1979). This classical method only requires information about the estuarine geometry and tidal range, but in many applications, P is also multiplied by a factor R TP : that represents the fraction of incoming water that is new ocean water rather than returning estuary water (Ketchum 1951;Sanford et al. 1992). This factor can be determined from the mean salinity of the incoming and outgoing water (Ketchum 1951;Fischer et al. 1979) or from geometric parameters describing the estuary and nearby coastal region (Stommel and Farmer 1952;Chadwick and Largier 1999). The tidal prism approach for determining an estuarine timescale assumes that tidal flow is dominant over river flow, that V is both well defined and much greater than P, and, consequently, that τ TP is much greater than T T . The assumption that V must be well defined is specifically included because V depends upon the selection of boundaries at the mouth and upriver extent of the estuary which are often arbitrarily imposed based on metrics such as tidal excursion or salinity intrusion that vary with forcing over time. This is a practical necessity for the application of this method to estuary environments. The freshwater fraction method calculates the timescale τ FW in estuaries where there is a significant difference between the average estuarine salinity (S AVG ) and oceanic salinity (S O ) by relating the fractional volume of freshwater in the estuary to Q R : where and (Ketchum 1950;Fischer et al. 1979;Pilson 1985;Huang 2007). This method also relies on defining V, which depends on the choice of upriver and downriver limits of the estuarine domain, which are not well defined theoretically. Similarly, the parameter S AVG is dependent on the chosen volume over which salinity is averaged, compounding the influence of how V is defined. The ratio Q R /R FW is an estimate of Q out , the flow rate of estuarine water leaving the estuary to the ocean. This can be shown from the mass balance and salt balance. If the net estuarine salinity is in steady state, the salt balance can be described by ( Fig. 1) and the mass balance by ( Fig. 1). Combining these equations yields Knudsen's relation: The ratio Q R /R FW is amplified relative to Q R because of mixing of oceanic and river water in the estuary. As mixing increases, the exchange between the estuary and the ocean increases and τ FW decreases for a given Q R .
All of the aforementioned timescales inherently rely on the condition that the estuary under question is completely mixed or that the timescale of mixing is much faster than the turnover timescale. This theoretical condition is not realistic for estuaries but is assumed in order to constrain the estimate of turnover time to a welldefined quantity. Under this assumption, if an amount of tracer, C 0 , was introduced into the estuary at an initial time, then the amount of tracer exiting the estuary over time would be determined by the volumetric flow rate out of the estuary, Q out , and the spatially homogeneous concentration of the tracer within the estuary, C/V. The rate of change of the amount of tracer, C, in the estuary could be described by The solution to this equation is an exponential decay, which has a characteristic timescale (in this case, τ FW ) commonly referred to as an e-folding timescale. Under the assumption of complete mixing, it can also be shown that τ TEF , τ FW , and τ TP are all equivalent representations of this characteristic timescale. If S in is estimated by S O , and the estuary volume is assumed to be well defined (such that V TEF =VS AVG /S O ), By making a substitution for Q in from (9) and for Q out from (10), it can be shown through algebra that (12) can be simplified to Similarly, by noting that the tidally averaged volumetric outflow rate is equal to the tidal prism volume divided by the tidal timescale and again assuming complete mixing, it can be shown by substitution that A consequence of the condition of complete mixing is that it is assumed that information about where a tracer (or a water parcel) is introduced into the estuary is immediately lost and that the locations of these sources relative to the sinks do not affect the turnover time (Bolin and Rodhe 1973;Takeoka 1984). If complete mixing is assumed, none of the timescales have any spatial dependence and, for any volume element within the estuary, the characteristic e-folding timescale is the same. For example, the timescale for the salinityweighted volume flushed by the inflow of salty water, τ TEF , is the same as the turnover time of the freshwater volume flushed by the inflow of freshwater, τ FW . There is clearly also an inherent assumption that the volume is well defined here such that the sum of V TEF and VR FW is ∼V. Although these volumes are defined differently for the application of each method, under the steady state and complete mixing assumptions and the theoretical limit that V is well defined, these definitions converge.
Recently, hydrodynamic models have become a more common tool for tracking the fate of tracers or particles released into estuary domains (Oliveira and Baptista 1997;Monsen et al. 2002;Yuan et al. 2007;Liu et al. 2011). These methods have the benefit that they can be used for point-source, regional, or estuary-wide analysis and to investigate spatial variability in mixing and transport. These results can also be used to determine a characteristic timescale for an estuary. Based on the condition of complete mixing that predicts an exponential decay of the amount of a tracer or the number of particles remaining in the estuary over time, this timescale has often been defined as the time that it takes until 1/e, or about 37 % of the initial quantity of a tracer (Yuan et al. 2007) or of the particles (Monsen et al. 2002;Shen and Haas 2004;Liu et al. 2011) remain in the estuary. This method does not make any a priori assumption that the estuary is completely mixed, but if that condition was met, then the number of particles remaining in the estuary over time would theoretically decay at an exponential rate. In which case, the time when 37 % of the particles remained in the estuary would be the same efolding time that was characteristic of the solution to (11). This definition of the characteristic timescale from particle tracking is applied here rather than, for example, a timescale determined from fitting the number of particles in the estuary to an exponential function since it is unrealistic to assume that the estuary would be completely mixed and therefore that an exponential curve would be the best representation of the number of particles remaining in the estuary over time. The 37 % limit is chosen however, because it provides a means to synthesize the particle transport data into a single well-defined timescale which also, under the assumption of complete mixing, converges to τ FW (11) and thus τ TEF and τ TP (13,14).
In this investigation, this particle tracking timescale method is applied to study the transit time from a localized source within the estuary, τ T , as well as the flushing time of particles uniformly distributed throughout the estuary, τ F . These two timescales are different from each other in practice, and the theoretical basis for the differences, which are due to the relative locations of the initial particle locations (sources) to the sink location (a, Fig. 2), is explained in the classical literature (Bolin and Rodhe 1973;Takeoka 1984;Zimmerman 1988). τ T , which is the transit time for particles whose source is far from the sink location (similar to τ FW ), should be greater than the flushing time in theory and is not described by an exponential function. The ways that τ T and τ F differ is of interest here in order to provide insight into the differences between the bulk turnover timescales.
Real estuaries are never completely mixed, estuary volume is ill-defined and changes with forcing, and they are rarely in steady state. Therefore, these timescales are expected to deviate from each other in application. The primary objective of this paper, then, is to compare these five timescales to each other in order to understand the differences in their dependence on river forcing and tidal amplitude as well as their sensitivities to changes in the boundaries used to define the estuary volume. A secondary objective is to compare bulk formula timescales to directly calculated timescales from particle tracking in order to assess the utility of the bulk formula timescales. Our final objective is to gain insight specifically into the applicability of the new turnover time, τ TEF .

Study Area
The Yaquina Bay (Fig. 2) is a macro-tidal drowned-river estuary located on the central Oregon coast of which about 60 % of the area is intertidal (Larned 2003). The estuary area is roughly 11.5 km 2 in a drainage basin about 445.6 km 2 reported by multiple sources (Choi 1975;Sigleo and Frick 2007). There are rock jetties where the mouth meets the Pacific Ocean at Newport, Oregon, through one main channel that is regularly dredged. For the calculations in this study, the oceanward and the landward extent of the estuary are defined by locations a and b, respectively (Fig. 2). Location b, which is about 16.7 km from a, is the upriver extent of the mean salinity intrusion. This distance was estimated by taking the average over all of the simulations of the maximum distance that wat er of salinity great er than or equal to .25 practical salinity units (PSU) was found in the estuary over a tidal cycle at any depth.
The maximum upriver extent of the salinity intrusion estimated by our modeling results varied by about a factor of 2 (from 12.9 to 21.8 km upriver for different river flow estimates). This distance can vary by more than 6 km over a tidal cycle. Due to the relatively small estuary size and high degree of seasonality to the forcing, the Yaquina Bay ranges from a salt-wedge estuary when there is high stratification from large freshwater flow events during the wet season to well mixed during the dry season.
The estuary volume is defined as the volume of water within the boundaries, a and b, of this domain (Fig. 2;Fischer et al. 1979). This volume varies in time, due mostly to tidal fluctuations in water level. The dashed lines in Fig. 2 show shifts in the oceanward and landward boundaries of the estuary that are used for timescale sensitivity analyses. To investigate the dependence of each timescale on the definition of the estuary-ocean boundary, timescales using a (Fig. 2) are compared to timescales calculated from shifting a by +/ −400 m. This boundary shift results in a tidally averaged change in volume of about 2-3 %. A similar analysis is performed by shifting only the upriver boundary (b, Fig. 2) by +/−1 km, which results in a tidally averaged change in volume of about 1-2 % because the channel is much shallower upriver.
The local tides are mixed semidiurnal with a mean range (mean high water to mean low water) of 1.9 m and a mean diurnal range (mean higher high water to mean lower low water) of 2.5 m as reported from NOS station 9435380 at South Beach, Oregon (tidesandcurrents.noaa.gov). The spring and neap tidal ranges are estimated to be 3.21 and 1.58 m, respectively ). The Yaquina river discharge, gaged near Chitwood, Oregon, displays the strong seasonal dependence of river flow as demonstrated in Table 1 (www. oregon.gov/owrd/Pages/sw/index.aspx, station 14306030). Almost 90 % of the annual flow accrues between November and April, which typically occurs in the form of episodic pulses lasting only a few days separated by longer periods of moderate-to-low flow. Two main tributaries converge, the Yaquina River and Elk Creek, and contribute to the total freshwater inflow to Yaquina Bay, Q R (m 3 /s). Elk Creek is not gaged but has a similar drainage area as the Yaquina River, so it is possible to estimate the total freshwater flow into the model domain (e.g., Table 1) by multiplying the Chitwood gaging station measurements by an amplification factor of 2. 52, which accounts for the contribution to the flow from Elk Creek as well as the additional drainage area between the river gage and the convergence Brown and Ozretich 2009).

Numerical Model
The Finite Volume Coastal Ocean Model (FVCOM) is used for this study because of its ability to resolve tidal elevations, water properties, and currents in areas with complex topography such as estuaries with intertidal regions (Chen and Beardsley 2003;Chen et al. 2008). The flexibility of an unstructured grid and the wetting and drying capability of FVCOM for modeling coastal areas and tidal inlets are important features for modeling the dynamics of the Yaquina Bay. FVCOM has been applied and validated in a wide variety of coastal and estuarine regions to study circulation, tidal mixing, and transport (Ralston et al. 2010;Ma et al. 2011;Shaha et al. 2013). The underlying hydrodynamics are based on the three-dimensional primitive equations in sigma coordinates.
The model domain extends 60 km north and south of the estuary and 50 km offshore of the mouth to minimize the effects of the open boundary on the estuarine dynamics. The horizontal resolution within the estuary is 50 m in order to resolve the narrow river channel. Beginning 750 m outside of the jetties, the grid scale increases slowly. At the ocean boundaries, the resolution is 3 km. There are 20 evenly distributed sigma layers in the vertical dimension. FVCOM distinguishes active, wet grid cells and dry cells in the model calculations by a minimum depth criterion for active cells, which is set to .05 m.
The model bathymetry within the estuary is based on the National Oceanic and Atmospheric Administration (NOAA) National Geophysical Data Center (NGDC) Central Oregon Coast digital elevation model (DEM) that was built in 2008 with 1/3 arc second resolution coverage over the Yaquina Bay estuary and coastal region (http://www.ngdc.noaa.gov/dem/ squareCellGrid/download/320). The DEM resolves channel and shoal gradients well enough to fit to the 50-m FVCOM grid, and correlated well with US Army Corps of Engineers surveys collected between 2008 and 2011 at different sites within the estuary. Bathymetry for the offshore region of the grid was gathered from the NOAA NGDC coastal relief model (http://www.ngdc.noaa.gov/mgg/coastal/crm.html).
The General Ocean Turbulence Model (GOTM), a generic turbulence model that can be adapted to reconstruct commonly used turbulence closures within the same format, is incorporated with FVCOM (Burchard 2002;Umlauf and Burchard 2003). For this study, GOTM was applied with the k-ε turbulence closure (Rodi 1976;Canuto et al. 2001). The horizontal eddy viscosity was set to 0, and the background vertical velocity was set to 10 −7 . A spatially uniform bottom roughness length scale was applied to the whole domain. A value of .005 m was chosen, which was found by Frick et al. (2007) to provide the best match to predicted tidal elevations and current data in a previous study of the Yaquina Bay using FVCOM.
The boundary conditions for the model include offshore tidal forcing and freshwater discharge input at the head of the domain. The magnitudes of each of these parameters varied for specific model simulations, which are discussed below. For a set of base model runs, an initial oceanic salinity of 32 PSU is applied to the entire domain, which estimates the mean salinity offshore from the Yaquina Bay (www.ndbc. noaa.gov//station_page.php?station=46094). Riverine water of 0 PSU is input at the head. The model was run for a spinup period until the spatial extent of the estuarine salinity reached tidal equilibrium. The initial spin-up period for higher river flow conditions was on the order of a few days, while it could take weeks for equilibrium to be reached from the initial conditions during low river flow. The model results from previous base runs were used as the initial conditions for subsequent model runs to reduce this spin-up time. Results reported here are from times after the spin-up period for each simulation using model output that is saved for every 414 s of simulation time (108 time steps over an M2 semidiurnal tidal cycle).
Using this model, 20 idealized simulations were run that were identical except for the tidal and river forcing applied in order to investigate the dependence of each timescale on these two parameters. These simulations were idealized in the sense that time-varying external forces such as winds and offshore stratification were excluded. Only the M2 tidal constituent was applied at the offshore boundary as a constant amplitude sine function and river flow was applied at the upriver end of the domain as a steady flow for each simulation. Across simulations, the tidal amplitude and river flow were varied. We used unique combinations of five different river flows, Q R =2, 5, 10, 50, and 100 m 3 /s, and four tidal amplitudes, A T = 1, 1.25, 1.5, and 1.75 m, equal to 1/2 of the tidal range, for a total of 20 model runs. The river flows were chosen to be representative of the seasonal range of the Yaquina Bay from low-flow summertime conditions to large flood events (Table 1). The tidal amplitudes chosen delineate the range across the spring-neap cycle of the local tides.
The model simulations were validated against observations from moored instruments that were deployed between the jetties on September-December 2010 (location c, Fig. 2). At the surface, a conductivity, temperature and depth instrument (CTD) was deployed on 21 September to 04 December to record temperature (T), pressure, depth, and salinity (S). A bottom-mounted acoustic Doppler current profiler (ADCP) recorded the water column currents and sea level from 11 September to 21 December 2010. A CTD was also deployed on the bottom lander, but the conductivity cell provided inconsistent data during the deployment so bottom salinity was calculated based on bottom temperature and the T-S relationship calculated for each day using the time series from the surface CTD. Bottom pressure observations were processed using the t_tide statistics package (Pawlowicz et al. 2002) to fill in data gaps and remove the mean water height from the tidal signal (Fig. 3).
The model results for the idealized scenarios are compared to tidal phase averages of the observed time series during periods of similar tidal and river forcing. A more precise validation, which would require running a realistic version of the model, is beyond the scope of this paper. A comparison of surface and bottom salinity, stratification, and velocity shear in the water column during low river flow (a-b) and high river flow (c-d) is presented in Fig. 4. In order to compare the idealized model results to the observations, periods where river flow of roughly 2 and 50 m 3 /s persisted during spring and neap tides were identified (shaded time periods in Fig. 3). Within each of these periods, there were two to five full higher-high water to higher-high water tidal ranges. The mean of the data over these ranges (with time normalized to tidal cycles) is compared to two tidal cycles from the idealized simulation with a corresponding river flow and the closest tidal amplitude to the higher-high tidal range. Averaging over the diurnal tide demonstrates the comparison across similar tidal and river forcing (the first tidal cycle of each  Fig. 4, which compare the results from periods of spring tides (darker shading) and neap tides (lighter shading) for periods where the river flow was roughly 2 and 50 m 3 /s to the idealized model simulations that match each of these forcing scenarios comparison), and an example of the difference between the idealized sinusoidal tidal forcing and the semi-diurnal forcing. Under low-flow conditions (Fig. 4a, b), the modeled shear and stratification align very closely with the observations for similar forcing. During the second half of the neap semi-diurnal tidal cycle in Fig. 4b, there is slightly less agreement in the shear which can be attributed to the difference in tidal forcing. Under high-flow conditions (Fig. 4c, d), the modeled simulations capture the same order of magnitude of variation in the shear and stratification. The changes in stratification appear to be more abrupt in the model results, but there is still a close agreement, especially during spring tides (Fig. 4d).

Cross-Method Comparisons
We use a power law multiple regression analysis in order to quantify the differences in the timescales calculated in this study. The regression is normalized for each method by a representative river flow, Q 0 =10 m 3 /s, and tidal amplitude, A 0 =125 cm, for the Yaquina Bay estuary such that regression parameters α, β, and γ are calculated by Thus, α is the timescale at the representative river flow and tidal amplitude, and β and γ capture the dependence on river flow and tidal amplitude, respectively. Compared to an exponential or linear regression, the power law regression had higher R 2 values and also displayed a closer statistical fit to all of the timescales across forcing. This power law relation has also been found in previous estuarine timescale analyses (Shen and Haas 2004;Huang 2007). Power law regression analysis has also been used to investigate the influence of tidal and river forcing on other estuarine characteristics such as estuary length Ralston et al. 2010).

Total Exchange Flow Method
Model results were sampled to find the total exchange flow quantities, Q in,out , and S in,out averaged over a tidal cycle for each simulation. Salinity and along-channel velocity were interpolated to 20 evenly spaced locations along a cross section from model results at the estuary mouth (location a in Fig. 2) at roughly .25-m depth intervals. The volume flux through each area element was binned by salinity class. To acquire tidally averaged totals, this process was completed every saved model time step and the mean was calculated over a tidal cycle. Results were found to be sensitive to salinity binning, especially at low river flows where the range of S at the mouth was small. To address this, salinity was divided into 160 bins evenly distributed across the range of salinity encountered at the mouth for a particular simulation. From the total exchange flow quantities, τ TEF was calculated using (1). The sensitivity of τ TEF to the location of the estuary-ocean boundary was tested by calculating the variability between τ TEF using location a and cross sections flanking the channel that are located +/−400 m upriver, which roughly correspond to the up-estuary and down-estuary ends of the jetty constriction (Fig. 2). τ TEF was also compared to a timescale τ TEF Qout , determined from the ratio V TEF /Q out in order to compare its dependence on river flow to τ TEF since exchange flow Q out was observed to have greater river dependence than Q in .

Tidal Prism Method
The tidal prism timescale, τ TP , is calculated from the model results using (2). We define V by (3) and P as the low and high tide volume difference between locations a and b (Fig. 2). Using these definitions, P∼V for all of the simulations, so the criterion that P<<V is not met for the Yaquina Bay. However, it is still possible to determine exchange factors for R TP and τ TP , and a brief analysis was performed to provide some insight into the exchange characteristics of the estuary and the response of τ TP to tidal and river forcing. Using a modified version of the method presented by Fischer et al. (1979), an exchange factor, R TEF , can be calculated using the TEF parameters: Another exchange flow factor for the Yaquina Bay can be calculated from the geometry of the estuary and coastal region using where W is the width and d is the mean depth of the channel at cross section a (Fig. 2), and w is the fraction of a circle that encompasses water from the estuary in the coastal region (Chadwick and Largier 1999). R TEF ranges from about 75 % during low flows to about 90 % during the two highest flow conditions but shows no pattern with tidal forcing, while R geo is mostly constant with river flow but varies from 85 to 90 % with increasing tidal amplitude. During higher river flow conditions, the two exchange flow factors are in greater agreement.
The predicted timescale, τ TP , from (2) is about one tidal cycle (α=1.27) across all forcing scenarios. This result is much lower than predicted by the other methods and is an unreasonable estimate based on the slow response of salinity in the estuary during low discharge. A power law regression analysis of τ TP using (15) shows little dependence on river flow (β=0.02), but as expected, the dependence on tidal amplitude is nearly an inverse relationship (γ=−0.77). When R TP is included (4), γ=−0.83 and −0.75 and α=1.47 and 1.61 from using R TEF and R geo , respectively, but there is almost no change in the weak river flow dependence.

Freshwater Fraction Method
τ FW was determined from (5) and (6) (Ketchum 1950;Fischer et al. 1979;Pilson 1985). Using model results from each simulation, it is possible to calculate the tidally averaged freshwater volume precisely from (5). A value of 32 PSU is used for S o . The mean tidally averaged difference between S 0 and S ranged from 6 to 29 %. Using model output, where values for the freshwater volume and salinity of each grid element throughout a tidal cycle are readily available, the computation of τ FW is straightforward. However, in practice, collecting salinity information with both temporal and spatial coverage requires extensive sampling.
In order to make a comparison to a freshwater fraction timescale that might be estimated from observation, we make an estimate of τ FW obs . Rather than taking full advantage of the extensive spatial and temporal resolution provided by the model fields, τ FW obs is estimated from (5) using parameters that would be available from existing data sources or from bulk estimates. R FW is calculated from (6) where S AVG is the tidally averaged surface salinity for each model run at d (Fig. 2), the location of a Land/Ocean Biogeochemical Observatory (LOBO) from which hourly salinity data is publicly available (http://yaquina.loboviz.com/loboviz.shtml). A tidally averaged mean value V=3.2×10 7 m 3 from all of the estuary runs between downriver and upriver boundaries, a and b (Fig.  2), is used in (5) to determine τ FW obs . The mean value of V varies by about 1 % across runs. The mean tidally averaged difference between S 0 and S AVG ranged from 1 to 50 % across runs; at low river flows, there is not a large difference between estuarine and oceanic salinities, a condition that brings into question the applicability of the freshwater fraction method.
A sensitivity analysis of the selection of an estuary volume was also performed. This was done by first repeating the calculation of τ FW for each run, but using upriver boundaries that are +/−1 km from location b (Fig. 2), and then by repeating the calculation of τ FW for each run, but using downriver boundaries that are +/−400 m from location a (Fig. 2). Fig. 4 A comparison of the mean of an ensemble of observations (narrow lines) to model results (bold lines) at the Yaquina Bay estuary mouth. The ensemble means are found by taking the average results from higher-high water to higher-high water with the timescales adjusted to tidal cycles. Each panel (a-d) includes a comparison over two tidal cycles of river flow (Q R , m 3 /s), stratification between surface and bottom water at the mouth (Δ S, PSU), along-channel velocity shear between the surface and bottom waters (Δ u, m/s), and sea level (η, m). The two panels on the left correspond to spring tides for low (a) and high (c) river flow conditions from Fig. 3, and the panels to the right correspond to neap tides during similar river flow conditions Particle Tracking Timescales Three-dimensional particle tracking was computed using the model-generated hydrodynamic fields. The particle tracking code can use the saved model time step or a given multiple of submodel time steps to advance particles around the model domain using model fields that are linearly interpolated between model time steps. Particle motions are calculated using the Runge-Kutta time stepping scheme that determines integrated particle trajectories to fourth-order accuracy. A vertical random displacement due to vertical turbulent mixing is also included in the particle trajectory that accounts for the spatial non-uniformity of the eddy diffusivity (Visser 1997;North et al. 2006). This vertical displacement is based on the diffusivity generated by the numerical model. To simulate the transport of waterborne particles while accounting for the wetting and drying treatment in FVCOM, elements are considered inactive during a time step if they are dry and particles are required to always remain in active model grid elements. In order to meet this criterion, if a particle is transported from an active cell into an inactive cell, it is reflected back across the inactive cell boundary that it entered from. If a particle remains inside an element that becomes inactive, its position is translated to the center of the nearest wet element.
Two transport timescales determined from this Lagrangian particle tracking scheme are herein referred to as transit time (τ T ) and flushing time (τ F ). In this study, particles are released from their initial positions at eight stages of the tide. Statistics are calculated from the combined particles from all eight of these release times by considering each release time as t 0 =0. The percentage of particles that remain in the estuary is determined for each model time step, and the transit or flushing time is defined as the time when 37 % of the ensemble of particles remain inside of the estuary (Fig. 5). Similar versions of this method of estimating a timescale from modeled tracers have been applied in recent studies (Miller and McPherson 1991;Monsen et al. 2002;Liu et al. 2011).
The difference between the τ T and τ F methods is the initial number and distribution of particles when they are released. To determine τ T , 250 particles were released for each of the eight tidal stages (N=2000) from uniformly distributed locations over an upriver cross section at location b (Fig. 2). τ F was determined from particles initialized at locations within the estuary between locations a and b (Fig. 2) across a uniform 150 m horizontal grid. At each grid location, particles were vertically spaced by 1 m, with the number of particles released dependent on the total water depth at that location and time. Over the eight tidal stages, N ranged from 11,072 to 11,974.
A time series of the number of particles remaining in the estuary over time for the flushing time particle simulations is compared to an exponential curve that is fit to this time series, defined by P(t)=P 0 exp(−t/τ F mix ). P 0 is the number of particles released over eight stages of the tide. P is determined by an exponential fit to the number of particles remaining in the estuary over time while there were greater than 20 % of the particles still within the estuary. The coefficient of determination (R 2 ) is only greater than .1 (.69 to .97) for the runs with 50 and 100 m 3 /s river flow. In the cases with less river forcing, the exponential fit does not describe the pattern of particles exiting the estuary into the ocean. For all of the runs, the number of particles remaining in the estuary initially drops off more rapidly than the exponential fit predicts but flattens out more quickly, crossing the exponential fit line at some point. This comparison for a subset of the runs is shown in Fig. 5a.
A sensitivity analysis was performed for both flushing time and transit time. For this analysis, four forcing simulations (Q R =2 m 3 /s, A T =100 cm; Q R =5 m 3 /s, A T =150 cm; Q R = 50 m 3 /s, A T =125 cm; Q R =100 m 3 /s, A T =175 cm) ranging  (N=11,107; 11,200; 11,287; and 11,366) remaining in the estuary from four flushing time simulations with 125 cm tides and 100, 50, 10, and 2 m 3 /s river flows from left to right (a; thick black lines). The thin gray lines in a show an exponential fit to each of the four flushing time simulations. The percentage of particles (N=2000) remaining in the estuary from four transit time simulations with the same tidal and river forcing as above are shown in b. The dashed gray lines mark 37 % remaining in the estuary. The filled circles in the legend mark τ TEF , τ FW , τ TP , and τ F , on the curves. τ FW for the low-flow case in a is not shown because it is outside of the axis limits, which were selected to show the separation between the four curves from low to high flow across different tidal amplitudes were chosen. For each of these, the particle tracking was repeated 10 times in the exact same way. However, due to the random walk associated with vertical mixing, the resulting particle paths were not identical. Therefore, we verified how sensitive the flushing time and transit time calculations were to this random walk mixing and thus the limited number of particles used in the simulations. Transit time was tested in this way with between N=800 and N=3200 particles. For the strongest forcing case (Q R =100 m 3 /s, A T =175 cm), the standard deviations across 10 repetitions of each of the four forcing simulations varied from about 1 to 2 % of the mean, and for the weakest forcing case (Q R =2 m 3 /s, A T =100 cm), the standard deviations ranged from about 2 to 4 % of the mean depending on the number of particles released. N=2000 was chosen because it was nearly as stable as the simulations with N= 3200, but required much less processing time. For τ T and τ F , the standard deviations across 10 realizations of each of the four simulations tested were within 1-2 % of the mean. From this, we concluded that the timescales are not sensitive to the finite number of particles used.
For all of the simulations, the sensitivity of τ T and τ F to V is determined by comparing the calculation of each of these timescales, to the result from shifting the estuary-ocean boundary by +/−400 m. For τ F , a similar analysis is repeated but involved shifting the upriver boundary by +/−1 km.
Also, a brief investigation into the spatial timescale characteristics was performed. Particle tracking was run until all but 5 % of uniformly distributed particles had left the estuary. Local residence time was determined by the average time that it took for each of the particles released from each grid location to leave the estuary. If a particle had not left the estuary, the maximum number of tidal cycles that the simulation was run for was used in the location average, so the reported timescale is always less than or equal to the local flushing timescale.

Total Exchange Flow Timescale
The TEF terms Q in,out and S in,out demonstrate a greater amount of mixing between ocean and river water within the estuary for high river flows (Fig. 6). Although for all of the simulations the incoming water is almost entirely of oceanic salinity, S out varied with tidal and river forcing. For low-flow conditions, the water leaving the estuary is almost of oceanic salinity (Fig. 6c), while for the same tidal conditions but high river flow, the outgoing salinity classes have a broad range of salinities (Fig. 6a, b). The flow-weighted mean of incoming and outgoing salinities across all of the simulations is plotted in Fig. 6d. τ TEF , relative to a tidal cycle, is reported in Fig. 7a. These results demonstrate dependence on both tidal and river forcing. τ TEF decreases with increasing tidal forcing for a given river flow and also decreases with increasing river flow for a given tidal amplitude. A sensitivity analysis was performed by comparing the reported τ TEF to the values of τ TEF calculated by shifting the estuary-ocean boundary (location a, Fig. 2) by +/−400 m. The differences between these timescales and the reported τ TEF timescale range from roughly −1 to 14 % change and have no apparent trend with forcing. τ TEF Qout has slightly greater river dependence than τ TEF , but the two timescales are comparable ( Table 2).
Examination of the components of the τ TEF equation provides some insight into the source of the variability across forcing scenarios (Fig. 8). The exchange flow Q in increases with tidal amplitude and captures 78-91 % of the variability between simulations with different tidal forcing for a given Q R , but is independent of river forcing (Fig. 8a). The salinityweighted estuary volume, V TEF (1), decreases with both increasing tidal and river forcing (Fig. 8b).

Freshwater Fraction Time
The freshwater fraction timescale, τ FW , has a similar dependence on river and tidal forcing as τ TEF , but there is a greater range from low flow to high flow values and a greater relative dependence on Q R than on A T (Fig. 7b). The freshwater fraction timescales from the upper estuary boundaries defined as 1 km upriver and downriver of b (Fig. 2) vary from the reported τ FW by ∼3-12 %. This sensitivity is greatest for low river flows and for larger tidal amplitudes. The ratio of freshwater to oceanic water accounts for much of the variability because S AVG decreases as the volume is increased. Shifting the estuary-ocean boundary by +/−400 m only impacted τ FW by ∼0.2-0.6 %. For comparison, we also relate τ FW to τ FW obs , which was estimated using values for S AVG (7) that could be readily collected from ongoing observations of the Yaquina Bay estuary (gray lines, Fig. 7b). τ FW and τ FW obs vary by ∼12-140 %.

Transit Time
The transit time, τ T , displays the same qualitative dependence on river and tidal forcing as τ TEF and τ FW although the resulting timescales are considerably larger (Fig. 7c).
There is a relatively small sensitivity across realizations of this method for determining transit time. Across 10 realizations, the standard deviations of the resulting transit times were less than 1-2 % of the mean for each of the four simulations tested. Also, for all 20 of the forcing simulations, transit time sensitivity was calculated from shifting the

Flushing Time
Flushing time results also show an inverse relationship with tidal and river forcing (Fig. 7d). Across 10 realizations of the particle tracking simulation for four forcing simulations, the resulting timescales varied from the mean by ∼1-2 %. The sensitivity of τ F to shifting the upriver boundary by +/−1 km ranges from approximately 0 to 15 % with no evident pattern with forcing. The sensitivity to shifting the estuary-ocean boundary by +/−400 m is roughly half this, ranging from 0 to 6 % of τ F .
Using particle tracking to investigate flushing time, it is also possible to make assessments of the spatial variability of flushing within the estuary. Although thorough analysis of spatial variability is beyond the scope of this study, some of the trends observed are useful to consider in order to understand the bulk estuarine turnover time. For example, it is evident that particles released closer to the estuary mouth remain in the estuary for less time than particles released from further upriver (Fig. 9). Near the mouth, the local flushing time for low flow is similar to the flushing time during high river flow conditions. In the river channel, the difference in local flushing time between low and high river flow is greater and increases further upriver. Figure 9 illustrates this pattern, including for particles released upriver of b, which were not included in the bulk timescale analyses (Fig. 2).

Cross-Method Comparisons
The results of the power law multiple regression analysis (15) are shown in Fig. 7 and Table 2. For all of the methods for determining an estuary timescale applied to the Yaquina Bay in this analysis, the inverse trends with both tidal and river forcing are consistent. The most striking differences are apparent in the details of how each timescale varies with forcing. τ T has the largest variability across forcing conditions. While τ TEF and τ F are similar at representative forcing Q 0 and A 0 , τ F has a much greater range. For all the three methods τ FW , τ T , and τ F , there is a greater range across forcing evident in the timescales than for τ TEF , which varies across all forcing conditions by less than a factor of 3. Most of the timescales are of O(1) at representative forcing, except τ FW and τ T , which have higher values of α.
The TEF timescale, τ TEF , has less dependence on Q R than the predictive timescale, τ F . Both τ FW and τ T have greater dependence on Q R than the flushing timescale, τ F . The greatest river dependence is observed in the transit timescale, τ T . While much of the variability in the timescales is captured by (15) (R 2 values are all greater than 0.75 except in the case of τ FW obs ; Table 2), interactions between tidal amplitude and river discharge are apparent that are not captured by the regression (Fig. 7). For example, for all of the timescales, dependence on tidal amplitude increases slightly as river The artificial skill for this regression model is 0.3 with 95 % confidence and the probability of obtaining a value of R 2 greater than .8 if the null hypothesis that the model has 0 skill was true is 1×10 −6 . An R 2 value of .58 is also very unlikely under the null hypothesis, with a p value probability of 1×10 −3 Fig. 8 The variability of each of the two main parameters in the calculation of τ TEF , the salinityweighted volume Q in (a) and V TEF (b), with tidal amplitude (A T ) and river flow (Q R ) for the 20 idealized simulations. Much of the tidal variability is captured by the exchange flow term Q in , while V TEF represents most of the dependence on river forcing discharge decreases; for τ F , γ ranges from −0.75 to −3.3 from high to low discharge.
Only τ FW has a greater dependence on Q R than on A T , but it also has the least dependence on tidal forcing than all of the timescales. The dependence on A T for most of the runs is almost directly an inverse relationship (γ ∼ −1) except that there is a weaker dependence for the cases when the freshwater fraction method is applied and τ F has a stronger dependence on A T .
A comparison of the bulk timescales to the time series of particles remaining in the estuary for the particle tracking simulations provides another means to assess the forcing dependence of all of the timescales. For example, τ TEF aligns with times when about 20-50 % of the flushing time simulation particles remain within the estuary. τ FW aligns with times when about 15-30 % of the particles remain. An example of this comparison for a subset of the runs is shown in Fig. 5.

Particle Tracking Timescales
There are significant differences between the transit and flushing timescales that were determined from particle tracking. For example, τ T is much greater than τ F , as can be seen by the fitting parameters, α, which are ∼33 and 5 tidal cycles, respectively. For every run, τ T is at least two times as large as τ F and this difference is greatest during low river flow. This is due to the fact that τ F averages particles released over the entire domain (including close to the mouth where time spent in the estuary is small), while τ T is calculated for particles that are released far from the mouth, in the narrow upriver portion of the domain. When the source of particles is far from the region where the particles exit the domain and mixing is incomplete, the local timescale is longer than would be expected from turnover time theory ( Fig. 9; Bolin and Rodhe 1973).
τ T has nearly three times as much dependence on Q R than does τ F (β=−0.84 and −0.31, respectively). This is likely due to differences in transport and dispersion mechanisms in different regions of the estuary (Vallino and Hopkinson 1998). In the narrow, upriver channel, dispersion processes such as shear dispersion and tidal trapping and stirring are inhibited, causing the timescale in this portion of the domain to be dominated by advective transport by Q R . As the estuary widens and become more bathymetrically complex toward the mouth, tidal dispersion mechanisms likely dominate over transport by Q R .
Both τ T and τ F have strong dependences on tides (γ=−0.92 and −1.85, respectively). τ F has nearly two times the tidal dependence as τ T , which we also attribute to the different transport and dispersion mechanisms dominant over different portions of the domain. What is surprising is the very strong dependence of τ F on tides, which exceeds the value for γ predicted from the tidal prism method (γ∼−0.77). This strong scaling with tides warrants further investigation, which the authors recommend as a topic for future study.
The differences in scaling and overall amplitude between τ T and τ F demonstrate the importance of selecting a timescale that is closely associated with the ecological or water quality application being considered. Using particle tracking, there is a benefit to being able to consider a specific release location from within the region that is being flushed, such as in the difference between flushing time and transit time. For example, flushing time might be more descriptive of a phytoplankton bloom that is introduced throughout the estuary at t 0 , while Fig. 9 The spatial variability of flushing time, determined by the average time for all of the particles released from each grid location to exit the estuary, for the simulation with Q R =100 m 3 /s, A T =125 cm (a) and Q R = 5 m 3 /s, A T =125 cm (b) transit time might be more representative of the release of larvae or a pollutant from a specific location within the estuary.

Bulk Timescales
In a comparison of the bulk turnover timescales with the particle tracking timescales, we would expect a closer agreement between the bulk timescales and τ F , which integrates over the volume of the estuary, and we focus on this comparison.
The tidal prism method predicts a timescale that is about one tidal cycle across all forcing conditions, scales roughly as the inverse of tidal forcing, and is very weakly dependent on river flow. This estimate of a timescale is not realistic for the Yaquina Bay, as we would expect based on the assumptions underlying this method (namely that P<<V).
There is a significant difference between τ FW and τ FW obs as can be seen by the significant differences between all of the fit parameters in the regression models for each ( Table 2). The fitting parameter α is nearly three times greater for τ FW , which has nearly five times greater dependence on Q R and half the dependence on tidal amplitude as τ FW obs . The magnitude of these differences highlights the ambiguity in the operational implementation of the relatively straightforward equations for calculating bulk timescales. S AVG is not necessarily well defined, since V can be difficult to define in estuaries. In addition, due to temporal and spatial variability in the salinity structure of estuaries, S AVG can be difficult to measure or accurately estimate in practice. Choosing the tidally averaged salinity from a single location may not be a useful representative measure in an estuary that is not completely mixed.

TEF Timescale
The TEF results were consistent with results found by MacCready (2011) for the Columbia River estuary; S in stays close to oceanic values for all of the simulations and Q in increases with tidal amplitude. The lack of dependence of Q in with Q R was also found in TEF analyses from numerical simulations of both the Merrimack and Columbia Rivers, short salt-wedge-type estuaries where salinity intrusion is close to the length of the tidal excursion (MacCready 2011; Chen et al. 2012). Chen et al. (2012) noted the opposite trend where Q in was greater for neap tides for the Hudson River estuary, a long partially mixed estuary, indicating that the tidal dependence of τ TEF is possibly not consistent across types of estuaries.
Both Q in and V TEF are well constrained and variable under different tidal and river forcing. τ TEF is therefore straightforward to calculate for an estuary from numerical simulations or observations. τ TEF predicts the most similar tidal dependence to τ F as all of the bulk methods, but predicts much smaller river dependence and tidal forcing dependencies as τ F . Based on this study, it underestimated the directly measured timescales (from particle tracking) and, more importantly, did not capture the strong dependence on river discharge predicted by τ F . Overall, this method does not provide any clear advantages or disadvantages over the other bulk methods.
Application to "Real" Estuaries A number of idealized assumptions are implicit in the turnover time theory-i.e., forcing is in a steady state and the time for complete mixing in the estuary is much shorter than the turnover timescale. Under these assumptions, flushing of water to the ocean is purely a diffusive process that can be described by a single timescale. In theory, under these idealized conditions, τ TEF , τ TP , τ FW , τ T , and τ F are all equivalent timescales.
However, estuaries are complex forms of reservoirs where transport and mixing can be spatially heterogeneous and are forcing dependent, the boundaries and volume are often poorly constrained and can vary over time, and steady state is usually an unrealistic condition. Even under the simplified forcing with sinusoidal tides and steady river flow imposed in each of the model simulations used in this study, the timescales analyzed here varied significantly in magnitude and responded differently to river and tidal forcing. If the estuary was completely mixed, all of these timescales could be interpreted as the characteristic e-folding time of an exponential decay function. The difference between the decay of particles remaining in the estuary during the flushing time simulations and the exponential fit (Fig. 5a) means that it is unclear how to interpret each of the bulk timescales and how to interpret all of the timescales in relation to each other.

Temporal Variability in River Flow and Tidal Amplitude
In the Yaquina Bay, low to moderate flows can persist for timescales greater than the turnover times reported here, but the tides are mixed semi-diurnal with a spring-neap cycle that varies at the same timescale or more rapidly than the steady state transit and flushing times. The estuary responds more slowly than the spring-neap cycle and a steady state for a particular tidal forcing may not be reached.
Many studies using various methods across a range of estuaries have found the same pattern that estuarine timescales decrease with increasing river flow (Miller and McPherson 1991;Vallino and Hopkinson 1998;Ferguson et al. 2004;Wang et al. 2004;Huang 2007;Liu et al. 2011;Barcena et al. 2012). During high river flow conditions, all of the methods applied to the Yaquina Bay indicate that the turnover time is on the order of a few tidal cycles. It is worth noting, however, that within the Yaquina Bay a significant fraction of the annual flow occurs in the form of episodic pulses over the same timescale, and flows greater than 100 m 3 /s are generally not sustained for longer than about a day. In this idealized study, Q R was held constant for each simulation and the effect of pulses of freshwater were not specifically considered. We leave a study of turnover time in the Yaquina Bay estuary under realistic, unsteady conditions as a topic for future research.

Defining an Estuarine Volume
In order to calculate a timescale for a reservoir, the boundaries of that reservoir have to be well defined. For estuaries, the boundary and therefore volume are often difficult to define, and all of the methods explored in this paper rely on these definitions. For a specific application, the choice of volume may be obvious based on a specific ecological parameter, but for this theoretical study, there was no such imposed boundary. Instead, we relied on salinity and the physical constriction between the jetties as guides. We performed some analysis to determine the sensitivity of each timescale method to changes in the upriver and/or downriver cross sections. However, the scale by which we changed the estuary volume is small compared to the range of salinity intrusion over a tidal cycle and over the range of river flow conditions. While each timescale may not be sensitive to small perturbations of the definition of estuary volume, this underestimates the amount by which these timescales could vary if volume was defined by an extreme of the salinity intrusion length, for example. The range of salinity intrusion in the Yaquina Bay can mean that the estuary length is almost doubled between high and low river flow conditions. We also see V TEF double over the range of forcing scenarios in Fig. 8b. This sensitivity could also be even greater for other estuaries, such as the Hudson River, where its length defined by salinity intrusion could vary by a factor of 3 or more (Lerczak et al. 2009).
Although b was chosen by the annual mean salinity intrusion over all of the forcing scenarios, it is not clear to the authors that that is an appropriate choice for defining the estuary volume for all of the forcing scenarios. For example, for the freshwater fraction method, one could use a salinity intrusion length that varied depending on forcing. In this case, as Q R increased, V might decrease and τ FW may have been more sensitive to Q R . In order to make a more direct comparison, the flushing time would be calculated with this volume definition and would also likely have a greater dependence on Q R . A comparison to this forcing-dependent definition of estuary volume based on salinity intrusion is beyond the scope of this study.

Spatial Considerations
As in many previous studies of estuarine timescales, the Yaquina Bay has a large amount of spatial variability in local flushing time (Oliveira and Baptista 1997;Vallino and Hopkinson 1998;Shen and Haas 2004;Shen and Lin 2006;Yuan et al. 2007;Warner et al. 2010;Barcena et al. 2012). Therefore, depending on the specific application, a bulk timescale for the entire estuary may not be appropriate.
Near the mouth of the estuary, the local flushing time does not vary as much with river discharge as it does upriver, where the differences in local flushing time between the lower and higher river flow examined could be nearly a month (Fig. 9). This large upriver dependence on Q R has also been found in studies of other estuaries (Oliveira and Baptista 1997;Vallino and Hopkinson 1998;Shen and Haas 2004;Yuan et al. 2007). This indicates that the river forcing may be dominant over tidal forcing upriver, which is supported by the fact that τ T , which was determined from only releasing particles upriver, has greater dependence on Q R than does τ F . Given that τ TEF has very little dependence on Q R compared to τ F (β=−0.14 compared to β=−0.31), it is possible that for the Yaquina Bay, this timescale is more descriptive of the turnover time in the lower portion of the estuary near the mouth.

Conclusion
Particle tracking methods are a good tool for quantifying volume average and local timescales over a range of forcing conditions and for assessing the utility of bulk timescales for particular estuaries and applications. The strong scaling of τ F with tides (γ=−1.85) warrants further investigation. In this comparison of particle tracking timescales to those obtained from bulk formulas, we find deficiencies in the ability of the bulk formulas to capture the amplitude of these timescales and their variations with forcing. Ambiguities in the bulk formula parameters and in the application of these timescales further the need to be careful in their use. This analysis puts into question the utility of these bulk formulas in estuaries such as Yaquina Bay.
τ TEF captures the same general trends with forcing as do other methods of calculating an estuarine timescale, but it shows significantly less dependence on Q R than do the other timescales which indicates that it does not capture seasonal and shorter timescale variations of turnover time in this estuary associated with variations in river discharge. This method, as well as the other bulk methods under investigation, does not appear to be applicable to the Yaquina Bay.