Factors Controlling Seasonal Phytoplankton Dynamics in the Delaware River Estuary: an Idealized Model Study

Phytoplankton biomass in estuaries is controlled by complex biological and chemical processes that control growth and mortality, and physical processes that control transport and dilution. The effects of these processes on phytoplankton blooms were systematically analyzed, focusing on identifying the dominant controlling factors out of river-induced flushing, tidal dispersion, nutrient limitation, and light limitation. To capture the physical processes related to flow and sediment dynamics, we used the idealized width-averaged iFlow model. The model was extended with a nutrient-phytoplankton module to capture the essential biological-chemical processes. The model was applied to the Delaware River Estuary for the productive months of March to November. Model results were compared with field observations. It was found that phytoplankton blooms cannot form in the lower bay due to tidal dispersion, as water from the estuary and coastal ocean mix in early spring, and due to local effects of nitrogen limitation in summer. In the middle to upper bay, sediment-induced deterioration of the light climate limits the growth but allows for blooms in the mid bay, while no blooms can form in the turbidity maximum zone in the upper estuary. Further upstream in the tidal river, the effects of river-induced flushing dominate in early spring and prevent bloom formation. In the summer and fall, lower river discharges and higher growth rates at higher temperatures allow blooms to form and persist. Analysis of the connectivity between mid bay and tidal river blooms showed that coastal ocean phytoplankton may contribute to mid bay blooms, but do not penetrate beyond the turbidity maximum zone.


Introduction
Phytoplankton biomass is considered one of the main indicators of the health of an estuarine ecosystem, because of its key role in the food web and the oxygen cycle. However, understanding the dynamics of phytoplankton in Communicated by Mark J. Brush the light penetration depth (or euphotic depth). In estuaries, the euphotic depth is often controlled by the amount of suspended sediment (Wofsy 1983;Peterson and Festa 1984). Using local biological-chemical models together with observations of sediment concentrations, critical depth theory has been successfully used to qualitatively describe phytoplankton dynamics in some estuaries (e.g., Colijn 1982;Cloern 1987, and references therein). However, bloom dynamics may not obey the critical depth theory. In oceanic and coastal ecosystems, this is for example observed when greater growth rates of phytoplankton than those of zooplankton initiate phytoplankton spring blooms in the absence of thermal stratification, such that the spring bloom occurs very early in the season (February in the North Atlantic) (Behrenfeld 2010;George et al. 2015). In estuaries, critical depth theory may be violated due to salinity stratification and turbulence, which control the exchange between the euphotic and aphotic layers due to mixing and sinking (Lucas et al. 1998) or due to non-local effects (Lucas et al. 1999b), most notably river flushing (e.g., Filardo and Dunstan 1985;Zakardjian et al. 2000;Liu and De Swart 2015). When conversely assuming that the phytoplankton biomass is fully controlled by flow and dilution, residence time is a useful predictor of estuarine primary production (e.g., Howarth et al. 2000). High residence times are typically associated with higher biomass, as the phytoplankton has time to grow. However, Lucas et al. (2009) explained that the converse may be true if local losses exceed growth or that there may be no apparent relation if local losses and growth are balanced. Clearly, a full understanding of phytoplankton dynamics in estuaries requires a combined insight into local and non-local processes. While both classes of processes are generally built into complex numerical simulation models, the high complexity of such models and the variability on many timescales makes it difficult to distinguish between the effects of local and non-local processes and evaluate their relative importance.
In this study, we developed and analyzed a method to assess and compare the relative importance of local and nonlocal processes to the control of phytoplankton dynamics in an estuary in terms of equivalent growth rates. The method is formulated generally and can be used to compare local and non-local processes in various modelling frameworks. The goal of this modelling study is to provide insight into the main processes that govern phytoplankton blooms on the scale of the entire estuary and on a seasonal timescale. As the focus is on understanding of large-scale dynamics, we used an idealized width-averaged model that extends the iFlow model (Dijkstra et al. 2017a), resolving the dynamics of the transport of water, sediment, phytoplankton, and nutrients. The strengths of the model are the ability to make a further decomposition of the local and non-local processes into individual biological and physical processes and the models' computational speed, taking only seconds to compute long-term equilibrium solutions.
The model is applied to the Delaware River Estuary, and results are calibrated against an extensive set of observations. The current understanding of phytoplankton dynamics in the Delaware River Estuary is based on observations and conceptual local models (Pennock 1985;Pennock and Sharp 1994), where the effect of non-local flow-related processes has not been considered. Therefore, in this study, we focus on the relative importance of local and non-local processes during various seasons.
A brief introduction to the study area and the model and analysis methods developed in this study are described in "Model and Site Description". Month-to-month model results covering the entire year are presented in comparison to observations in "Year-Round Results for Default Settings" and analyzed in the context of local and non-local processes in "Balances and Limiting Factors" to "Synthesis". The model results and assumptions are discussed in context of other literature on the Delaware Estuary and in context of its general implications in "Discussion". Finally, the conclusions are summarized in "Summary".

Study Area: Delaware River Estuary
The Delaware River Estuary is located on the east coast of the USA (Fig. 1). Tides in the Delaware propagate from the mouth at Cape May and Cape Henlopen to Trenton, 215 km upstream, beyond which the tidal influence disappears. Several tributaries flow into the Delaware Estuary, of which the Schuylkill River (km 149) is the most significant. The Delaware River is well monitored with long-term data on tidal elevation available from eight NOAA tide gauge stations, information on the river discharge available from USGS, and data on suspended sediment and biochemical quantities available from several sources. Biochemical data at the surface has been gathered in several cruises by researchers from the university of Delaware in the 1980s and 1990s (see Sharp et al. 2009, for an overview). McSweeney et al. (2016b) present data from several cruises in 2010-2011, where the distribution of sediment, oxygen, chlorophyll-a, and nitrate over the depth has been measured along the length of the channel. The most extensive long-term data set is collected by the Delaware River Basin Commission (DRBC). The data set consists of measurements taken approximately every month since 1967 at the surface at 22 stations along the estuary and tidal river and includes observations of suspended sediment concentration, temperature, salinity, chlorophyll-a, The zone indication is adapted from Sharp et al. (2009) and several nutrients. In our study, we use the DRBC observations from 2000 to 2016, which are well documented and available online. 1

