Simple relationships between residence time and annual nutrient retention, export, and loading for estuaries

Simple mathematical models are derived from mass balances for water and transported substance to provide insight into the relationships between import, export, transport, and internal removal for nonconservative substances in an estuary. Extending previous work, our models explicitly include water and substance inputs from the ocean and are expressed in terms of timescales (i.e., mean residence time and the timescale for net removal). Steady‐state, timescale‐based expressions for ratios of export to import, retention to import, and net export to loading, as well as for loading and annually averaged concentration, are provided. The net export:loading model explains the underlying mechanisms for a well‐known empirical relationship between fractional net export and residence time derived by other authors. Although our simplified models are first‐order approximations, the relative importance of physical and biochemical processes influencing export or retention of a substance can be assessed using mean residence time and the timescale for net removal. Assumptions employed in deriving the simplified models (e.g., well‐mixed, dynamic steady state) may not be met for real estuaries. However, model application to Chesapeake Bay for 1985–2012 demonstrates that interannual variations in total nitrogen (TN) net export:loading can be evaluated, and annual nutrient loadings can be well estimated using numerically modeled time‐varying mean residence time, observation‐based mean concentration, freshwater inflow, and an appropriately estimated removal timescale. Our model shows that net fractional export of TN loading ranges from 0.3 to 0.5 over the 28‐yr period. The models can be employed for other substances and water bodies if the underlying assumptions are applicable.

Steady-state, timescale-based expressions for ratios of export to import, retention to import, and net export to loading, as well as for loading and annually averaged concentration, are provided. The net export:loading model explains the underlying mechanisms for a well-known empirical relationship between fractional net export and residence time derived by other authors. Although our simplified models are first-order approximations, the relative importance of physical and biochemical processes influencing export or retention of a substance can be assessed using mean residence time and the timescale for net removal. Assumptions employed in deriving the simplified models (e.g., well-mixed, dynamic steady state) may not be met for real estuaries. However, model application to Chesapeake Bay for 1985-2012 demonstrates that interannual variations in total nitrogen (TN) net export:loading can be evaluated, and annual nutrient loadings can be well estimated using numerically modeled time-varying mean residence time, observation-based mean concentration, freshwater inflow, and an appropriately estimated removal timescale. Our model shows that net fractional export of TN loading ranges from 0.3 to 0.5 over the 28-yr period. The models can be employed for other substances and water bodies if the underlying assumptions are applicable.
Anthropogenic input of excessive nutrients is the main cause of eutrophication over the past few decades (Nixon 1995;Carpenter et al. 1998;Diaz 2001). Both the amount of nutrient input and the nutrient retention time contribute to the formation and severity of eutrophication in a waterbody. Because nutrients are carried by the fluid in a waterbody, retention of nutrients depends, in part, on the transport processes within the waterbody. Time spent during transport within an estuary determines time of exposure of substances to source-sink processes occurring inside the estuary, thus also influencing how much substance is available for export (Lucas et al. 2009). A first-order characterization of transport processes is provided by the water residence time (Monsen et al. 2002). Residence time is a key timescale for estimating the time spent by constituents or pollutants in a water body (Delhez et al. 2004) and is one of the most widely used concepts for quantifying the renewal of the water. Residence time is one of several "transport timescales" (e.g., flushing time, age, turnover time, freshwater replacement time) that can be estimated to distill complex hydrodynamic processes (Takeoka 1984;Monsen et al. 2002).
A number of authors have previously explored relationships between the nutrient balance and residence time or other transport timescales (Lucas and Deleersnijder 2020). For example, Vollenweider (1975) and Dillon and Rigler (1974) used simple mass balance models to show that nutrient retention in lakes is a function, respectively, of "mean residence time" or "replenishment coefficient." Nixon et al. (1996) studied the empirical relationship between "residence time" and nutrients exported from estuaries and lakes. They found that the net export:import ratios for total nitrogen (TN) and total phosphorus (TP) can be described well by empirical linear-log functions of residence time. Dettmann (2001) derived a simplified mass balance model for TN, proposing relationships between freshwater replacement time and annually exported nitrogen and denitrification in estuaries. Using that simplified model, observations, and estimated transport timescales for several estuaries, he showed that the ratio of net TN export to upland loading and the fraction of upland TN loading denitrified can both be estimated based on the net annual removal rate and (surrogates for) freshwater replacement time (also referred to by that author as "freshwater residence time"). It should be noted that the previous mass balance-derived models of Dettmann (2001) and Vollenweider (1975) did not explicitly account for or delineate inputs from a downstream adjacent water body (in the case of an estuary, the ocean). Moreover, the definitions of "residence time" implemented by the above authors are not all consistent.
Although these previous models advanced our understanding of how residence time influences export and retention of reactive substances in aquatic systems, some questions still need to be addressed: (1) Can the simple mass balance approach for predicting export ratios be extended to explicitly account for material import from the ocean, and how would such a new model differ from previous ones? (2) Are the empirical and extended mass balance approaches for predicting export ratios consistent with each other? (3) Are there general timescale-based expressions relating constituent loading, ocean input, internal loss, mean concentration, and export that apply not only to nutrients, but potentially also to other water quality constituents for estuaries or other waterbodies (within assumptions)?
One of the challenges in applying these simple models is the determination of residence time. It is frequently unclear what kind of "residence time" or other transport timescale should be used, as the values associated with different transport timescales and even with different definitions for "residence time" can vary substantially (Monsen et al. 2002;Lucas and Deleersnijder 2020). A large error may be expected if an inappropriate timescale is used. In addition, residence time can vary significantly with freshwater discharge and other dynamic conditions in an estuary, so usually there is not a single value applicable for all time. A clear definition and an accurate estimation of residence time is therefore needed.
For this study, we adapted the concepts and approaches of Vollenweider (1975) and Dettmann (2001) to derive our simplified models. Like their previous models, our new model for exported (or retained) substance depends on both the net removal timescale and the residence time. One difference is that the simple models herein are based on two mass balances (water, substance); whereas those authors built their models using only conservation of substance mass. Because the conservation of water and substance are both incorporated into our models, the contribution of import at the seaward boundary can be explicitly delineated and accounted for (a key difference relative to the previous models).
The ultimate objective of this study is to provide insight into the relationship between the transport and net removal timescales and the mass balance of transported, reactive material in an estuary, to better understand the dynamics of material retention and export. Application of our models to the Chesapeake Bay demonstrates that they are useful for estimating annual TN net export:import and total annual loadings for TN and TP. The models can also be applied to lakes and other waterbodies for other substances if the underlying assumptions are applicable.

The method
The goal in this section is to derive basic relationships relating fluxes, loads, and concentrations of substances to transport and loss processes in an estuary, and to express these relationships in terms of timescales. (Timescales are useful encapsulators of complexity, sometimes leading to simple but useful mathematical models; Lucas and Deleersnijder 2020.) We derive these relationships with the ultimate objective of describing the relevant balances and processes as averages over the estuary and over a chosen period of time (e.g., 1 yr), while acknowledging that many processes contributing to distributions and transport of materials in an estuary are inherently spatially and temporally variable over finer scales.

Balance of water and mass
As the foundation for our derivations, we begin with mass balances for water volume, and then for mass of the transported constituent. We assume the estuary or coastal embayment of interest is a partially enclosed waterbody with tidal inflow and outflow through one principal opening. The timeaveraged mass balance of water volume over a specified averaging period T (e.g., T = tidal cycle, day, month, or year) can be written as follows: where each variable is an average over the averaging period T but may vary between periods. A detailed derivation is presented in Supporting Information Section S1. (For simplicity, we have replaced variables with over-bars [used in Supporting Information Section S1 to denote time-averages] with upper case variables here. Equation 1 [Eq. 2] is therefore equivalent to Supporting Information Eq. S5 [Eq. S7].) Q in (m 3 d À1 ) is the time-averaged total (mixed new ocean and ebb) water that enters the estuary through the ocean boundary; Q e (m 3 d À1 ) is the time-averaged total mixed water that leaves the estuary; Q f (m 3 d À1 ) is the average total daily freshwater input; and V is the mean volume of the waterbody (m 3 ).
For material transport, we assume that the time-averaged total mass of the substance in the estuary over the averaging period T can be expressed by a spatial and temporal mean concentration C (g m À3 ) times the mean estuary volume. We further assume that (1) the mass of the substance exported due to transport is proportional to this mean concentration, and (2) the internal net removal (or regeneration) of the substance from the water column due to biochemical processes or settling is also proportional to the mean concentration (i.e., a first-order reaction). The first-order time-averaged mass balance for a transported, reactive substance can then be written as follows (see Supporting Information Section S1 for derivation): where each variable is an average over the averaging period T. L is time-averaged total mass loading discharged into the waterbody (g d À1 ). The coefficient k (d À1 ) is a simplified representation of the net removal rate that accounts for biochemical processes, settling, and other processes (e.g., bottom flux). For example, the removal of TN includes both settling to the bottom and denitrification. For the cases considered herein, the net removal rate is positive; however, Eq. 2 may be applicable for situations where that rate is negative (generation or growth). In the time-averaged substance mass balance equation (Eq. 2), the net removal rate is the mean value over the averaging period T. Q e C is an estimate of the time-averaged mass exported due to transport at the ocean boundary. Expression of the outward flux as Q e C indicates an implicit assumption that the estuary is well mixed (i.e., that the outgoing concentration is the same as the mean internal concentration; see discussion of this later in the "Assumptions and Caveats" section). C in is the mean concentration of the substance entering at the seaward boundary and represents the combination of mixed new substance from the ocean and ebb water over the period of concern.
To further simplify our derivations and understand the relationship between mass transport and transport timescales, we next express the water and constituent mass balances in steady-state form. Under the steady-state assumption, Eqs. 1 and 2, respectively, can be written as follows: Although the assumptions employed to derive Eqs. 3 and 4 neglect spatial variation within the estuary and temporal variation within the averaging period, it will be shown that these expressions can be applied to study the mean condition and variability between averaging periods of length T. See "Assumptions and Caveats" for further discussion of these assumptions.

Mean residence time
Both the average inflow from the sea (Q in ) and outflow to the sea (Q e ) in Eqs. 3 and 4 are often unknown. In this section, we will show that they can be estimated using the mean residence time for an estuary. Following many previous authors, residence time has been defined as the time taken by a particle of a constituent "to leave a water body or defined region of interest" (Lucas and Deleersnijder 2020). Here, we follow the definition of "particle" provided by Lucas and Deleersnijder (2020), that is, a material point possessing zero volume but nonzero mass of a single transported constituent, thus containing potentially many like molecules, cells, grains, and so on that share the same history. By the above definition, residence time varies spatially and temporally, as particles located in different parts of the domain, or discharged at different times or tidal phases, will likely take different amounts of time to leave. A "strict" definition of residence time employed by some authors (Monsen et al. 2002;Delhez and Deleersnijder 2006;De Brauwere et al. 2011) is the time taken by a particle to exit the domain for the first time; this is a critical distinction since, in tidal environments, a particle may leave on one tidal phase and return (and exit again) on subsequent tidal phases (Lucas and Deleersnijder 2020). The strict definition is employed herein.
In the following, we define "τ r " as the spatial mean residence time; the "mean residence time" is also temporally averaged. For a well-mixed estuary, the seaward outflow (Q e ) can be estimated using the residence time (see derivation in Supporting Information Section S1, "Relation between residence time and outflow"). If we assume that the averaging period is long enough to approximate infinity relative to the timescale for flushing of the estuary, then the seaward outflow can be estimated by the following equation: where all three variables are assumed constant over the averaging period and, in practice herein, are represented as averages over that period. The exceedance of the mean residence time by the averaging period appears to represent a minimum requirement for applicability of Eq. 5. Equation 5 here is the same as Supporting Information Eq. S13b. The mean seaward outflow (Q e ) can thus be derived from Eq. 5, the space-and time-averaged residence time, and mean water volume. Note that the form of Eq. 5 is analogous to other authors' definitions for "mean residence time" (Vollenweider 1975), "flushing time" (Monsen et al. 2002), or "turnover time" (Sheldon and Alber 2006), though with flow rates and/or volumes defined differently in some cases. Equation 5 is assumed to hold within an averaging period; however, its component parameters may vary between averaging periods. Once we know the mean seaward outflow, the mean inflow from the ocean (Q in ) can be obtained from Eq. 3 under the hydrodynamic steadystate condition, assuming freshwater flow rate Q f is available.