iFlow Model and Application to the Delaware
The nutrient-phytoplankton model used in this study extends the iFlow model (Dijkstra et al. 2017a). This is a process-based width-averaged idealized model for tidal hydrodynamics and sediment dynamics in estuaries and tidal rivers. The geometry of the estuary is represented in the model by a smooth width and depth profile capturing the estuary-scale features. A perturbation method is used to obtain an approximate solution to the nonlinear equations for water motion and sediment dynamics. The approximations resulting from the perturbation method lead to short computation times, allowing for sensitivity studies investigating the sensitivity of model results to uncertain or variable model parameters. Moreover, the model allows for an explicit decomposition of the water motion, sediment concentration, and sediment transport into different components resulting from various physical processes, thus allowing for a mechanistic interpretation of the results. Finally, the model immediately solves for dynamic equilibrium conditions, i.e., the condition that develops after a long time of constant forcing, and therefore quickly gives insight into long-term trends. Hence, numerical timestepping routines are not needed. The model uses entirely analytical solutions in the vertical direction. Solutions in the horizontal direction are numerical on an equidistant grid 1 http://www.state.nj.us/drbc/quality/datum/boat-run.html with 250 computational cells. In this section, we give a short overview of the physics included in the model and the application to the Delaware Estuary. For a detailed description of iFlow, we refer to Dijkstra et al. (2017a).

Geometry and Water Motion
The geometry and water motion is adopted from the work by Wei et al. (2016). Consistent with their approach, we approximate the estuary-scale geometry by a constant width-averaged depth of 8 m and an exponentially converging width B, according to B = B 0 e −x/L b , where x is the along-channel axis, B 0 the width at the seaward boundary of 39 km, and L b a convergence length of 42 km.
The water motion is described by approximations of the width-averaged momentum and continuity equations. The model resolves the M 2 , M 4 and tide-averaged water level and width-averaged flow velocity in dynamic equilibrium. The model is forced by a prescribed M 2 and M 4 water level amplitude at the estuary mouth with an amplitude of 0.72 m and 0.14 m, respectively, and a phase difference between the M 2 and M 4 water level of −152°. This forcing equals the year-averaged conditions measured at the NOAA tide gauge at Cape May in 2016. The water level inside the estuary is calibrated against seven other NOAA tide gauges by adjusting the bed roughness coefficient and eddy viscosity in the model. Fresh water enters the estuary at the landward boundary at Trenton and at the confluence with the Schuylkill River. The discharge is represented by the monthly average of the discharge measured by USGS at Trenton between 2000 and 2016. Observations by USGS of the Schuylkill River discharge show that this is 25% of the discharge at Trenton on average and we therefore represent the discharge of the Schuylkill River as 25% of the discharge at Trenton. The two discharges combined vary between 300 m 3 /s on average in July and 727 m 3 /s in March (see Table 1 for an overview per month).

Sediment Dynamics
The sediment model includes settling, resuspension, and along-channel advection and diffusion of a single sediment fraction. Sediment is represented as fine non-cohesive sediment with a prescribed settling velocity of 0.5 mm/s, based on model calibration and fitting in the range of observed settling velocities (Cook et al. 2007). The resuspension rate for sediment was chosen sufficiently high so that the muddy bottom pool is formed and depleted over the course of a tide and not growing on a subtidal timescale (i.e., availability limited conditions, for such conditions, the exact value of the resuspension parameter is irrelevant to the model (Brouwer et al. 2018)). The models return approximations to the M 2 , M 4 and tide-averaged signals No such rating curve is available for the Schuylkill River, but Delaware Estuary Regional Sediment Management Plan Workgroup (2013) shows that the long-term sediment discharge of the Schuylkill River between 1950 and 2009 is 70% of that of the Delaware River. Hence, we use the same rating curve as for Trenton, but with a coefficient of 1.8 · 10 −5 instead of 2.5 · 10 −5 .

Salinity
Salinity in the model is assumed to be well mixed in the vertical and constant over the tidal cycle. This simplified representation of the salinity profile is sufficient to capture the subtidal density gradient induced by the salinity gradient and hence the flow by gravitational circulation. Effects related to periodic stratification or strong stratification are not captured by the model. Dijkstra et al. (2017b) and Burchard and Hetland (2010) estimated that periodic stratification leads to an amplification of the gravitational circulation of approximately a factor of 2 for estuaries like the Delaware. Hence, to parameterize these effects, the gravitational circulation was increased by a factor of 2. The expression for the salinity distribution, s, along the distance of the estuary, x, reads: where model parameters s sea = 31 psu, x L = 42 km, and x c = 100 km follow from calibration to surface salinity data of DRBC collected between 2000 and 2006. This salinity distribution is related to the river discharge at Trenton using a power law ∼ Q −1/7 Trenton (Monismith et al. 2002;Aristizábal and Chant 2013).

Nutrient-Phytoplankton Model
In order to describe the phytoplankton dynamics, a width-averaged nutrient-phytoplankton module was added to iFlow. This nutrient-phytoplankton model includes a biological-chemical component and an advection-diffusion component, which are briefly introduced in this section. For an extensive and mathematical description, we refer to the Supplemental Material to this article. The model uses two classes of nutrients: dissolved inorganic nitrogen and phosphorus (DIN and DIP). The DIN fraction represents nitrate, nitrite, and ammonium, while the DIP fraction represents phosphate. The model uses one class of phytoplankton that represents the entire phytoplankton community found in the estuary. The response of this phytoplankton class to environmental conditions is controlled by representative aggregate parameters. Like the model for water motion and sediment concentration, the nutrient-phytoplankton model describes an equilibrium state, i.e., the state that is attained after a sufficiently long time of constant forcing, thus representing long-term trends instead of transient behavior.
Throughout this study, we express phytoplankton in terms of its Chl.-a content (in μg Chl.-a), DIN in mol N, and DIP in mol P. To convert between these units, we assumed constant conversion rates based on the Redfield ratio, i.e., C:N:P equals 106:16:1 (in mol) and a constant C:Chl.-a ratio of 50 g C/g Chl.-a. This ratio is a reasonable estimate for average conditions given the range of values reported for estuaries (e.g., Cloern et al. 1995). As a result, the Chl.-a:N ratio equals 1.6 μg Chl.-a/μmol N, which is within the range of 0.9-1.8 reported by Sharp et al. (2009) for the Delaware.

Biological-Chemical Component
The biological-chemical component of the model is sketched in Fig. 2 Fig. 2 Schematic representation of the biological-chemical model component with phytoplankton, nutrients and pools of organic nutrients. All sinks to the phytoplankton biomass are parametrised by a mortality rate. As the model computes equilibrium conditions, the uptake flux (1.) equals the phytoplankton sinks (2. and 3.) and nutrient remineralization (5. and 6.), so that the pools of organic nutrients do not have to be resolved using a growth rate μ that depends on temperature, light, and nutrient availability. All sinks in the phytoplankton biomass, including mortality, grazing, respiration, and sinking are parameterized by a loss rate m. In the remainder of this study, we refer to this loss rate as the mortality rate, as is conventional, even though m includes more than only mortality. This mortality parameter is treated as a calibration parameter. The organic nutrients originating from these phytoplankton sinks are represented by suspended and bottom pools, where they are remineralized to inorganic forms that can be taken up. Since the model describes an equilibrium state, all fluxes in the model are in balance. Hence, uptake (1. in Fig. 2) equals the phytoplankton sinks (2. and 3.) as well as the remineralization (5. and 6.). It is additionally assumed that transport of the organic bottom pool may be ignored. An important consequence of these assumptions is that the amount of suspended and bottom organic nutrients, the sediment nutrient flux (i.e., 5. in the figure) and the time required for remineralization are irrelevant and do not have to be resolved by the model. Hence, we only explicitly resolve DIN, DIP, and phytoplankton. It is known that time lags related to remineralization of the bottom nutrient pool may be important on the monthly timescale, so that the equilibrium assumption is only an approximation. This is discussed further in "Discussion".
The growth rate μ depends on the temperature, light availability, and nutrient availability according to: In this equation, μ max (T ) is the temperature-dependent maximum growth rate, N (1) is the DIN concentration, N (2) is the DIP concentration, E is the photosynthetically active radiation (PAR), and H (1) N , and H E are saturation parameters. The brackets · denote averaging over the tide and day. Hence, the minimum function in the equation is evaluated at each instance of time, taking the most limiting out of the instantaneous DIN, DIP, and PAR, and is then averaged over time.
The maximum growth rate μ max is described following Eppley (1972) and using average monthly water temperatures T (Table 1) The nitrogen and phosphorus limitations are described by Michaelis-Menten formulations with saturation coefficient H (1) N = 3 μmol N/l (Banas et al. 2009;Eppley et al. 1969;MacIsaac and Dugdale 1969) and H (2) N = 0.2 μmol P /l (using Redfield ratio). The light limitation is described by a different saturating function (e.g. Smith 1936). This function has a saturation parameter H E = 110 μmol photons/(m 2 s), based on community averaged insitu incubations in the Delaware by Harding et al. (1986) (see also Supplemental Material). The PAR depends on the vertical coordinate z, time t and the phytoplankton and sediment concentrations P and c according to Here, E 00 (t) is the seasonal variation of the maximum daily light availability and is determined using PAR measurements from 2016-2017 by the National Ecological Observatory Network (NEON 3 ) at the Smithsonian Environmental Research Centre (SERC), MD approximately 100 km from the Delaware Estuary (Table 1). The function d E (t) accounts for daily variations in PAR and seasonal variations in day length, also determined using NEON data. The function α describes light attenuation according to the Lambert-Beer law This includes a background attenuation coefficient k bg = 0.095 1/m (Pennock 1985), sediment-induced light attenuation with coefficient k c = 50 m 2 /kg and shading by phytoplankton (i.e. self-shading) with coefficient k p = 18 m 2 /mol N (Pennock 1985;Banas et al. 2009). Literature values for k c for the Delaware River were originally derived for light attenuation models driven by measured surface sediment concentrations. We corrected the value of k c using modelled vertical sediment profiles, such that the depthaveraged light attenuation is the same as when using surface concentration with a surface attenuation coefficient of 60 m 2 /kg (Cloern 1987;Sharp et al. 2009).

Advection-Diffusion Component
The model includes advection of nutrients and phytoplankton with the tidal and subtidal flow as well as diffusion with a prescribed diffusion coefficient K h of 100 m 2 /s (Wei et al. 2016). In the vertical direction, an eddy diffusivity was used to describe vertical mixing of nutrients and phytoplankton. Phytoplankton was additionally assigned a prescribed settling velocity w p equal to 1 m/day (Sarthou et al. 2005). At the bed, it was assumed that any live phytoplankton that settles is immediately resuspended. Phytoplankton losses by settling and burial or benthic grazing were parameterized by the mortality rate.
In order to force the model, depth-averaged timeaveraged nutrient and phytoplankton concentrations were prescribed at the seaward boundary. The nitrogen concentration at this boundary was set to 0, as nitrogen concentrations are negligible at the ocean compared with the concentrations inside the estuary. The phosphorus concentration at the boundary is set to 1 μmol P/l, informed by measurements. For phytoplankton, we imposed a small phytoplankton concentration, P sea = 0.1 μg/l Chl.-a. This is so small compared with typical bloom concentrations of 30-40 μg/l that it can be interpreted as a minimal background condition required for the model to develop phytoplankton growth. At the upstream boundary, fluxes of DIN Q (1) N (in mol N/s) and DIP Q (2) N (in mol P/s) into the estuary were based on the measured influx at Trenton, according to (Buxton et al. 1999): Based on annual average measured Chl.-a concentrations and the average river discharge at Trenton, the upstream input of phytoplankton, Q P , was fixed at a value of 1300 μg Chl.-a/s. Additional sources of nutrients were added to the model to obtain a reasonable representation of the measurements as input to the phytoplankton model. We chose to impose sources of DIN and DIP at the confluence of the Delaware and Schuylkill Rivers in Philadelphia, representing nutrients flowing into the estuary from the Schuylkill River and effluents from the city of Philadelphia (Lebo and Sharp 1993). The DIN source is represented as a constant point source of 30 mol N/s, irrespective of the discharge. The DIP source is represented as a constant point source of 3 mol P/s. Lebo and Sharp (1992) describe that a part of the phosphorus source is in the form of particulate material that is remineralized into DIP at a different location in the system. For simplicity and since this source is not quantified, we did not take such a particulate phosphorus source into account.