Exported fraction of material
Previous authors have explored empirically (Nixon et al. 1996) or mass-balance derived (Dettmann 2001) relationships between nutrient export, import, and hydrodynamic transport processes in estuaries. Next, we derive new simple mathematical relationships that account for physical processes previously neglected and that may be employed for substances other than nutrients if relevant assumptions apply. By using Eqs. 4 and 5, the export to import ratio (E:I) for an estuary (i. e., mass exported normalized by the mass imported) can be estimated as follows under the steady-state condition: Here we define τ k , the timescale for net removal of the substance, as the reciprocal of the net removal rate (k). By converting internal removal and transport to timescales, those different processes can be compared with a single comparable currency (Lucas et al. 2009;Lucas and Deleersnijder 2020). For example, if transport is faster than internal loss, then the residence time is less than the net removal timescale, and the ratio of residence time to the net removal timescale is less than unity. If loss is faster than transport, then the net removal timescale is less than the residence time and the ratio of residence time to the net removal timescale is greater than unity. The export:import ratio depends on the relative rates of internal removal and transport and is thus determined by the dimensionless ratio of residence time to the net removal timescale (τ r /τ k ; see Fig. 1a, which shows the relationship for positive net removal rate only). It can be seen in Fig. 1a that as transport becomes very fast compared to net removal and therefore the ratio residence time:net removal timescale approaches zero, export:import tends toward unity (transport is so fast that not very much removal occurs within the estuary and most imported material is exported to sea). On the other hand, as transport becomes very slow compared to net removal and therefore the ratio residence time:removal timescale tends toward infinity, export:import tends toward zero (transport is so slow and time spent within the estuary is so long that most material is removed within the estuary and never exported to sea). The general behavior of the relationship is consistent with the findings of Lucas et al. (2009), who showed that for a system that functions as a net sink for a transported substance, the longer the transport time (i.e., the further to the right in Fig. 1a), the lower the ratio of exported mass to imported mass. This concept has also been demonstrated in the context of tributary nutrient loads to lakes, with higher river flows (shorter transport times) producing higher concentrations in receiving waters because "higher flows moved phosphorus further, faster, with less attenuation" (Stow et al. 2019). Once export:import is known, the retention to import ratio for an estuary can be estimated as retention: import = 1export:import. It is interesting to note that retention: import thus calculated for estuaries based on Eq. 6 reduces to the time-scale based expression derived by Vollenweider (1975) for lakes (see his eq. 2.14); our relationship differs, however, in that "import" herein explicitly includes inputs from the ocean. Nixon et al. (1996) expressed exported nutrient (TN or TP) as the ratio of the net exported TN or TP ("NE") to total loading from land and the atmosphere (what we call L). Employing Eqs. 4 and 5 and our timescale definitions, this ratio can be expressed as: , which we call the "ocean exchange factor." Form (d) of this equation is the same as Eq. 6d, except that here the net removal rate k is modified by the ocean exchange factor β, which accounts for material input to the estuary through the ocean boundary. Eq. 7 is also very similar to Dettmann's (2001) eq. 9, and it would be the same (but for different transport timescales employed) if we assumed Q in C in < < Q e C, that is, substance mass input from the ocean is negligible relative to gross export (ocean exchange factor is approximately 1). Such an assumption would be consistent with that author's assumption of negligible influence of the sea on the nitrogen mass balance in an estuary. However, the ocean exchange factor β should not always be assumed to be 1. For example, estimates of the ocean exchange factor for TN in the Chesapeake Bay can vary from 1.2 to 2.2 with an average of about 1.5, based on a combination of model outputs (Du and Shen 2017) and measured TN concentrations by the Chesapeake Bay Program (https://www.chesapeakebay.net/ data). In that case, if the ocean exchange factor were incorrectly assumed to be 1 and if Eq. 7 were used to estimate net removal rate k, one could overestimate net removal rate by 20-120% because a portion of the reduction in net export:loading would be incorrectly attributed to net removal, instead of to mass input from the ocean. If we define the "adjusted net removal rate" K as the product of the ocean exchange factor and the net removal rate (i.e., K = βk), then the distributions of export:import (Eq. 6d,e) and net export:loading (Eq. 7d,e) are the same (see Fig. 1a), except in the latter case the timescale for net removal τ k should be replaced by the timescale for adjusted net removal τ K (i.e., 1/K). Equations 6 and 7 can be applied to any pollutant if the assumption of a first-order removal rate and other assumptions employed above are applicable.
Equation 7 indicates that net export:loading will be negative if material input at the ocean boundary is larger than outputs from the estuary and positive if gross material export to the ocean exceeds oceanic inputs. Equation 7e is depicted in Fig. 1b, which shows behavior of the relationship for positive and negative net removal rate and net export. Regime "1" (blue shaded area) represents those cases for which net removal rate is positive and 0 < net export:loading < 1. For this regime of net removal, it makes sense that any positive net export:loading ratio is less than unity, since internal removal diminishes mass loaded to the estuary before it is exported ( Fig. 1a). Regime "2" (green shaded area) represents cases for which net removal rate is negative (net generation or growth within the estuary, per definitions employed herein); Regime "2" accordingly reflects cases for which net export:loading is greater than unity, since mass loadings are amplified by net generation or growth before export from the estuary. The ocean exchange factor β is positive for both Regimes "1" and "2," signifying positive net export:loading. Regime "3" (red shaded area) represents the case of positive net removal and negative net export:loading (negative ocean exchange factor); physically, this case reflects positive net input from the ocean which, along with other (upland, atmospheric) loadings, must be balanced in steady state by net removal within the estuary (positive net removal rate). This regime corresponds to the case of TP in the Chesapeake Bay during 1985-1986(Boynton et al. 1995. Note that the net export:loading form in Eq. 7e does not work for the special case Q in C in ¼ Q e C for which net export is zero (since the ocean exchange factor β would equal infinity), or for the case where β τr τ k ¼ À1 (for which net export: loading would equal infinity). In theory, Eq. 7 works for all three regimes depicted in Fig. 1b; however, we have only applied and tested it herein for Regime "1" (positive net removal rate, 0 < net export:loading < 1).

Annual loading and mean concentration
Using Eqs. 4 and 5, the first-order approximation for total loading of a substance to an estuary can be obtained as follows: where τ f = V/Q f is the timescale characterizing replacement of the estuary volume by freshwater only (i.e., the "freshwater calculated using Eq. 6. (b) The ratio of net exported matter to loading (NE:L) as a function of the ratio of residence time τ r to the timescale for adjusted internal removal τ K ; calculated using Eq. 7. The blue shaded area (Regime "1") represents those cases for which net removal rate k is positive (net removal) and 0 < NE:L < 1. The green shaded area (Regime "2") represents cases for which k is negative (net growth within the estuary) and NE:L > 1. The red shaded area (Regime "3") represents the case of positive k (net removal) and negative NE:L (net import through the ocean boundary).
(c) The ratio of mean estuarine concentration to the annual "maximum loading concentration"; calculated using Eq. 10.
renewal timescale"). Note that terms (a) and (b) (i.e., those which capture the loading, transport within the water body, internal removal, and mean estuarine concentration) rearrange to Vollenweider's (1975) eq. 2.12 for lakes; however, our term (c) incorporates the added effect of input from the ocean. Interpretation of Eq. 8 is aided by nondimensionalization, which produces a simpler dimensionless form of the relation: where we define the dimensionless loading as L * = L Â τ r / (VÁC in ), the dimensionless removal timescale as τ Ã k = τ k /τ r , the dimensionless freshwater renewal timescale as τ Ã f = τ f /τ r , and the dimensionless mean estuarine concentration as C * = C/C in . Casting the equation in this way reduces a sevendimensional problem (one dependent variable as a function of six parameters, Eq. 8) to a four-dimensional problem (one dependent variable as a function of only three parameters, Eq. 9) and makes it easier to understand how the main processes interact. For example, the (dimensionless) mean concentration of the substance within the estuary depends linearly on the (dimensionless) loading, modulated by the (dimensionless) removal timescale. Moreover, large (dimensionless) freshwater flow (i.e., small [dimensionless] freshwater renewal timescale) decreases (dimensionless) concentration via dilution.
As a means of developing a simple estimate of mean concentration within an estuary-and understanding its dependence on residence time and the net removal timescale-let us first define the "maximum loading concentration" (C max ) as the concentration that would be attained for a substance if total inputs (from the watershed, atmosphere, and sea) over a specified period were added to an estuary volume V without internal removal or export (a hypothetical "lossless" estuary). The realized concentration C (with internal removal and export), expressed as a ratio with maximum loading concentration, can be described by: where C max ¼ Q in CinþL V T load and T load is the time over which inputs continue and accumulate within the hypothetical lossless estuary (i.e., the "loading period"). The timescale-based term (c) above is derived by substituting Eqs. 4 and 5 into Eq. 10b. If the loading period is 365 d, Eq. 10 describes the effects of an annual loading. The distribution of Eq. 10 is shown in Fig. 1c for a range of residence time and positive net removal timescale. The ratio of realized concentration to maximum loading concentration increases toward 1 as residence time and net removal time increase (i.e., as transport and net removal rate become slower, and the "real" estuary tends toward the hypothetical lossless estuary). We note that terms (a) and (c) above would provide a relationship equivalent to Vollenweider's (1975) eq. 2.11, if oceanic material input were incorporated into that author's loading parameter.

Computation of residence time
One of the key parameters in the above formulations is the average residence time τ r , which itself is a function of time. It is often computed numerically by injecting particles at a fixed time, following their paths, and recording the times when they leave the embayment for the first time (Monsen et al. 2002;Gong et al. 2008;Blaise et al. 2010). Another method for calculating residence time is to employ Takeoka's (1984) remnant function concept, integrating a model-calculated tracer concentration time-series over time (as in Supporting Information Eq. S10; Shen and Haas 2004;Wang et al. 2004;Wang and Yang 2015). Average residence time computed by each of these approaches depends on the particle (or tracer) release time, making it difficult to obtain long-term time-varying residence time for a large estuary because multiple simulations with particle or tracer release for each day or tidal phase would be required. As an alternative, Delhez et al. (2004) proposed an adjoint method to compute space-and time-varying residence time with a single backward model run. The approach of Delhez et al. (2004) was implemented by Du and Shen (2016) in a calibrated 3D hydrodynamic model of the Chesapeake Bay based on the Environmental Fluid Dynamics Code (EFDC, Hamrick 1992). The adjoint model was integrated backward from 2014 to 1980. The first 2 yr (2014 and 2013) were used for model spin-up. For a detailed description of the model set-up, readers are referred to Du and Shen (2016). The computed time series of annual mean residence time (averaged over the Chesapeake Bay) is shown in Fig. 2 and Supporting Information Table S1. Computed annual mean residence time for Chesapeake Bay ranges from 127 to 245 d, with a mean of 179 d.

Data sources
Two primary sets of data were used for this study. The first was published by Nixon et al. (1996) and used by Dettmann (2001). Nixon et al. (1996) compiled estimates of "residence time," TN and TP loading, net export, and loss for 11 estuaries and eight lakes; one of the estuaries included in this data set is the Chesapeake Bay, based on data from Boynton et al. (1995). Data used herein were obtained from a combination of values provided explicitly by Nixon et al. (1996), Dettmann (2001), or by sources referenced in those papers, as well as from values digitized from Nixon's figs. 1 and 2. Nixon et al. (1996) used these data to study the relationship between residence time and net export:loading. Their values of "residence time" were collected from different published sources and based on a variety of estimation methods; it is important to note that freshwater replacement time and e-folding flushing time (surrogates for residence time used in Nixon's compiled dataset) differ, potentially significantly, from the definition of residence time employed herein. These data were also used and supplemented by Dettmann (2001) to study annual export of TN and denitrification. For details regarding these data, readers are referred to Nixon et al. (1996) and Dettmann (2001) and the original publications cited therein.
The second primary dataset used is that collected by the Chesapeake Bay Program, which has conducted routine observations of water quality parameters including TN and TP in the Chesapeake Bay since 1984 (https://www.chesapeakebay. net/data). Chesapeake Bay Program data are collected bimonthly or monthly throughout the Bay. We used the Chesapeake Bay Program data at stations shown in Fig. 3 to compute volume-weighted annual mean concentrations of TN and TP in the Chesapeake mainstem from 1985 to 2012. We used the annually averaged bottom concentration of three stations (CB7.4, CB7.4N, and CB8.1E) to represent the TN and TP concentration associated with inflow (C in ) from the ocean. The annual mean TN concentrations at Stations CB7.4 and CB7.4N near the Bay mouth have almost the same value, which is slightly less than the annual mean value of TN at CB8.1E at the mouth.
In addition, the Chesapeake Bay Program estimated annual loadings of TN and TP from 1990 to 2012 based on the sum of watershed model results, point sources, and atmospheric deposition (https://www.chesapeakeprogress.com/clean-water/ water-quality). These data were used for verifying our simple model estimates of TN and TP loading and are listed in Supporting Information Tables S2 and S3 in Section S3.
A detailed procedure for using these data, estimating parameters, and implementing our simple models is summarized in Supporting Information S2. The values for each input variable including annual mean residence time, volume-weighted average concentration (C) and concentration at the mouth (C in ) for TN and TP, freshwater discharge Q f (http:// waterdata.usgs.gov/nwis/), average inflow from the sea Q in , and outflow to the sea Q e are listed in Supporting Information Table S1 in Section S3. Definitions, values, and sources for parameters used in our simple models are summarized in Table 1.

Results and discussion
Verification of NE:L model and estimation of removal rate Nixon et al.'s (1996) published nutrient net flux, loading, and residence time data were used for verifying our simple net export:loading model (Eq. 7). Dettmann (2001) used nonlinear regression to fit his own derived model to an adapted version of Nixon et al.'s (1996) TN dataset (estuaries only), estimating the net removal rate "k" as the fitting parameter. Because our Eq. 7 takes the same general mathematical form as Dettmann's but incorporates different assumptions and parameter definitions (including our explicit inclusion of ocean input), the unknown parameter of Eq. 7 is the adjusted net removal rate "K," not the net removal rate k. Using nonlinear regression to fit our own Eq. 7 to Nixon's estuary and lake TN data, we estimated the cross-system adjusted net removal rate K XÀsys TN to be 0.010 d À1 -the same value as Dettmann's (2001) net removal rate k TN .
For TP, we chose to be consistent with the approach of Nixon et al. (1996), thus excluding from our analysis estuaries with negative net export (Chesapeake Bay and Potomac River). Otherwise following the approach for TN described above, we found the estimated cross-system adjusted net removal rate K XÀsys TP to be 0.006 d À1 . We also compared goodness of fit for our simple model to that for Nixon et al.'s (1996) empirical Table 1. Parameters used in simple models and calculations presented herein. "CBP" refers to the Chesapeake Bay Program. "USGS" refers to the U.S. Geological Survey. "TN" is total nitrogen, and "TP" is total phosphorus.

Parameter
Value Description V (m 3 ) 7.5 Â 10 10 Mean estuary volume (Du and Shen 2016) 0.010 Cross-system value for TN adjusted net removal rate. Based on fit of Eq. 7 to Nixon et al. (1996) dataset K XÀsys TP (d À1 ) 0.006 Cross-system value for TP adjusted net removal rate. Based on fit of Eq. 7 to Nixon et al. (1996) dataset, positive net export systems only K Ches TN (d À1 ) 0.01 Chesapeake-specific TN adjusted net removal rate. Based on Eq. 7 herein and τ r and net export:loading from Nixon et al. (1996) and Boynton et al. (1995) Ocean exchange factor, which modifies k to account for material input to the estuary through the ocean boundary.  Nixon et al. (1996) 0.011 Chesapeake-specific net removal rate for TP. Back-estimated from Eq. 10 using estimated 1985-1986 mass input and output data (Boynton et al. 1995), CBP TP observations (for C and C in ), 3D model computed τ r (Du and Shen 2016), and estimated τ f based on 1985 USGS flow measurements γ (d À1 ) 0.0005 Factor used in expression for k Ches TP as a function of total suspended solids (TSS), k TP = γe (ε•TSS) . Value derived by fitting modified Eq. 10 (i.e., with above power function expression for k) to CBP loadings and TSS. Time-average of daily total mixed water that leaves the estuary through the ocean boundary over the averaging period T (1 yr herein). Based on Eq. 5, V, and computed τ r Q in (m 3 d À1 ) Time-average of daily total (mixed new ocean and ebb) water that enters the estuary through the ocean boundary over the averaging period T (1 yr herein). Based on Eq. 3, Q e and Q f Q f (m 3 d À1 ) Time-average of daily total freshwater input over the averaging period T (1 yr herein). Based on USGS measurements C (g m À3 ) Volume-weighted average substance concentration over the estuary and over the averaging period T (1 yr herein). Based on CBP measurements in the Chesapeake Bay mainstem Mean substance concentration at the seaward boundary over the averaging period T (1 yr herein). Based on CBP measurements near the mouth of Chesapeake Bay L (g d À1 ) Time-average of daily mass loading of transported substance from all nonoceanic sources over the averaging period T (1 yr herein) TSS (g m À3 ) Annually averaged total suspended solids concentration (arithmetic average of all stations) in upper and middle bay region based on CBP measurements regression of net export:loading against log (residence time), which we recreate herein. Graphical and statistical comparisons of our model (Eq. 7) to the Nixon empirical regression model are shown for TN and TP in Fig. 4a,b, respectively. Our model fits Nixon et al.'s (1996) data well and is characterized by model skill (measured by correlation coefficient R) nearly identical to that of Nixon et al.'s (1996) empirical models. Both of these comparisons (fitting parameter, goodness of fit) help to verify that our physically derived model and its application to Nixon et al.'s (1996) dataset aligns with the performance and results of previously published simple models.
Interannually varying NE:L for Chesapeake Bay As explained earlier, for a system influenced by nutrient inputs at the ocean boundary, the net removal rate (k) will be lower than the adjusted net removal rate (K = βk), assuming the ocean exchange factor β and net removal rate are both positive. Using Eq. 7d along with the 1985-1986 Chesapeake Bay net export:loading ratio for TN (0.3) and residence time (228 d) used by Nixon et al. (1996), we estimated a Chesapeake-specific adjusted net removal rate to be K Ches TN = 0.01 d À1 . A mean ocean exchange factor of 1.5 was estimated for the Chesapeake (see Supporting Information Section S2, Step 3c) and then was used to back out a (presumed) timeinvariant Chesapeake-specific value for net removal rate, k Ches TN = 0.0067 d À1 . Using Eq. 7d, interannually varying mean residence time and ocean exchange factor, and constant net removal rate k Ches TN , a time series of annual net export:loading was estimated for Chesapeake Bay TN between 1985 and 2012 (Fig. 5). Estimated annual net export:loading ranges from about 0.3 to 0.5, indicating that on average more than half of TN loading to the Chesapeake is taken up or lost during transit within the Bay, and consequently less than half of that loading experiences net flux to the sea. Note that the years with highest estimated net export:loading (1996,2003,2004,2011) correspond to years with particularly short residence times (Fig. 2), consistent with the notion that faster transport through a system (shorter mean residence time) corresponds with less time for a reactive transported constituent to be modified by internal processes (in this case, net loss), thus resulting in a greater portion of imported mass being exported (Lucas et al. 2009).

Interannually varying loading of TN
We used the Chesapeake Bay TN balance as a case study for implementing the simple model in Eq. 8 to calculate a firstorder estimate of annual constituent loading. Sensitivity tests for the loading model (Eq. 8) were conducted by varying the mean value for each independent parameter by AE 20%; sensitivity results are shown in Table 2. The parameter to which the loading model is most sensitive (at least for the range of parameters used herein to characterize the Chesapeake Bay) is mean concentration C; in that case, the percent change in loading is on the order of the percent change in mean concentration, a variable for which reliable measurements are often available for estuaries. The model is also sensitive to estuary volume V (a parameter also commonly available), as loading is linearly related to volume. The variable to which the model is next-most sensitive is the net removal timescale τ k (1/k); we expect this model parameter to be the most challenging to estimate for many estuaries. A AE 20% change in the net removal timescale resulted in a change in loading of +15.9% to À10.6% (Table 2). The relatively low sensitivity of loading to mean residence time is due to the fact that the terms C/τ r Fig. 4. Ratio of net export to loading (NE : L) vs. residence time for (a) TN and (b) TP. Red lines represent the best fit of Eq. 7 to the Nixon et al. (1996) estuary and lake data. Blue lines represent the Nixon et al. (1996) empirical linear-log regression models. Data points, including estuaries (red dots) and lakes (yellow triangles), are from Nixon et al. (1996), Dettmann (2001), and sources referenced therein. The reader is referred to those publications for details. and ÀC in /τ r can approximately cancel each other if the values of the estuary-averaged concentration C and the inflowing concentration at the estuary-ocean boundary C in are close (see Eq. 8).
Using an estimated Chesapeake-specific net removal rate k Ches TN (Table 1), annual mean Chesapeake Bay Program TN concentrations for the estuary and mouth, annually varying computed mean residence time τ r (Fig. 2), and U.S. Geological Survey (USGS) observed flow Q f (to obtain τ f ), TN loading (black circles) for each year from 1985 to 2012 was estimated ( Fig. 6; see Supporting Information Table S1 in Section S3 for input parameters.) Given the demonstrated sensitivity of the loading model to the net removal timescale τ k (1/k; Table 2) and our expectation that the net removal rate k may be difficult to reliably parameterize for many estuaries, we also use Fig. 6 (gray bars) to visually convey how estimated loading responds to a modest estimated variation in net removal rate k (see Step 4a in Supporting Information Section S2, Table 1).
Our estimated mean annual TN loading for the entire Chesapeake Bay is 175 Â 10 6 kg yr À1 , averaged over the period 1985-2012. The interannual variation during this period is significant (Fig. 6), with our estimated maximum loading approximately doubling the minimum loading. For the period 1990-2012, we compare the TN loading estimates from our simple model to those computed by the Chesapeake Bay Program (see Fig. 6). Our simple model predictions (black circles) agree well with the Chesapeake Bay Program's estimates (red triangles): the two estimates are significantly linearly correlated (coefficient of determination R 2 = 0.85), and the rootmean-square difference (RMSD) between estimates is 29 Â 10 6 kg yr À1 (17% of our estimated mean loading). Although the coefficient of determination is high, our model biases high compared to the Chesapeake Bay Program estimate. Our mean loading is about 13% larger than that estimated by the Chesapeake Bay Program. It is important to note that the loading estimated by Chesapeake Bay Program used a weighted regression method based on USGS monitoring data for river input (Moyer and Blomquist 2018) and a Hydrological Simulation Program-Fortran watershed model (Phase 5.3.2) for downstream nonpoint sources, which has since been updated to the Phase 6 version, and therefore should also not be taken as an exact representation of loads. Our loading estimation is sensitive to the specified net removal rate (see gray bars conveying loading range based on the estimated range for net removal rate k Ches TN ; see Table 1); this sensitivity underscores the importance of developing reliable estimates for the net removal rate for use with our simple model. Overall, our model estimates of TN loading are less variable over time than the Chesapeake Bay Program estimates. The main cause of this is likely our use of a constant net removal rate across all years. The internal net loss of nitrogen strongly depends on primary production, respiration, nutrient recycling, groundwater inputs, and other processes which add or remove TN to/from the water column within an estuary.  One specific process that may be contributing to the lower variability in our predicted TN loading is enhanced sediment ammonium release under higher salinity/lower flow conditions (Weston et al. 2010). This mechanism (and our neglect of it as a potential driver of interannual variability in the net removal rate k) could help explain why our model tends to overestimate TN loading (relative to Chesapeake Bay Program estimates) in low-flow years (see Supporting Information Table S1 in Section S3 for flow data). Our neglect of the interannual variation term d(VC)/dt under the steady-state assumption could also affect the results for some years.

Interannually varying loading of TP
We used estimated k Ches TP ( Table 1, Step 5a in Supporting Information Section S2) and Eq. 8 to estimate TP loading to the Chesapeake for the period 1985-2012 (Fig. 7, blue stars). We compared our predicted loading to Chesapeake Bay Program estimates for 1990-2012 (Fig. 7, red triangles). Correlation of our estimate with the Chesapeake Bay Program estimate for TP (coefficient of determination R 2 = 0.52) is much weaker than for TN. For example, estimates based on our simple model exceed Chesapeake Bay Program estimates by about 149% in the year 2000 but are 49% lower than Chesapeake Bay Program estimates in 2011.
We hypothesized that one contributor to this discrepancy in TP loading estimates is the use of a constant Chesapeakespecific net removal rate (k Ches TP ) with our model. Because the settling component of phosphorus removal depends strongly on total suspended solids (TSS) concentration due to sorption/ desorption processes (Chapra 1997;Cerco and Noel 2002), we sought to improve our model by assuming that the Chesapeake-specific net removal rate k Ches TP can be described as a power function of TSS (mg L À1 ), that is, k Ches TP = γe (ε•TSS) . Direct observations of settling loss are unavailable to estimate these parameters; therefore, to explore whether specifying k Ches TP as a function of TSS could improve our model, we fit our modified version of Eq. 10 to 12 yr (2001-2012) of annual Chesapeake Bay Program loadings and observed Chesapeake Bay Program TSS to estimate the model parameters γ and ε.
(We used mean annual observed TSS concentrations measured in the upper and middle Bay to represent the TSS concentration, as a large amount of TP will settle in the turbidity maximum zone in the upper Bay region. See Supporting Information Table S1 in Section S3.) The estimated TSS fitting parameter values are γ = 0.0005 and ε = 0.1458 (Fig. 8). For TSS ranging from 10 to 25 mg L À1 , the estimated Chesapeakespecific TP net removal rate k Ches TP ranges from 0.0021 to 0.0191 d À1 . We applied the removal rate-TSS equation to compute the TP removal rate and an improved estimate of TP loading for the entire 1985-2012 period. The years 1990-2000 are treated as validation. A comparison of TP loading estimated by our improved model (black circle symbols) and by the Chesapeake Bay Program is shown in Fig. 7. Correspondence between our model and the Chesapeake Bay Program estimate is improved overall with incorporation of the removal rate-TSS relationship (RMSD is 3.7 Â 10 6 kg yr À1 , coefficient of determination R 2 = 0.62) for the period of 1990-2012, but a few abnormal years were characterized by significant mismatch (e.g., 1996, 1998, 2005, 2011). Discrepancies between our simple model and Chesapeake Bay Program loading estimates could be associated with our neglect of intra-annual variability in the net removal rate and/or other biogeochemical processes influencing the net removal rate but not captured by the TSS-dependent function used herein. Overall, this Bay. The blue curve shows the loading estimation based on Eq. 8, data from 1985 to 2012, and a constant net TP removal rate of k Ches TP = 0.011 d À1 . The black curve shows the loading estimation using a modified Eq. 8 that includes a net removal rate expressed as a dynamic function of TSS concentration. The red curve is the loading estimated by the Chesapeake Bay Program (CBP), that is, the sum of watershed loading, point sources, and atmospheric deposition. The coefficient of determination R 2 values reflects agreement between our simple model-based loading estimates and the CBP estimate. exercise suggests that our simple model estimate of TP loading can be improved if we know the time-varying net removal rate (or an empirical relationship with a surrogate measure such as TSS).
One other possible explanation for mismatch between our modeled loadings and the Chesapeake Bay Program's estimates is our implicit assumption that the mean estuary nutrient concentration responds essentially instantaneously to loading. In reality, given the long residence times associated with the Chesapeake Bay, the mean estuary concentration could take 100 d or more to reflect a loading event (Du and Shen 2016). Therefore, a loading event late in the calendar year could influence real concentrations in the next calendar year, while our simple loading model using those elevated concentrations would associate the loading with the latter year. Such calendar year cut-off effects could lead to over-or underestimation of loading by our model.
The TP loading estimation is evidently very sensitive to the net removal rate. Because the TP net removal rate is highly variable (Chapra 1997), large uncertainty can be expected for the loading estimation because the relationship between removal rate and TSS is often unknown. More studies are thus warranted for estimating net removal rate for TP in estuaries and lakes.

Assumptions and caveats
Three major assumptions employed in deriving the simple models herein are: (1) the mass flux of a substance out of an estuary is proportional to the mean concentration of the substance within the estuary; (2) the internal net rate of removal of a substance within the estuary is proportional to the mean concentration of the substance in the estuary; and (3) the system is in a dynamic steady state. Let us step through these assumptions one by one.
The first translates into the assumption that the concentration of substance exiting the system at the ocean boundary is the same as the average concentration within the system (i.e., the "well-mixed" assumption). This assumption was implemented in Eq. 2 and was needed to translate Takeoka's (1984) expression for averaged residence time (our Supporting Information Eq. S10) into the simple relationship in Eq. 5. That expression for mean residence time, which is analogous in form to other authors' definitions of flushing time (Monsen et al. 2002;Delhez et al. 2004) or turn over time (Sheldon and Alber 2006), was further used in developing the simple timescale-based models herein. Depending on the constituent and the estuary, the well-mixed assumption is frequently violated in reality because substances are often not uniformly distributed in space. This assumption is nonetheless employed widely, for example, whenever the commonly used e-folding or tidal prism flushing times (which are also based on the well-mixed assumption) are estimated (Sanford et al. 1992;Monsen et al. 2002). Despite the frequent realworld violation of the well-mixed assumption (Du and Shen 2016), our aim here was to, first, develop simplified models to reveal the underlying relationships in the constituent mass balance for complex estuaries (which inherently suggests that our models will not be entirely consistent with reality). Second, we aimed to evaluate whether and how well the simple models work despite incongruities between their foundational assumptions and reality. The Chesapeake Bay is the largest estuary in the United States and is not well mixed, thus violating this first key assumption. It experiences persistent stratification and exhibits large vertical (Li et al. 2017) and horizontal (Kemp et al. 2005) variations in nutrient concentration. Nonetheless, the comparisons herein of our model equations with available data suggest that our simple models work quite well as first-order approximations on the annual scale, despite the violation of this key assumption. Readers should note that, the further their real estuary departs from the assumptions employed here, the less well the models should be expected to characterize the real estuary. We would expect our simple models to perform better when applied to small, shallow waterbodies, which may better satisfy the wellmixed assumption. Potential refinement of our approach could employ a two-(or more) box model.
The second fundamental assumption above essentially means that a single value for the specific net removal rate of the substance (k) applies at all locations throughout the estuary. This assumption, too, may depart substantially from reality as physical, chemical, and biological conditions influencing internal sources and losses of a constituent can vary significantly in space. This assumption has nonetheless been employed by other authors (Vollenweider 1975;Dettmann 2001). As with the first assumption, we know it is violated in the Chesapeake; processes influencing net internal removal (e.g., bottom nutrient fluxes, primary production, denitrification, burial) vary spatially (Boynton et al. 1995;Kemp et al. 2005;Son et al. 2014), likely resulting in a spatially variable net removal rate. Regardless, our simple approach employing this assumption appears to work reasonably well for describing the nutrient balance in the Chesapeake Bay on an annual timescale, justifying the usefulness of the simple models as a first-order approximation.
The third primary assumption is that of a dynamic steady state, that is, that the variable of interest oscillates in the short term around a mean value, which remains constant over longer timescales (Bierman and Montgomery 2014). In the present case, we assumed that estuary volume V and estuaryintegrated substance mass VC-both means over the period of 1 yr-do not vary between periods (e.g., interannually), despite likely oscillations about the mean value at shorter timescales. Because we assume volume is a constant, an assumed dynamic steady state for integrated substance mass VC reduces to the same assumption for concentration C. As the change in mean water volume for the Chesapeake Bay is relatively minor between years (Boon et al. 2010), the dynamic steady state for water volume may be considered to be satisfied on annual scales in that system. However, because of the large interannual variations in hydrodynamic conditions, the individual transport processes contributing to the water balance, as well as residence time, vary from year to year, thus resulting in estimated inflows and outflows with large interannual variations as well. The applicability of the dynamic steady state for material depends on the substance and may be appropriate for certain time periods. For example, Boynton et al. (1995) and Dettmann (2001) assumed that nutrient mass fluxes and net removal together satisfy the steady-state condition for timescales on the order of one or more years, and they successfully employed this assumption for investigating nutrient loss and transport over those timescales. Their approaches imply that nutrient mass integrated over the estuary and averaged over a period T of one or more years may not vary significantly between periods. For the Chesapeake Bay, interannual variability in mean TN concentration is much less than monthly variations (Du and Shen 2017). We calculated the time derivative for integrated substance mass, d(VC)/dt, each year for all simulated years and compared its magnitude to the substance mass export (Q e C) term (data not shown). For 70% of those years, the interannual variation term was less than 5% of the magnitude of the export term, suggesting that the dynamic steady-state assumption was largely satisfied for most years. Therefore, applying our simple models at the annual scale to TN in a system like the Chesapeake is more appropriate than for shorter timescales. When interannual variation of mean concentration is large, the unsteady constituent mass balance equation (Eq. 2) can be used and interannual differences in mean concentration can be incorporated. Readers should therefore be advised that, although seasonal dynamics, for example, of quantities such as TN may be of particular interest, if a dynamic steady state (i.e., relatively constant volume and mean concentration) cannot be reasonably expected to apply between seasons, then the models presented herein may not be appropriate for such an application.
Given the Chesapeake's violation of most of these three fundamental assumptions, their application to that particular system provides a challenging test for our models. Application to other systems-and evaluation of model performancewould be a worthwhile endeavor to better understand the limits of the models' applicability as well as the cases for which they are best suited.
As explained in Supplementary Information S1, an additional requirement underlying our derivations is that the model averaging period T is, at minimum, greater than the residence time. Because the mean residence time for the Chesapeake is 179 d (Table 1; Du and Shen 2016), a seasonal averaging period of a few months would be significantly shorter than the residence time, clearly violating this requirement and motivating our model applications herein to an averaging period of 1 yr. On the other hand, the model could be applied seasonally for a small waterbody that has a short residence time and that reasonably meets the dynamic steady-state assumption at the seasonal scale.
It should be noted that our simple equations also implicitly assume that tidal dispersion is unimportant, and that advection dominates scalar transport. For systems where concentration gradients are large and dispersion-induced transport is significant (Lucas et al. 2006;Martin et al. 2007), the error associated with this assumption could be substantial.
Although the models are derived for an estuary, they are general and can be applied for other waterbodies, such as lakes and river reaches. It is likely in such cases that mean inflow from the adjacent downstream water body (Q in ) tends toward 0, outflow to the adjacent downstream water body Q e tends toward freshwater inflow Q f , and the ocean exchange factor β tends toward unity, thus reducing to the form proposed by Dettmann (2001) for net export:loading and by Vollenweider (1975) for loading (though differing with respect to the specific definitions for residence time employed). Our models can be used for different substances if the assumptions are applicable and parameters can be appropriately determined.
It is further emphasized that the transport time employed herein is the spatial mean residence time, which may differ from freshwater replacement time, freshwater residence time, and flushing time, depending on how these timescales are computed (Dyer 1973;Officer 1976;Monsen et al. 2002). Although a rigorous computation of residence time according to the definitions employed here may not always be feasible, readers are advised that implementation of other transport timescales could lead to significant errors.
Any error in the calculation of estuary mean concentrations (such as that which may have been introduced herein by incorporation of Chesapeake mainstem data only) could affect values for estimated net removal rate k and ocean exchange factor β and thus estimates of net export:loading or loading (e.g., see Fig. 6 for the sensitivity of TN loading to removal rate). The main idea of this paper was not to provide definitive values for parameters and rates, however. Instead, our goals were to develop simple, accessible models that explain underlying processes and to provide a proof of concept with the examples shown. Given the sensitivity of results to net removal rate k, future efforts toward improved estimation of that parameter are warranted.

Implementation of models
The simplest of the models presented herein is the expression for export:import (Eq. 6), where "import" includes input from the ocean in addition to loading from the watershed, atmosphere and other sources, and "export" is gross transport to the ocean. Readers who wish to estimate export:import need only two timescales: mean residence time and the net removal timescale. In this paper, we used a detailed 3D numerical model implementing an adjoint method (Delhez et al. 2004;Du and Shen 2016) to obtain mean residence time, but there are other approaches for using numerical models to compute residence time in a way that is consistent with the definition employed in the present paper (Takeoka 1984;Monsen et al. 2002;Lucas and Deleersnijder 2020). The net removal rate k (from which the net removal timescale is estimated) may or may not be difficult to quantify for a particular estuary and substance. For example, the net removal rate has been estimated for lakes with the joint application of complex 3D and simple mass balance models (Bocaniov and Scavia 2018). Alternatively, measured or estimated rates of nutrient loss processes such as denitrification and burial may be combined and scaled appropriately by (depending on their units) nutrient concentration, volume, and/or other quantities to obtain a specific net removal rate in the units 1/time. We applied that method to the sum of denitrification, burial, and fisheries harvest data (in 10 6 kg N yr À1 ) presented in the 1985-1986 Chesapeake Bay TN mass balance of Boynton et al. (1995;their fig. 11), dividing by estuary volume (Table 1 herein) and by the mean CBP 1985-1986 TN concentration (Supporting Information Table S1). The resulting Chesapeake-specific TN net removal rate, k Ches TN , is 0.0063 d À1 -only 6% different from the estimate developed by the method described in Steps 3a-d in Supporting Information Section S2 (k Ches TN = 0.0067 d À1 ; Table 1). The two calculations are not entirely independent, as they both employ data from Boynton et al. (1995); however, the first method used herein relied on several additional observed and modeled variables.
Estimation of the ratio of net export to loading (Eq. 7) likewise involves the mean residence time and the net removal rate but also requires the ocean exchange factor . To estimate the ocean exchange factor, mean outflow and inflow at the ocean boundary, Q e and Q in , respectively, may be obtained from mean residence time and Eq. 5, and freshwater discharge Q f and Eq. 3 (alternatively, mean inflow and outflow at the ocean boundary could also be extracted from hydrodynamic measurements or model output); estuary-averaged and downstream inflowing concentration can be derived from measurements. As mentioned above, the net removal timescale τ k (reciprocal of the net removal rate) may be difficult to quantify, for example, due to inadequate measurements. If concentration measurements and estimates of net export:loading, mean residence time, and inflow exist, but net removal rate k does not, then the net removal rate may be backed out from adjusted net removal rate K and the ocean exchange factor β using Eqs. 3, 5, and 7, as was done here for Chesapeake TN. Alternatively, Eq. 10 may be used to estimate net removal rate if loading and other parameters are relatively easily obtained by other methods. The relationship we developed to estimate loading (Eq. 8) requires volume, mean residence time, the net removal timescale, and the freshwater renewal timescale (τ f ), as well as the average and inflowing downstream concentrations. Freshwater inflow Q f may be based on measurements or hydrologic models, and all other parameters may be obtained as described above.

Conclusions
We derived simple mathematical relationships providing linkages between import, export, transport, and net removal for a non-conservative substance within an estuary. The relationships presented expand on previous work that explored simple empirical or substance mass-balance based relationships between those terms. For example, similar expressions have been derived previously for lakes and estuaries. However, the mass contribution from the adjacent downstream water body (e.g., the ocean) was either neglected (Dettmann 2001) or not explicitly accounted for (Vollenweider 1975) in those previous studies. The relationships in the present paper were instead derived from two steady-state mass balances-one for substance mass and the other for water. Incorporation of the water mass balance allowed for the explicit inclusion and estimation of water and substance inputs from the ocean to the estuary, which further explains the underlying mechanisms captured by Nixon et al.'s (1996) well-known empirical relationship for net export:loading as a function of "residence time." We expressed export:import, net export:loading, retention:import, loading, and estuary averaged concentration in terms of timescales.
The key timescales used herein are mean residence time and the timescale for net removal (i.e., the reciprocal of the net removal rate for a substance within the estuary). Mean residence time is the parameter that reflects the underlying transport processes, while the net removal timescale represents net biochemical removal and settling processes. With the use of timescales, the contributions of physical transport and biochemical removal and settling processes to material retention and export in an estuary can be evaluated, as both processes are expressed in terms of a single comparable currency (Lucas et al. 2009). These relationships demonstrate that material retention, export, loading, and mean concentration can be estimated using the mean residence time and the net removal timescale (in some cases also requiring ancillary parameters that are commonly available for estuarine systems).
There are limitations associated with the simple models, as some assumptions employed in the model derivations (e.g., steady state, well mixed) may not always be satisfied in reality. However, we showed that the derived expressions can still be useful in providing first-order approximations for a real estuary. We applied the simple net export:loading model to Nixon et al.'s (1996) previously published lake and estuary data set, achieving performance comparable to those authors' empirical relationships with mean residence time. This indicates that it is feasible to estimate the estuarine net export:loading ratio (which may be otherwise difficult to quantify) based on two timescales and a couple additional commonly collected parameters needed to estimate the ocean exchange factor (i.e., concentration, and freshwater flow). For estimating export:import, only the residence time and net removal timescale are needed. We also applied our simple loading model to the Chesapeake Bay and compared our estimates of TN and TP loading over more than two decades to independent estimates produced by the Chesapeake Bay Program. Our basic TN loading model explained 85% of the variability in the Chesapeake Bay Program's estimated loadings, while an improved TP loading model (casting net removal rate as a function of TSS) explained 62% of the variability in Chesapeake Bay Program's estimated loadings over the period 1990-2012. An accurate estimation of loading may be challenging for large estuaries or other aquatic systems because of complex land use practices, multiple point and diffuse sources (Falconer et al. 2018), atmospheric and groundwater contributions, combined natural and anthropogenic processes (Bricker et al. 2008), and strong temporal variability over multiple timescales (Boynton et al. 1995). Monitoring programs are not always designed to quantify all of these, or even the most important, components (Falconer et al. 2018). The results herein demonstrate that our simplified loading model can be useful for estimating interannual variations in nutrient loadings based on changes in water quality conditions and a few timescales representing flows, transport and removal inside an estuary. Because the derivations of these simple models are based on mass-balances for water and transported material, the results are general and can be applied to different waterbodies and substances if underlying assumptions are applicable.