Method for Analyzing Local and Non-local Processes
In order to distinguish between local and non-local mechanisms, we consider the cross-sectionally integrated tide-averaged phytoplankton dynamics equation, which reads (Qin and Shen 2017): where A is the cross-sectional area, u is the alongchannel flow velocity, K h the horizontal eddy diffusivity, x and z denote the along-channel and vertical coordinates, and t denotes time. Angular brackets · are used for time averaging. The first term in the equation denotes the variation of P on a long timescale and equals 0 as we considered equilibrium conditions. The remainder of the left-hand side denotes the non-local terms related to advection and diffusion. The terms on the right-hand side denote the local terms describing the biologicalchemical component of the model. The equation thus describes a balance between the local and non-local terms when considered in equilibrium. The local and non-local processes scale with the phytoplankton concentration, so that these terms are large in phytoplankton blooms and much smaller outside of these blooms. Moreover, they scale with the cross-sectional area, which can vary strongly along the estuary. For interpretation purposes, it is therefore more practical to convert the local and non-local processes into an equivalent growth rate. This is done by dividing the equation by the depth-averaged tidally averaged phytoplankton concentration P and cross-sectional area. This yields the following equivalent growth rates (with units 1/day): This decomposition is the same as that used by Qin and Shen (2017). For further analysis, the local growth rate μ is separated into several effects. First, it is noted that the growth is equal to 0 at night due to light limitation. As this is a rather trivial limitation, light limitation at night is not accounted for in the decomposition. This is achieved by averaging the equation for the growth rate (Eq. 3) over daytime and nighttime, rewriting it as: The light limitation at night is now accounted for in the factor τ day τ (i.e., time between sunrise and sunset divided by total day length) and · day indicates averaging over daytime conditions. The limitation to the growth rate is decomposed by evaluating each of the growth-limiting terms M N (1) , M N (2) , M E following the example by e.g. Cerco and Cole (1994). These terms were defined as (cf. Eq. 3): Up to leading order, the functions for N (1) and N (2) are constant over a tidal or daily cycle (see Supplemental Material). The function for light limitation varies over the tidal and daily cycles due to the day-night cycle and tidally varying sediment concentration.
The term M E is furthermore separated into contributions by the daily variation of light and different light attenuation factors: background attenuation, sediment shading, and self-shading. Due to the strong non-linearity of the light limitation function, there is no unique way of making such a separation. The method used here is presented in the Supplemental Material. As a result of the perturbation method used in iFlow, the non-local terms in Eq. 10 can also be separated further into contributions by the tide, river discharge, diffusion, and several non-linear effects, as explained by Dijkstra et al. (2017a).
Below, results for the sediment, DIN, DIP, and phytoplankton concentrations in March to November are compared with DRBC observations from 2000 to 2016. It has to be noted that the calibration data set is the same as the data used for comparison of the results. This is acceptable for our purposes as we focus on a qualitative comparison of patterns and on the underlying balance of local and non-local processes. For the qualitative comparison of patterns, we focus on the relative importance of different estuary-scale phytoplankton blooms and on month-to-month variations, which do not follow trivially from the calibration procedure. This is discussed in "Year-Round Results for Default Settings". The underlying balance between local and non-local processes is discussed in "Balances and Limiting Factors". The model is additionally used to draw some conclusions about the connectivity of blooms in the Delaware Estuary in "Sources of Phytoplankton". The results are presented in synthesis in "Synthesis".

Sediment
Model results of the surface concentrations per month are shown in Fig. 3. The month-to-month variations result from a difference in river discharge, which in turn affects the salinity field and the sediment discharge. The model results are compared with surface sediment data collected by the DRBC between 2000 and 2016. The model produces a clear ETM between km 90 and 120 with a magnitude of 25 to 40 mg/l. The ETM location is very well captured by the model and the concentrations reflect the overall seasonality of the ETM, although they tend to underestimate the median of the observed concentrations in the ETM. The approximate magnitude of the median concentrations up-and downstream of the ETM is captured by the model as well. Only in March and April are the upstream concentrations high compared to the DRBC measurements. However, these high concentrations do correspond to the values measured at Trenton by USGS. It thus seems there is a discrepancy between the data observed near Trenton by DRBC and USGS, with USGS observing The model results vary month-to-month due to differences in the river discharge higher concentrations than DRBC during high discharge periods.

Nutrients
DIN concentrations are high in the estuary upstream from the lower bay and throughout the year, with median concentrations up to 150 μmol N/l occurring between approximately km 80 and 160 (see Fig. 4, left panel). As the month-to-month variation is relatively small, we present the observations of an entire year in one figure. This is plotted together with a band of model solutions representing the spread in model results from month to month. The overall patterns are captured by the model, showing a fairly constant DIN concentration from the head of the estuary to the confluence with the Schuylkill River, then a rapid increase in concentration and a gradual decrease toward the mouth of the estuary. As the concentrations are much higher than the saturation coefficient H N (1) in most of the estuary, small differences between the model results and observations will have little effect on the phytoplankton concentration.
DIP concentrations are also high throughout the estuary, with median concentrations up to 10 μmol P/l but with a significant spread in the observations (see Fig. 4, right panel). On the whole, the model captures the qualitative trend of relatively higher DIP concentrations in the ETM zone (km 70-115) and relatively lower concentrations upand downstream of this through the entire year. As the observations and model both show DIP concentrations much larger than the saturation coefficient H N (2) , the scatter in measurements and differences with the model results are of little influence on the phytoplankton concentration. The exception to this is found in the lower bay (< km 25). Here, the model forces the DIP concentration to a small value in the model, whereas observations show DIP concentrations anywhere between 0 and 5 μmol P/l. The model value at the mouth typically represents a lower estimate of the DIP concentration, although lower values have been observed.

Phytoplankton
Measurements of phytoplankton concentrations (Fig. 5) show two predominant bloom locations. The first is in the mid bay around km 25-50 and is present during the entire span of the data from March to November. The second bloom is located in the tidal river between km 120 and 180. It appears in April or May and disappears in August or September. The relative importance of the two maxima changes over the year. In March, the mid bay bloom is at its maximum, while the tidal river bloom is absent. Toward the summer, the mid bay bloom becomes less pronounced, while the tidal river bloom develops. After August, the chlorophyll-a concentration in both blooms becomes smaller and has nearly disappeared by October. Few measurements from November are available, so that the magnitude of the blooms in this month is unknown. The modelled chl.-a concentration (red lines in Fig. 5) also shows two bloom locations, which qualitatively capture observed locations and seasonality. The first bloom occurs in the lower-mid bay around km 20-30 from March to November and resembles the observed mid bay bloom. The second occurs in the tidal river around km 160-170 from May to September. This bloom is narrower than the observed bloom. From April to June, the magnitude of the tidal river bloom is underestimated by the model, while the magnitude is approximated well from July to October. The model shows a very strong minimum in the phytoplankton concentration between km 80 and 100 in all months except for March. This minimum is much more pronounced than in the measurements, which show a minimum around the same location but still with chlorophyll-a concentrations typically around 5 μg/l.
Focusing on the relative magnitude of the two blooms in the model, the upstream bloom grows relative to the mid bay bloom between March and July. While the same seasonal behavior appears in the observations, this is exaggerated and delayed in the model result. From August to October, both modeled blooms show similar maximum concentrations decreasing with time, similar to the observations.

Balance of Equivalent Growth Rates
The physical and biological-chemical processes underlying the results are analyzed by expressing them in terms of equivalent growth rates (cf. Eqs. 10-11). As the simulation represents steady-state conditions, the sum of all contributions to the equivalent growth rate equals 0. It is therefore not the absolute magnitude of the contributions that matters, but their relative magnitude. This relative importance of different processes gives information about the main factors that control blooms and the sensitivity of the results to model parameters.
A useful way of viewing the balance between non-local and local processes is in terms of the residence timescale versus the growth-mortality timescale (Lucas et al. 2009). The residence timescale is a complex function of the nonlocal flow-related processes, as the flow may either lead to flushing of phytoplankton or accumulation at different locations and times. The growth-mortality timescale depends on all the factors that affect growth and mortality and is a complex non-linear function of nutrient and light availability. The equivalent growth rates investigated here provide an insight into the relative importance of these timescales. If non-local processes are of a similar magnitude or large compared with local processes, the residence time can be important compared to the growth-mortality timescale. If in this case the local growth rate > loss rate, bloom concentrations can be restricted by flushing or reinforced by an inflow of phytoplankton from elsewhere. If conversely loss rate > growth rate, the phytoplankton may still persist because of the throughflow of phytoplankton from elsewhere. If local processes dominate over non-local processes, an equilibrium state develops primarily because of a balance between local growth and losses. This balance is controlled by the phytoplankton concentration. If growth rate > loss rate, an increase of the phytoplankton concentration leads to a depletion of nutrients and self-shading, so that the growth rate can balance the loss rate. Additionally, if the increase of the phytoplankton concentration is local, a large along-channel gradient of the phytoplankton concentration is created. This leads to an increasing non-local transport of phytoplankton biomass out of the bloom zone. If conversely loss rate > growth rate, a decrease of the phytoplankton concentration leads to a decrease in nutrient consumption and self-shading until a balance is achieved or until all phytoplankton have disappeared.
The balance of equivalent growth rates is illustrated for March and July (Fig. 6). The main balance in March is qualitatively representative for early spring conditions with low temperature and a high discharge, while July Fig. 6 Most important terms in the decomposition of the phytoplankton balance into equivalent growth rates for March and July for the phytoplankton. The growth and mortality terms are local terms; the other terms are non-local terms. The peaks around km 150 are artifacts of the localized input of water from the Schuylkill River and should not be considered is representative for summer and fall, with moderate to high temperatures and a low discharge. The figure shows the local growth and mortality (green and orange) and several non-local processes. The term "Diffusion"' (blue) represents effects of parameterized horizontal diffusivity, the "River" (red) term represents flushing by the river discharge, "Tidal return flow" (purple) is the combined effect of Stokes drift, and the resulting return flow and "Tide" (brown) is the net effect of dispersion by the M 2 tide. Positive values in the figure denote contributions to the growth, while negative contributions reduce the growth rate. The figures show some spikes in the equivalent growth rates near km 150, due to the local point source discharging water of the Schuylkill River. We will not further consider these peaks in our analysis.
The March equilibrium phytoplankton concentration is mainly established by a balance of positive local growth (green line) versus mortality (orange line) and river flushing (red line). The latter two are of the same order of magnitude in much of the upper bay and tidal river (>70 km). This means that the residence and growth-mortality timescales are similar. The result of this balance is an equilibrium concentration that does not allow for bloom formation in the ETM zone and tidal river, even though local growth > mortality. In the mid bay, local processes are dominant as the riverine influence decreases, thus allowing for the formation of a bloom. In the lower bay (< km 25), the nonlocal processes again are more in balance with the local processes. This is mainly expressed in the tidal and diffusive terms which act to reduce the phytoplankton concentration by mixing water from the lower bay with phytoplanktonpoor water from the coastal ocean. The river contribution opposes this by delivering water rich in phytoplankton from the mid bay bloom to the lower bay.
In July, the phytoplankton concentration results from a balance that is dominated by local processes along the entire estuary. The dominance of local processes is caused by a combination of a small river discharge and a high growth rate due to a high temperature. Hence, the residence time is large relative to the growth-mortality time. As a result, phytoplankton blooms are found wherever local growth rate > mortality rate, i.e., in the tidal river and mid bay. Consequently, phytoplankton also almost completely vanishes wherever local growth rate < mortality rate, i.e., in the ETM zone. Nevertheless, it would be a misconception to conclude that non-local terms can be omitted from the model. If non-local processes were switched off, phytoplankton concentrations of 100 to 200 μmol/l would be attained, which are not realistic. The blooms have set up along-channel gradients in phytoplankton concentration, which lead to some non-local transport. As local growth and mortality almost equilibrate, this small non-local transport closes the balance.

Limiting Factors to the Local Growth Rate
We further study the mechanisms underlying the local processes using the decomposition of the depth-averaged net growth rate (see Method for Analyzing Local and Non-local Processes). This is illustrated in Fig. 7 for March (left panels) and July (right panels). The numbers in the figure should be interpreted as the fraction of reduction of the growth rate μ max during daytime (dawn to dusk): 1 means no reduction, 0 means reduced to zero growth. These reduction factors are related to the limitation by DIN, DIP, light, and mortality. The reduction factor caused by the mortality rate computed relative to the product of the maximum growth rate and day length μ max τ day τ . The light limitation varies during the day due to tidal variations of water depth, sediment concentration, and daily variation of solar irradiance. This is visualized in the top panels by the red-shaded area showing the light limitation that occurs between 1 and 99 percentiles of the time. The lower panel shows the same results, but with a further decomposition of the light limitation for time-averaged conditions into effects of background shading, sediment, self-shading, and the daily cycle (i.e., variation of irradiance between sunrise and sunset; light limitation at night is not included).
In March, the growth rate is dominantly limited by light availability in most of the estuary. In the top panel, even the upper edge of the red band is more limiting than any of the other factors for x > km 15, indicating that light is limiting at each instance of time. Hence, light limitation still dominates in the midday around slack tide, when sediment concentrations are relatively low. The bottom panel shows that the light limitation is mostly caused by the daily light availability (i.e., variation of irradiance between sunrise and sunset) and sediment shading. The effect of sediment shading alone already exceeds the effects of mortality and nutrient limitation upstream from the lower bay. Furthermore, self-shading is nearly as important as sediment shading in the mid and lower bay. The phytoplankton growth becomes dominantly nitrogen limited only in the most downstream 10 km in the lower bay.
In July, the effect of light limitation still dominates the growth rate in most of the estuary and during the entire day, even though the light limitation is smaller in absolute sense than in March. The light limitation is smaller than in March mainly because the daily light availability is larger (i.e., larger maximum irradiance). Sediment shading is of similar importance as in March in the ETM zone but less in the tidal river due to lower sediment concentrations (cf. Fig. 3). The mortality rate has become a more important limitation to the net growth rate compared with that in March, indicating that processes including grazing and respiration, not explicitly resolved by the model, have become relatively more important. Nitrogen limitation remains the dominant limitation in the most downstream part of the estuary.
The resulting net growth rates, μ − m, and their structure in the vertical direction at three locations along the estuary are plotted in Fig. 8. The net growth rates are positive near the surface as the light limitation by sediment and self-shading vanish near the surface. As the light limitation increases with depth, the net growth rate decreases to a value below 0. The depth of zero net growth (horizontal dotted lines) is around 2 m in the ETM zone and 3-4 m outside the ETM zone, which implies that net growth can only occur in less than half of the water column in the entire estuary. Nevertheless, the net growth averaged over the water column (closed circles in the figure) is positive in most locations. Even in the ETM zone in March, the net growth rate is just positive. In July, the net depthaveraged growth in the ETM zone is negative, leading to decay of the phytoplankton concentration. Figure 8 also shows the vertical profiles of the Chl.-a concentration. This is well mixed over the water column, with a slightly smaller concentration near the surface. This is a consequence of Fig. 8 Vertical profiles of the net growth rate μ − m (top) and phytoplankton concentration relative to the depth-averaged phytoplankton concentration P (z)/P (bottom) for March and July at three locations along the estuary: in the mid bay bloom (km 30), in the ETM zone (km 100), and in the tidal river bloom (km 170). In the top row, the horizontal dotted lines indicate the zero crossings, i.e., the transition from net growth to net decay. The dots indicate the value of the depth-averaged net growth rate. Note that the plotted vertical profiles of the phytoplankton concentration are identical for March and July the relatively large vertical mixing, combined with a small settling velocity. As the growth rate is much smaller than the vertical mixing rate, vertical variations in the growth rate were estimated to have little effect and were not taken into account in the computation of the vertical Chl.-a profile. Hence, the larger growth rate near the surface than near the bed is not reflected in the vertical phytoplankton distribution.

Sources of Phytoplankton
The phytoplankton model includes two sources: at the seaward and landward boundaries. A better understanding of the pathways of the phytoplankton through the estuary is obtained by further investigating these sources separately. This is done in the model by setting the phytoplankton concentrations at one of the boundaries equal to 0. Figure 9 shows the phytoplankton concentration in March and July when only the source at the seaward boundary is taken into account (i.e., upstream phytoplankton concentration equal to 0). The model now only reproduces the mid bay bloom. Similar results are found for all months and this is insensitive to the magnitude of the seaward phytoplankton concentration (not shown). This result is found because there is either a net local loss in the ETM zone or a small net local growth combined with a small residence time of the marine phytoplankton. This indicates that the marine phytoplankton hardly penetrates into the fresh tidal river, regardless of its species-specific characteristics (salt tolerance, growth rates, nutrient uptake, etc.) and whether it would be outcompeted by more specialized freshwater phytoplankton.
When conversely only using the source at the landward boundary (i.e., downstream phytoplankton concentration equal to 0), the same result is obtained as in the default case (i.e., Fig. 5). The phytoplankton from upstream provides the initial source for growth of the tidal river bloom, then flows downstream and largely dies in the ETM zone. However, a sufficient population makes it through the ETM to also contribute to the mid bay bloom. These results only hold when assuming that the freshwater phytoplankton from upstream is sufficiently salt tolerant and is not outcompeted by more specialized species. These assumptions should at least be called questionable. Therefore, the main conclusion is that it is important to account for speciesspecific characteristics when investigating the spreading of freshwater phytoplankton species in the estuary, while marine phytoplankton are prevented from spreading into the freshwater zone by the ETM and river flow. Insufficient data on species composition is available to the authors to verify these conclusions with data.

Synthesis
Combining the results from Figs. 6 and 7, we propose a two-step analysis for gaining insight into phytoplankton distribution. As a first step, the equivalent growth rates are used to determine whether local or non-local terms are important. If local terms are important, a bloom may develop if light and nutrient conditions are favorable. In that case, it is interesting to investigate the limiting factors to the local growth rate (cf. Fig. 7) as a second step. The two steps for early spring (March) and summer (July) are summarized in Table 2, where the second step is printed in italics if it is not considered to be very important to the end result. Focusing on the second step of the analysis, it is concluded that sediment shading is the most important factor throughout the estuary. Self-shading is also important at the locations of the blooms and, together with non-local dispersion processes, acts to restrict the maximum-occurring phytoplankton concentration. Limitation by nitrogen occurs in the lower bay. However, in spring, this limitation is dominated by the effects of non-local tidal dispersion that mix water from the estuary and coastal ocean. Hence, nitrogen limitation is actually only relevant in summer and fall when the temperature is high and hence local growth is large compared to the effects of tidal dispersion.
Another representation of the bloom dynamics is sketched conceptually in Fig. 10. This figure summarizes the main limitations mentioned in Table 2 and additionally sketches the pathway of phytoplankton through the estuary. Fig. 9 Surface phytoplankton concentration expressed as chlorophyll-a concentration (in μg/l) for March and July as in Fig. 5, but only using a downstream source of phytoplankton

Discussion in Context of Delaware Bay
Nutrient Limitation in Delaware Bay Pennock and Sharp (1994) used observations of nitrogen, phosphorus and carbon ratios, enrichment experiments, and a simple light limitation model to assess the most limiting factor to the local growth rate in Delware Bay. They concluded that nutrient limitation can sometimes dominate over light limitation. In the mid bay, they indicate that phosphorus may be limiting in early-late spring and nitrogen might be limiting in summer. Yoshiyama and Sharp (2006) also analyzed observations and found that light limitation is usually dominant in the mid bay, except for during the spring bloom. In our results, we used lower estimates of the nitrogen and phosphorus levels by assuming little to no nutrients in the ocean, yet our results indicate that sediment-induced light limitation is much more important Fig. 10 Schematic representation of long-term bloom conditions in Delaware River Estuary. Characteristic are two blooms: in the mid bay (40-60 km) and in the tidal river (> 115 km). The ETM is typically located around 100 km. At I, fluvial phytoplankton enters the tidal river. If the temperatures are sufficiently high and river discharge is sufficiently low, this phytoplankton resides in the tidal river long enough to grow and bloom. Much of the phytoplankton that flushes downstream through the ETM (II) dies due to sediment-induced light limitation, unless the river discharge is high and temperature is low.
Whether the fluvial phytoplankton can survive downstream in the mid bay depends on its salt tolerance and other characteristics that are not explicitly modelled in this study. In the mid bay, tolerant phytoplankton from the tidal river may grow again (III). Their growth location coincides with the growth location of the marine phytoplankton that is brought into the estuary with the tides (IV).The dynamics in the lower bay is dominated by tidal dispersion mixing water from the estuary and the coastal ocean in early spring and nitrogen limitation in summer than nutrient limitation in the mid bay. Moreover, it is found that light limitation is still dominant during the middle of the day during times with the lowest sediment concentrations. Hence, the results of Pennock and Sharp (1994) do not agree with those found here, and the nutrient limitation during the spring bloom found by Yoshiyama and Sharp (2006) is not identified. As one possible explanation for this disagreement, Yoshiyama and Sharp (2006) discuss that the enrichment experiments may underestimate light limitation due to sediment settling during the experiment. Furthermore, our study focuses on monthly averages over many years, which may not be representative of the rapidly varying conditions during individual spring blooms.
In the lower bay, Pennock and Sharp (1994) indicate that nitrogen is limiting in early spring, nitrogen and phosphorus are both limiting in late spring, and possibly nitrogen is limiting in summer. Our results support the potential occurrence of nitrogen limitation in the lower bay. Phosphorus is not limiting on average, but observations show significant variation in phosphorus concentrations, which could result in phosphorus limiting during some part of time. A further source of uncertainty is the ratio between nitrogen and phosphorus taken up by phytoplankton (Pennock and Sharp 1994), which we assumed to be constant. In summer, we found nitrogen limitation to be the main limitation. In early spring, we found that, while nitrogen is limiting, tidal dispersion and low temperature are the dominant reasons for not finding a phytoplankton bloom in the lower bay.

Stratification and Three-Dimensionality
Over the seasons and within a spring-neap cycle, the mid bay and ETM zone in the Delaware switch between wellmixed and partially stratified states (Aristizábal and Chant 2014). Stratification in estuaries is generally associated with high growth as it prevents sediment from mixing high up in the water column and prevents phytoplankton from mixing to the bed and restricts sediment to the lower part of the water column (Cloern 1987). For the Delaware, Pennock (1985) and Sharp et al. (1986) hypothesize that the variation in bloom magnitude in the mid bay is due to varying salinity stratification, but the importance of stratification to the phytoplankton dynamics in the Delaware has not been proven.
It is not evident that stratification is important to phytoplankton dynamics in the Delaware. Stratification occurs mainly in the narrow channel, while most of the surface area of the estuary consists of shallow and well-mixed flanks. Hence, the focus on the channel is probably not justified and likely overemphasizes the role of stratification. Our results, not accounting for salinity stratification, show that the average variability of the mid bay phytoplankton bloom from high values in spring to lower values in summer can be explained just by accounting for the differences in the temperature-related local growth and mortality.
Shallow flanks and lateral exchange of water, sediment, and phytoplankton are likely to be important in Delaware Bay. Pennock (1985) shows observations of distinct lateral patterns for chlorophyll-a in summer and McSweeney et al. (2016a) show that lateral processes are important for the sediment transport. Additionally, the lateral dynamics has been shown to be important to understanding the phytoplankton dynamics in other estuaries (e.g., Lucas et al. 1999b). Hence, the three-dimensional dynamics is worth further investigation. Shallow areas on the flanks of the estuary are likely less light limited and can therefore serve as areas of positive growth, whereas the channels may suffer from net loss of phytoplankton if stratification is weak. Hence, a model for this three-dimensional behavior at least needs to capture the large-scale lateral circulation, the lateral distribution of sediment, and formation of stratification in the channel. An idealized three-dimensional model like iFlow has already been developed for sediment and salinity by Kumar et al. (2017) and Wei et al. (2018). Like iFlow, the scaling and solution method in this model is developed for well-mixed estuaries and would require substantial changes to account for stronger stratification.

Model Parameterizations and the Mortality Coefficient
In our model, we have chosen a simple representation of the biological-chemical processes by including only phytoplankton and two nutrients. The phytoplankton mortality is represented by a simple linear formulation mP , which accounts for several processes including respiration, pelagic grazing, and benthic grazing. Studies that explicitly include one or more of the processes parameterized by the mortality rate use formulations with a number of highly uncertain parameters and model formulations (see, e.g., Franks 2002;Gentleman et al. 2003;Brush and Nixon 2017). Moreover, it has been demonstrated that a model may display one or multiple equilibrium solutions or autonomous periodic solutions depending on the choice of model formulation for the grazing formulation (Steele and Henderson 1992;Edwards and Brindley 1996). Therefore, these processes can only be resolved reliably when these parameters and formulations can be constrained by data. This is a problem for grazing, where the only study known to the authors on zooplankton numbers and species along the entire Delaware Estuary is by Cronin et al. (1962). A downside of the linear mortality formulation with spatially constant, m, is that spatial gradients in mortality, e.g., due to spatial gradients in grazers, are not covered. As it is expected that grazers are more abundant where phytoplankton concentrations are highest, we could have overestimated the mortality in the phytoplankton-poor areas, such as the ETM zone. This effect is more quantitative than qualitative and does not affect the qualitative conclusions of this study.
In the current model, the mortality coefficient has been calibrated for each month to give the best fit to the median of the observed phytoplankton biomass. The resulting mortality coefficients are plotted against the water temperature in Fig. 11. The mortality coefficient shows what seems like an exponential function, approximated by: Even though this result simply follows from a calibration, there is a strong dependency on the temperature and no such clear seasonality seems to exist with other input variables, such as the discharge or the light intensity. Also, assuming grazing is important, the trend shown in Eq. 13 seems to be supported by formulations used in the size-structured NPZ model of Taniguchi et al. (2014) and Cloern (2017). They use μ ∼ 1.049 T and m ∼ 1.095 T , resulting in a ratio μ/m ∼ 0.96 T . We find the same ratio when combining Eqs. 4 and 13, i.e., both formulations yield the same ratio of the growth and mortality rate. This observation suggests that a relationship like Eq. 13 between m and T might be more universally applicable and sets a fixed ratio between the maximum growth and mortality. However, this observation needs further investigation. The magnitude of the mortality found in this study is of the same order of magnitude as grazing rates measured by Sun et al. (2007) for Delaware bay (only microzooplankton) and White and Roman (1992) for Chesapeake Bay. Sun et al. (2007) found an average grazing rate of 0.46 day −1 based on measurements in one location in the lower Delaware bay at the end of April. Our mortality numbers for April or May are somewhat lower than this measured grazing rate, especially considering that grazing is only one factor that contributes to mortality. White and Roman (1992) on Fig. 11 Base growth coefficient μ max (red) and calibrated mortality coefficients m (blue) versus the temperature. The blue dotted line shows an indicative trend of the calibrated values of m with the temperature the other hand found depth-averaged grazing rates varying between 0.01 and 0.1 day −1 based on measurements at one location during various seasons, which are lower than our mortality rates. Given the variation in measured grazing rates, our calibrated mortality rates at least seem within a realistic range.

The Equilibrium Assumption
The model assumes dynamic equilibrium conditions. When compared with observations, this assumption implies that transient behavior is negligible. Using the example of the Upper James River, Qin and Shen (2017) show that transient behavior becomes less important compared to local and non-local processes when averaging over longer timescales, starting from a spring-neap timescale. At the monthly scale investigated in this study, Qin and Shen (2017) indicate that transient behavior constitutes less than onethird of the effects of non-local and local processes. This means that transient behavior is small, but not completely negligible at this scale. This means that the equilibrium assumption gives a good first estimate of the processes, but cannot give a full detailed description. Nevertheless, the equilibrium assumption provides a large advantage to the model interpretation as time lags in the remineralization of lost phytoplankton to nutrients become irrelevant (cf. Fig. 2), thus highly simplifying the model.

General Applicability of the Model
The model developed in this study is generally applicable to tide-dominated well-mixed estuaries that largely consist of a single main channel. For example, in studies by Brouwer et al. (2018) and Dijkstra et al. (2019), it has already been shown that observed hydrodynamics and sediment dynamics could be qualitatively reproduced using iFlow in the Scheldt and Ems Rivers and their models could be extended to include phytoplankton. As the model contains few parameters and is computationally inexpensive, it is fast to set up, calibrate, and use as a first assessment tool. As illustrated here using the example of the Delaware River Estuary, the model quickly gives insight into the importance of different aspects including light limitation, nutrient limitation, mortality, and flow throughout the estuary. This basic understanding of phytoplankton dynamics is essential for the development and interpretation of more complex models, as it indicates which parts of the more complex model are most important and should therefore receive more attention. As has been remarked and demonstrated that models with different degrees of complexity can lead to a similar skill in reproducing observations (Franks 2002;Friedrichs et al. 2007), such a systematic approach of using simple models to motivate where and why increasing complexity is necessary can provide a promising strategy to a better choice of model complexity and improved model interpretation.

Summary
We applied a method of explicitly distinguishing between local biological-chemical processes and non-local flowrelated processes governing phytoplankton dynamics in well-mixed estuaries. This method was combined with a newly developed nutrient-phytoplankton module for iFlow, which allows for a further decomposition of the local and non-local processes into specific limiting factors such as nutrients, light, tidal dispersion, and river flushing. The model was used to study phytoplankton dynamics in the Delaware River Estuary as a function of the flow, temperature, nutrient availability, irradiance, and light attenuation due to suspended sediments and self-shading. Average monthly conditions for March through November were simulated and compared with observations collected by the Delaware River Basin Commission (DRBC) between 2000 and 2016.
Model results show that, in early spring, the lower bay (< km 25) is dominated by non-local processes due to tidal dispersion. This leads to mixing of the phytoplankton in the estuary with phytoplankton-poor coastal ocean water, preventing bloom formation. Whereas nitrogen is potentially limiting to the local growth rate, this is of relatively little importance compared to the effect of tidal dispersion. In summer, however, the local growth exceeds the effects of tidal dispersion due to higher temperatures and nitrogen limitation is the main factor limiting phytoplankton growth. In the mid bay (km 25-70), local processes are dominant and allow for the formation of a bloom during the entire period from March to November. The growth rate is limited by sediment shading and self-shading. The local growth rate is also dominant in the ETM zone (km 70-115) during the entire year, but the sediment-induced light limitation is so strong that mortality almost equals or even exceeds growth, hence not allowing the formation of blooms. Finally, in the tidal river (> km 115), phytoplankton dynamics vary with the seasons. In early spring, non-local processes dominate due to a low temperature (i.e., low growth rate) and high river discharge (i.e., short residence time). As a result, blooms are absent. In late spring to early fall, the temperature, and hence the growth rate, is higher and the river discharge is lower, leading to a more dominant role of the local growth rate and bloom formation. The main factors controlling the growth rate are again the sediment concentration and self-shading.
To study the connectivity between the mid bay and tidal river blooms, it was investigated whether sources of phytoplankton from the sea and from the tidal river could contribute to both blooms. The marine source of phytoplankton only appears in the mid bay bloom and cannot penetrate beyond the ETM zone due to the ETM and river flow, regardless of its species-specific characteristics. The tidal river source of phytoplankton is not constrained by the flow and ETM zone and can contribute to both the mid bay and tidal river blooms. However, whether it does contribute to both blooms depends on its species-specific characteristics and could not be verified due to lack of data on species composition.
Based on calibration of monthly values for the mortality rate in the model, an apparent power-relation between the mortality and temperature emerged. Similar relations between mortality and temperature have been obtained using other models and studying other estuaries. Whether this emergent property applies more universally is unclear but is worth further investigation.