Quantifying the Effect of Salinity Stratification on Phytoplankton Density Patterns in Estuaries

To quantify the effect of salinity stratification on phytoplankton density (denoted as P) patterns, experiments were conducted with an idealised model that couples physical and biological processes. Results show that the idealised model is capable of capturing the main features of observed P patterns in the Columbia River estuary during the spring season: during weak stratification, P is almost vertically uniform with values decreasing towards the estuary mouth, whereas during strong stratification, high values of P extend further seawards but are confined to the upper layer. Sensitivity studies reveal that the strong vertical gradients of P can only occur if the intensity of turbulence (measured by depth-averaged values of vertical eddy viscosity and eddy diffusivity) is weak. The advection of P by subtidal currents is important in obtaining a smaller along-estuary gradient of P during weak stratification and in obtaining a smaller horizontal gradient and a larger vertical gradient of P during strong stratification. Accounting for stratification controlled vertical distribution of vertical eddy viscosity and eddy diffusivity is necessary for obtaining realistic P patterns if stratification is strong, but not if stratification is weak. A higher osmotic stress, which leads to faster loss of phytoplankton in salt water, results in a larger along-estuary gradient of P if stratification is weak and in a larger vertical gradient of P if stratification is strong.


Introduction
Phytoplankton density (viz, phytoplankton cell number density, hereafter denoted as P ) in estuaries shows distinct spatial patterns for different strength of salinity stratification. It has been frequently observed that the vertical distribution of P is almost uniform under weakly stratified conditions, whereas if the water column is highly stratified, high densities are only found in the upper layer (Indrebø et al. 1979;Cloern 1996;Watanabe et al. 2014). Also, the horizontal distribution of P is remarkably different under different stratification conditions. A specific example is the Columbia River estuary, where in spring season P during spring tides (weak salinity stratification) significantly decreases along the estuary towards the mouth, whereas during neap tides (strong salinity stratification) high densities extend seawards in the upper layer (Roegner et al. 2011).
Salinity stratification affects vertical turbulent exchange processes, subtidal currents as well as the lysis rate of phytoplankton due to osmotic stress. The latter three aspects are likely to be closely linked to the observed P patterns in the Columbia River estuary. Cloern (1996) pointed out that in estuarine and coastal systems, a strong tidal stirring effectively mixes phytoplankton downward into the aphotic zone where net growth rate is negative, whereas a vertical stratification can isolate phytoplankton in the euphotic zone such that it allows P to increase. Roegner et al. (2011) related the isolation of phytoplankton in the upper low-salinity water during neap tides to the extension of the fluvial conditions to the estuary mouth during strong stratification. Lara-Lara et al. (1990) suggested that estuarine salinity variation also directly influences the P pattern, since it has been observed that a large number of freshwater phytoplankton lyse in low salinity (3-5 psu) regions.
Previous studies did not quantitatively explain the relative importance of the above-mentioned aspects in determining estuarine P patterns under different stratification conditions. In order to study this, account should be taken of the fact that with increasing stratification, the maximum values of eddy viscosity occur closer to the bottom (Geyer et al. 2000) because stratification limits the development of turbulence generated at the estuary bottom (Stacey and Ralston 2005). According to Munk and Anderson (1948), the qualitative behaviour of eddy diffusivity is similar to that of eddy viscosity, i.e. with increasing stratification the values of eddy diffusivity become smaller and their maxima are attained closer to the bed. The impact of the vertical distributions of turbulent exchange coefficients on P patterns, in relation to that of other factors, has not been quantified yet. Moreover, the vertical distribution of eddy viscosity directly impacts the magnitude and vertical structure of the subtidal current (Cheng et al. 2013), which affects the along-estuary distribution of P (Liu and de Swart 2015).
The fact that salinity stratification, vertical turbulent exchange processes and subtidal currents are mutually coupled leads to difficulty in understanding observed P patterns in stratified estuaries and in interpreting results of sophisticated numerical models. To reduce this gap in knowledge, this study is dedicated to quantify the effect of each of the above mentioned aspects on P patterns. Specifically, this is done by fulfilling four aims. The first is to examine the effect of vertical stratification, which largely affects the magnitudes and vertical distributions of eddy viscosity and eddy diffusivity, on the overall characteristics of P patterns, in particular their typical vertical and along-estuary gradients. The second aim is to quantify the importance of advection of phytoplankton and nutrients by subtidal currents on the spatial patterns of P . The third and fourth aims are to assess the sensitivity of the overall characteristics of P patterns to the vertical distributions of vertical eddy viscosity and eddy diffusivity and to the salinity-dependent loss rate, respectively.
To achieve these aims, an idealised model that couples physical and biological processes (Liu and de Swart 2015) was adapted in this study. The motivation to choose an idealised model is that it allows investigation of individual processes in isolation, thereby yielding fundamental insight into the dynamics of the system. Thus, an idealised model is suitable to identify and to analyse mechanisms through which the above-mentioned aspects influence P patterns in stratified estuaries. Moreover, idealised models, also referred to as exploratory models, are fast and flexible, and therefore suitable for sensitivity studies (Murray and Thieler 2004). Furthermore, the Columbia River estuary was chosen as the prototype estuary because sufficient data were available to check the performance of the model.
The hydrodynamic module, turbulence closure and biological module are introduced in "Material and Methods" section. Here, also two parameters are introduced that quantify the vertical and along-estuary gradients of P . In "Results" section, results are first shown for weakly and strongly stratified conditions and they are compared to field data of Roegner et al. (2011) to demonstrate that the model is capable of capturing main characteristics of the observed P patterns. Next, results of sensitivity experiments that relate to the last three research aims are presented. In "Discussion" section, key processes controlling phytoplankton density patterns are identified and results of sensitivity experiments are further analysed. Also, results of additional sensitivity studies (which concern other parameters and boundary conditions, see the Electronic Supplement) are briefly summarised and limitations of the model are discussed. Conclusions are given in "Conclusions" section.

Model
The model used in this study consists of different modules and is an extended version of that of Liu and de Swart (2015). Estuarine subtidal currents are described by a hydrodynamic module. A turbulence closure module (details see "Turbulence Closure" section) that accounts for the dependence of the magnitudes and vertical distributions of eddy viscosity and eddy diffusivity on salinity stratification was adopted, which is one of the new aspects of the present model. Phytoplankton and nutrient dynamics are described in a biological module, which is coupled to the previous two modules.

Domain
The simplified estuarine geometry considered in this study is illustrated in Fig. 1, whose length is denoted by L. The cross section is assumed rectangular with a constant total water depth H and a constant width. A Cartesian coordinate system is used with the x-axis pointing from the estuary mouth (x = 0) to the riverine boundary (x = L) and the z-axis from the water surface (z = 0) downward to the bottom (z = H ). Fig. 1 The configuration of the idealised estuary with constant width in the model. Here, L is the length of the estuary, H is the total water depth and 0 is the origin of the Cartesian coordinate system

Hydrodynamic Module
In this study, the estuarine subtidal current model of Burchard and Hetland (2010) is adapted by extending the formulation of eddy viscosity such that it explicitly depends on the strength of salinity stratification. The width-averaged along-estuary velocity u(x, z) in the present study accounts for the river flow u r induced by freshwater discharge and density-driven flow u d due to along-estuary salinity gradient. The expressions of u r and u d are In these expressions, U r is the depth-averaged velocity of river flow, A v is the tidally-averaged vertical eddy viscosity, ρ * = 10 3 kg m −3 is the density of fresh water and ρ is the density. The latter is related to salinity s by the equation of state: where β is the coefficient of isohaline contraction. Salinity s can be split into depth-averaged s and depthvarying s parts, i.e. s = s + s . Assuming that ∂s /∂x ∂s/∂x (Pritchard 1952), only s is accounted for when calculating the density-driven flow u d . The formulation fors follows that of Warner et al. (2005): Here, s 0 is seawater salinity, x c denotes the position where the salinity is 0.5 s 0 , and x δ is the length scale over whichs varies. At x = x c , the along-estuary salinity gradient ds/dx reaches its maximum. Furthermore, (x c + x δ ) is an estimate of tidally-averaged salt intrusion length. The width-averaged vertical subtidal velocity w(x, z) is obtained by solving the continuity equation and it reads Turbulence Closure The formulation of eddy viscosity A v is that of Chen and de Swart (2016) and it reads in which κ = 0.41 is the von Karman's constant, u * is the tidally-averaged friction velocity, whose value is input into this model, and A z describes the vertical distribution of A v .
Here, A z is modelled as a piecewise function: The vertical shape of A z calculated by the above function is plotted in Fig. 2a. At z = z t , which separates the two sub-domains, continuity of A z requires A I = The value of z t is determined by requiring dA z /dz to be continuous at z = z t . The function of A z contains three independent parameters. Parameter A S is proportional to the value of A v at the water surface, i.e. A v | z=0 = κu * H A S . Furthermore, z 0 is the bottom roughness length and h b is a measure of the height of the turbulent bottom boundary layer. As shown in Fig. 2b, if h b decreases, the maximum of A z decreases and its location shifts closer to the bottom. The turbulent bottom boundary layer is represented by (H − h b ) ≤ z ≤ H (as shown by the grey areas in Fig. 2).
Note that if A S equals the maximum of A z and z t = H , Eq. 9 yields a vertically constant eddy viscosity as was used in e.g. Hansen and Rattray (1965). Furthermore, the choice of A S = 0 and z t = 0 yields a parabolic shape of eddy viscosity used by Burchard and Hetland (2010). The latter two shapes are plotted in Fig. 2a as well.
In this study, A S is tuned such that the amplitude of the modelled subtidal current agrees well with observations, and z 0 is an input parameter whose value is obtained from previous studies. The turbulent bottom boundary layer height h b is calculated from Here, R is the bulk Richardson number (Dyer 1997), g = 9.81 m s −2 is the gravitational acceleration, s is the tidally-averaged bottom-to-surface salinity difference and U T is the amplitude of depth-averaged tidal current. Equation 10 was determined by Chen and de Swart (2016) by fitting Eq. 9 to field data of several estuaries. In this study, U T and s are derived from field data. The formulation of the vertical eddy diffusivity κ v reads The above expression is based on Munk and Anderson (1948), except that here the bulk Richardson number rather than the gradient Richardson number is used. Finally, the along-estuary turbulent diffusivity κ h was assumed spatially constant.

Biological Module
The equations govern phytoplankton and nutrient dynamics are In the above equations, P = P (x, z, t) denotes the widthaveraged phytoplankton density (cells m −3 ) and N = N(x, z, t) is the width-averaged nutrient (e.g., phosphorus, nitrogen or silicon) concentration (mmol m −3 ). Parameter μ is the specific growth rate of phytoplankton, m is the loss rate, v is the sinking velocity, α is the nutrient amount in each cell and is the proportion of dead phytoplankton that is recycled into nutrients. The specific growth rate μ is constrained by local nutrient concentration N and light intensity I . The parameterisation for μ follows Huisman et al. (2006): where μ max is the maximum specific growth rate, H N and H I are the half-saturation constant for nutrient-limited and light-limited growth, respectively, and min denotes the minimum function. Based on the above equation, a limitation index is defined as N H N +N − I H I +I . If this index has negative values, phytoplankton growth is nutrient-limited, whilst if the limitation index has positive values, the growth is light-limited.
The light is supplied at the water surface and its intensity decreases exponentially with depth according to Lambert-Beer's law: Here, I in is the incident light intensity, k bg is the light extinction coefficient due to absorption by water, coloured dissolved organic matter and suspended particulate material and k is the light absorption coefficient of phytoplankton. The loss rate m of phytoplankton is assumed to depend on salinity. For freshwater phytoplankton, saltier water implies higher loss rate because of osmotic stress. In this study, m is modelled as In this expression, m L and m 0 are the loss rate of phytoplankton in fresh and salt water, respectively, s c is the salinity where m = (m 0 + m L )/2 and s δ defines the salinity 0 1 1 Fig. 3 The shape of loss rate m normalised by loss rate m L of phytoplankton in fresh water. The horizontal axis shows normalised salinity, in which s 0 is sea water salinity, s c is the salinity where m = (m 0 +m L )/2 and s δ defines the salinity scale over which m varies scale over which m varies. Input parameters are m L , m 0 , s c and s δ . An analytical shape of Eq. 17 is plotted in Fig. 3.
The boundary conditions at the water surface, the bottom and the estuary mouth are set as At the surface and the bottom, zero-flux conditions are applied to P and N Eq. 18. At the seaward boundary, zero diffusive flux is assumed for both P and N Eq. 19. At the riverine boundary, the depth-averaged density P L = (1/H ) H 0 P | x=L (z) dz is imposed, and the vertical distribution of P | x=L (z) obeys vP − κ v ∂P /∂z = 0. This condition is motivated by field data in the Columbia River estuary that show freshwater phytoplankton are transported into the estuary from the river compartment. Furthermore, a vertically-constant nutrient concentration is imposed at x = L. In the initial state, P and N at each location along the estuary are set the same with their riverine boundary conditions, respectively.

Model Implementation and Verification
To solve the two nonlinear differential equations in the biological module, central finite difference schemes are employed for spatial discretisation, and a fourth-order Runge-Kutta scheme is used to discretise the time derivatives (Press et al. 1992). This numerical scheme was also used by Liu and de Swart (2015) and it has been verified against data from published studies, yielding good agreement. The hydrodynamic module is solved analytically.

Research Aim 1
Two experiments were carried out to mimic the observed phytoplankton density patterns in the Columbia River estuary during weak stratification and strong stratification, respectively. The effect of vertical stratification (which affects the magnitudes and vertical distributions of vertical eddy viscosity and eddy diffusivity) on P patterns follows from the analysis of the differences between the results of these two cases. In the field data of Roegner et al. (2011), the stratification conditions mainly vary due to variation in tidal mixing during the spring-neap cycle, rather than due to variation in river discharge. These cases are referred to as DF-W and DF-S, i.e. the default cases representing weakly stratified conditions (spring tides) and strongly stratified conditions (neap tides), respectively. Figure 4 displays the map of the Columbia River estuary, in which sampling stations are indicated. The distance Fig. 4 The map of the Columbia River estuary, reproduced from Roegner et al. (2011). The horizontal and vertical axis shows longitude and latitude, respectively. The full circles represent the sampling stations in the south channel, and the inset displays the location of the estuary from sampling station SC16 to SC01 is 30 km. The length of the model domain is chosen as L = 45 km. The reason for selecting this length is that it allows an area for the impact of boundary conditions on P and N dynamics to diminish such that the modelled P patterns in 0 < x < 30 km are determined by internal dynamics.
At the riverine boundary x = 45 km, the same value of the depth-averaged density P L has been used for both DF-W and DF-S. Moreover, for both cases, a vertically constant high value for eddy diffusivity κ v = 10 −2 m 2 s −1 at x = 45 km has been chosen such that P | x=L (z) has little vertical variation. These choices for the riverine boundary condition of P guarantee that P patterns in DF-W and DF-S are not dictated by boundary conditions at the riverine boundary.
The total water depth H =16 m, which is an estimate of the mean water depth in the area 0 < x < 30 km of the Columbia River estuary. A rectangular uniform grid system was used, with 41 × 201 grid nodes in the vertical and along-estuary direction, respectively. The time steps are 5 × 10 −5 day and 2 × 10 −4 day for DF-W and DF-S, respectively.
The values of parameters in the hydrodynamic module are listed in Table 1, including their sources. Here, x c and x δ were determined by fitting Eq. 6 to salinity data of Roegner et al. (2011). Furthermore, motivated by the success of earlier studies in explaining observed features of estuarine circulation (for instance, Hansen and Rattray 1965;Guha and Lawrence 2013), vertical eddy viscosity A v and diffusivity κ v are assumed longitudinally uniform. Accordingly, constant values of u * , z 0 , U T and s are used. Specifically, the same values of friction velocity u * and bottom roughness length z 0 have been employed in DF-W and DF-S. These choices will be further discussed in "Insights from Additional Sensitivity Studies in the Electronic Supplement" section. Parameter U T was estimated as the amplitude of depth-averaged tidal current at x = 11 km using field data of Chawla et al. (2008). The location x = 11 km (that is Station red26 in Chawla et al. 2008) was chosen to estimate U T because its depth is approximately 16 m. Using the field data of Roegner et al. (2011), the tidally-averaged bottom-to-surface salinity difference s was also estimated at x = 11 km. The isohalines have steep slopes under weakly stratified conditions. It is a fair estimate to assume that in this case the loss rate m depends on the depth-averaged salinitiess, i.e. m only varies in the along-estuary direction. On the contrary, isohalines are almost horizontal during strong stratification and m is assumed to vary only in the vertical direction. Specifically, in the latter case, m is a function of the lengthaveraged salinity that was obtained by averaging the field data of salinity along the estuary. For both DF-W and DF-S, s c = 4 psu and s δ = 0.8 psu were prescribed such that m increases from m L to m 0 when the depth-averaged (length-averaged) salinity increases from 3 to 5 psu.
The position wheres is 0.5s 0 km 18 15 The length scale over whichs varies Parameter that is proportional to the values of vertical eddy viscosity and eddy diffusivity at the water surface Tidal-averaged bottom-to-surface salinity difference psu 10 28 κ h | DF Along-estuary turbulent diffusivity m 2 s −1 500 The Light extinction coefficient due to background turbidity m −1 1.6 k| DF (b) Light absorption coefficient of phytoplankton m 2 cell −1 6 × 10 −10 μ max | DF (c) Maximum specific growth rate Loss rate of phytoplankton in fresh water Loss rate of phytoplankton in salt water The salinity scale over which m varies psu 0.8 Sinking  2011), (e) Sarthou et al. (2005) In this study, phytoplankton considered in the biological module represent freshwater diatoms (Lara-Lara et al. 1990) and the depth-averaged density P L at the riverine boundary was set to 10 9 cells m −3 . Furthermore, based on the above-cited study, the nutrient that potentially limits phytoplankton growth in the Columbia River estuary is phosphorus, and the nutrient concentration at the riverine boundary N | x=L = 0.6 mmol m −3 . Finally, the values of other parameters in the biological module are listed in Table 2, which were adopted from other studies that are cited in the footnote of this table.
The following sensitivity experiments were further conducted to achieve the other three research aims. The set-ups of all the experiments are summarised in Table 3. Aim 2: Effect of advection by subtidal currents on P patterns. as DF-W (DF-S), but for the subtidal current completely switched off, i.e. u = 0. as DF-W (DF-S), but for the subtidal current only driven by river discharge, i.e., u = u r . as DF-W (DF-S), but for the subtidal current only driven by along-estuary salinity gradient, i.e., u = u d .
Aim 3: Sensitivity of P patterns to vertical distributions of A v and κ v . as DF-W (DF-S), but for constant vertical shape A z | const . as DF-W (DF-S), but for parabolic vertical shape A z | para .

Research Aim 2
To quantify the impact of subtidal currents u (which consists of river flow u r and density-driven flow u d ) on P patterns, experiments were conducted in which the subtidal currents were completely switched off (u = 0). Because the structure of subtidal currents vary with stratification conditions, additional experiments were conducted in which u = u r or u = u d .

Research Aim 3
To assess the sensitivity of P patterns to the vertical distribution of vertical eddy viscosity A v and eddy diffusivity κ v , experiments were carried out with constant and parabolic vertical shapes of A z , which are denoted by A z | const and A z | para , respectively. Note that a change in A z affects A v through Eq. 8 and κ v through Eq. 12. Compared to the default vertical shape that depends on stratification (denoted by A z | DF ), the former two shapes assume that the height of the turbulent bottom boundary layer equals the total water depth. Furthermore, the depth-averaged values of A z | const and A z | para are kept the same as that of A z | DF . Moreover, the value of A S for A z | para is identical to that of A z | DF , as is illustrated in Fig. 2a.

Research Aim 4
To assess the influence of salinity-dependent loss on P patterns, experiments were conducted with the value of the loss rate m 0 of phytoplankton in salt water halved and doubled, respectively, with respect to that of the default cases (m 0 | DF = 0.5 day −1 ). The former case mimics the behaviour of phytoplankton that are tolerant to the change of salinity, and the latter accounts for the high osmotic stress that freshwater species experience in salt water. Besides, experiments were conducted with a loss rate m| const that is constant in space. Here, m| const = 0.4 day −1 is an approximation of the domain-averaged values of m in both DF-W and DF-S. A spatially constant m| const mimics the case that the osmotic stress of freshwater phytoplankton is roughly accounted for by increasing the loss rate in the whole domain rather than by parameterising it as a function of salinity.
Finally, to test whether model results are qualitatively affected by the choices of boundary conditions in the biological module or by values of other model parameters (for instance, along-estuary turbulent diffusivity κ h , nutrient input N| x=L , maximum specific growth rate μ max , halfsaturation constant H N for nutrient-limited growth, etc.), additional sensitivity experiments were carried out.

Analysis of Output
To facilitate comparison between model results and field data, spatial distributions of model results are shown in the area 0 ≤ x ≤ 30 km. Furthermore, two parameters are defined to quantify the overall vertical and along-estuary gradient of P in the area 0 ≤ x ≤ 30 km: Here, P u max | x=0 and P l min | x=0 are the maximum P in the upper layer and minimum P in the lower layer at the estuary mouth, respectively, and P u max | x=30 km is the density at x = 30 km at the same depth as P u max | x=0 . The definitions of the upper and lower layer are given in "Results" section. The above parameterisations are motivated by field data of Roegner et al. (2011). Furthermore, φ z , being the fraction of the relative change of P between the upper and lower layer, is analogous to the salinity stratification indicator proposed by Guha and Lawrence (2013). Parameter φ x is the fraction of the relative change of P along the estuary over a distance of 30 km. Higher values of φ z (φ x ) indicate larger vertical (along-estuary) gradients of P . Comparison among the results of sensitivity experiments in this study are made in terms of φ z and φ x .

Results
In this section, first the output of the default cases DF-W and DF-S, for the weakly and strongly stratified conditions, respectively, are presented in "Default Cases Under Weakly and Highly Stratified Conditions" section. Next, results of sensitivity experiments listed in Table 3 are presented in "Sensitivity Experiments" section. For all default cases and sensitivity studies, equilibrium states of P are reached after several days starting from the given initial condition (the time evolution of P is not shown here). Results presented of sensitivity studies concern the values of φ z and φ x at equilibrium. Figure 5 contains physical conditions estimated for DF-W and DF-S. The observed and parameterised along-estuary profiles of depth-averaged salinitys are shown in Fig. 5a, b for DF-W and DF-S, respectively. Clearly, the overall behaviour of the observeds is captured by the parameterisation ofs used in this study. In case DF-W, the modelled Colour-contour plots of width-averaged along-estuary subtidal current for DF-W and DF-S, respectively salt intrusion length (x c + x δ ) = 26 km, and in case DF-S (x c +x δ ) = 63 km. These estimates for tidally-averaged salt intrusion length are consistent with those observed by Jay and Smith (1990). Interestingly, data observed by Roegner et al. (2011) show an oscillation in the salinity pattern during high stratification. The potential reason for the presence of this oscillation will be discussed in "Discussion" section. The modelled vertical profiles of eddy viscosity A v and diffusivity κ v are displayed in Fig. 5c. The heights of turbulent boundary layer h b are 8.6 and 3.6 m in DF-W and DF-S, respectively. The maximum of A v for DF-W is 1.4 × 10 −2 m 2 s −1 , which is six times larger than that for DF-S. Similar to A v , κ v for DF-W has a significantly larger magnitude than that for DF-S. Contour plots of the width-averaged along-estuary subtidal current u for DF-W and DF-S are shown in Fig. 5d, e, respectively. Clearly, for DF-W, maximum velocities occur at x = 18 km, where the along-estuary salinity gradient ds/dx reaches its maximum. For DF-S, the magnitude of u is large throughout the domain. Both magnitudes and vertical profiles of the modelled u agree well with those of the observed subtidal current (Chawla et al. 2008). Fig. 6 a,

Default Cases Under Weakly and Highly Stratified Conditions
The spatial distributions of phytoplankton density P at equilibrium for DF-W and DF-S are shown in Fig. 6a, b, respectively. For DF-W, the spatial pattern of P is characterised by little vertical variation and by a gradual longitudinal decline of density from about 10 9 cells m −3 at the riverine boundary to about 0.24×10 9 cells m −3 at the estuary mouth. In contrast, the spatial distribution of P for DF-S shows a two-layer vertical structure. Within a thickness of 6 m below the water surface, the values of P are significantly higher compared to those in the lower water column, and the maxima of P are located with in 1 < z < 2 m. Moreover, P shows less variability in the along-estuary direction for DF-S.
The modelled spatial distributions of P described above capture the main features of those observed by Roegner et al. (2011) (displayed in Fig. 6c, d). Specifically, in the field data, high values of P also occur in the region 0 < z < 6 m. In this study, "the upper layer" refers to 0 ≤ z ≤ 6 m, whereas "the lower layer" is 6 ≤ z ≤ 16 m. Locations of P u max | x=0 , P l min | x=0 and P u max | x=30 km are marked in Fig. 6a, b. Values of φ z and φ x are calculated using Eqs. 20 and 21. Here, φ z = 0.01, φ x = 0.64 in DF-W, and φ z = 0.67, The values of the limitation index (which is defined in "Biological Module" section) suggest that the growth of phytoplankton for both DF-W and DF-S is mostly lightlimited (result not shown here). Nutrient limitation occurs only within 0 < z < 1 m where light intensity is high. The results of nutrient dynamics for DF-W and DF-S and those of sensitivity experiments concerning parameters directly related to nutrient dynamics are presented and discussed in the Electronic Supplement.

Advection by Subtidal Currents
The values of φ z and φ x for the cases in which different subtidal current components are switched off are plotted in Fig. 7a. Under weakly stratified conditions, completely switching off advection of phytoplankton by the subtidal current u results in a negligible decrease in φ z but a noticeable increase of φ x (from 0.64 to 0.88). Furthermore, when the subtidal current only consists of river flow, i.e. u = u r , φ x = 0.68 is close to that of DF-W. In contrast, if only the density-driven flow component is accounted for, that is u = u d , φ x = 0.83 is evidently larger than that of DF-W.
Under strongly stratified conditions, switching off both u r and u d leads to large changes in both φ z and φ x . Specifically, φ z decreases from 0.67 to 0.11 and φ x increases from 0.2 to 0.79. If u = u r , the value of φ z is still much smaller than that of DF-S and the value of φ x is much larger than that of DF-S. However, when u = u d , values of both φ z and φ x are close to those of DF-S.

Vertical Shape of A z
The values of φ z and φ x obtained for different vertical shapes of A z are shown in Fig. 7b, c. For both weakly and strongly stratified conditions, when the default vertical shape A z | DF is changed to a constant shape A z | const or to a parabolic shape A z | para , values of φ z decrease and those of φ x increase. Furthermore, values of φ z and φ x for A z | const differ marginally from those for A z | para . Moreover, Fig. 7 a Scatter plot of φ x and φ z values for different subtidal current components being switched off under weakly stratified conditions (indicated by open circles) and strongly stratified conditions (full circles). b As (a), but for different vertical shapes of A z . Here, the dashed-lined box indicates the plot area of (c). c As (b), but for weakly stratified conditions only. d As (a), but for different loss rates m of phytoplankton. Note the differences in ranges of horizontal and vertical axes in the panels  (d) the ranges over which φ z and φ x vary are much larger during strong stratification than those during weak stratification.

Osmotic Stress
The values of φ z and φ x for different values of loss rate m 0 of phytoplankton in salt water and for a spatially-constant loss rate m| const are plotted in Fig. 7d. Under weakly stratified conditions, as m 0 increases, the value of φ x increases from 0.39 to 0.85, whereas that of φ z hardly changes. Furthermore, when the salinity-dependent loss rate in DF-W is replaced with m| const , the value of P u max | x=0 (defined in "Analysis of Output" section) hardly change, whereas the value of P u max | x=30 km slightly decreases (results not presented). Hence, according to Eq. 21, φ x for m| const is smaller than that of DF-W.
During strong stratification, as m 0 increases, the value of φ z increases from 0.48 to 0.76, and that of φ x slightly increases from 0.16 to 0.21. When m| const is used, φ z = 0.57 and φ x = 0.25 for m| const .

Discussion
To understand the mechanisms underlying phytoplankton dynamics under different stratification conditions, at first, the observed P patterns in the Columbia River estuary are explained by quantifying the contribution of each factor to the local change of P ("Factors Contributing to Local Change of P " section) for cases DF-W and DF-S. Next, the results of sensitivity studies are further analysed in "Sensitivity Experiments" section. Finally, results of additional experiments (that are described in the Electronic Supplement) and the model limitations are discussed in "Insights from Additional Sensitivity Studies in the Electronic Supplement" and "Limitations" sections, respectively.

Factors Contributing to Local Change of P
In the equation governing phytoplankton dynamics (13), the left-hand side term ∂P /∂t is the accumulation rate of phyto-plankton density at a fixed position. The latter is a summation of contributions related to the net growth (viz., gross growth and loss), sinking, along-estuary and vertical advection, along-estuary and vertical mixing. Figure 8 contains the spatial distributions of these terms in the equilibrium state, when ∂P /∂t = 0.

Case DF-W
Due to light limitation, the net specific growth rate (μ − m) (results not shown) of phytoplankton is positive within a thickness of 2 m below the water surface. As a results, positive values of net growth (μ−m)P only occur in a thin layer close to the water surface, as is shown in Fig. 8a. If the vertical mixing is strong, phytoplankton cells produced in the euphotic zone are rapidly moved downward to the aphotic zone (Cloern 1996). This statement is confirmed by Fig. 8c, in which the vertical turbulent diffusion term persistently has negative values in the upper layer. As a result, P is vertically almost uniform (see Fig. 6a). Moreover, the strong vertical mixing results in negligible impact of the sinking of phytoplankton on the P pattern, as is illustrated by Fig. 8e.
High intensity of vertical turbulent mixing, together with light limitation on phytoplankton growth, results in an unproductive estuary (in which the domain integrated value of (μ − m)P is negative) regarding to pelagic phytoplankton production, as was concluded from field studies (Cloern 1987;Lara-Lara et al. 1990;Ragueneau et al. 1996). Consequently, phytoplankton observed in such estuaries are mostly imported from the river compartment rather than are in situ produced (Frey et al. 1983). This point is proven by the fact that the along-estuary advection term (see Fig. 8g) is positive in the upper layer and it generally compensates for the vertical mixing term. Since the depth-averaged (μ−m)P is negative, P gradually decreases as phytoplankton are advected seawards, as is seen in both model results and field data (Lara-Lara et al. 1990;Roegner et al. 2011).
The along-estuary turbulent diffusion term (Fig. 8i) only has large values in the vicinity of the estuary mouth, where the along-estuary advection term becomes zero due to imposed zero diffusive flux boundary condition (Eq. 19), and in the area 20 < x < 25 km, in which the increased loss rate causes increased ∂P /∂x. Moreover, the magnitude of longitudinal turbulent diffusion is much smaller than that of along-estuary advection. Thus, the along-estuary advection of P by subtidal currents has a much larger impact on P patterns than the longitudinal turbulent diffusion of P . The vertical advection term (see Fig. 8k) has a much smaller magnitude compared to that of the vertical turbulent mixing term and it only has noticeable impact in the area where vertical velocity w is large.

Case DF-S
Under strongly stratified conditions, the positive net growth of P is still restricted to a thin layer below the water surface (see Fig. 8b) although the loss rate of freshwater species is low in the surface fresh water. As shown in Fig. 8d, a weak vertical turbulent mixing is not efficient in reducing P in the surface layer. Therefore, high values of P are possible in the upper layer, which is consistent with field observations (such as Ragueneau et al. 1996;Roegner et al. 2011). Note that Huisman and Sommeijer (2002) argued that bloom development requires a minimal vertical turbulent mixing in the euphotic zone to prevent phytoplankton to sink downward into the aphotic zone. The vertical mixing in the upper layer in case DF-S is weak such that sinking of phytoplankton effectively removes phytoplankton from the surface productive layer, as is demonstrated in Fig. 8f.
Because of light limitation, weak vertical turbulent mixing and sinking of phytoplankton, during strong stratification, the prototype estuary is unproductive. Hence, phytoplankton abundance in the estuary domain again relies on the seaward advection of freshwater species by subtidal current, as is demonstrated in Fig. 8h. Furthermore, P slightly decreases towards the estuary mouth as phytoplankton gradually sink out of the upper layer.
The along-estuary turbulent diffusion term (Fig. 8j) is smaller than the longitudinal advection term in the most area of the domain. The former term only show high values in a narrow area close to the estuary mouth, where zero diffusive flux of P is imposed, and in the area x > 25 km, where ∂P /∂x is relatively high due to the constant phytoplankton density assumed at the riverine boundary (see "Biological Module" section). High values of the vertical advection term only occur in the areas where w is large, as is shown in Fig. 8l.

Advection by Subtidal Currents
The importance of the subtidal currents in forming specific spatial distributions of P has been revealed in the previous section, and has been illustrated by the results of sensitivity experiments (see "Advection by Subtidal Currents" section) where different subtidal current components are switched off. During both weak and strong stratification, when the subtidal current u is completely switched off, the seaward transport of phytoplankton from the riverine source is governed only by the along-estuary turbulent diffusive process. Consequently, it takes significantly longer time for the riverine phytoplankton to be moved to the estuary mouth, and Fig. 8 Results shown are the spatial distributions of the terms on the right-hand side of Eq. 13 in the equilibrium state, at which the accumulation rate ∂P ∂t = 0. Left panels are for a the net growth (μ − m)P , c vertical turbulent mixing ∂ ∂z (κ v ∂P ∂z ), e sinking −v ∂P ∂z , g along-estuary advection −u ∂P ∂x in case DF-W. Right panels are as left panels, but for case DF-S. Left panels are for i along-estuary turbulent diffusion ∂ ∂x (κ h ∂P ∂x ) and k vertical advection −w ∂P ∂z in case DF-W. Right panels are as left panels, but for case DF-S for P to be subject to loss and/or sinking processes. Hence, a large along-estuary gradient of P occurs when u = 0, as is indicated by Fig. 7a. These results agree well with Lucas et al. (2009), who used a conceptual model to reveal that in unproductive estuaries, an increase in the transport time will cause the domain integrated loss to increase. Note that under strongly stratified conditions, even if the subtidal currents are completely switched off, a noticeable vertical gradient (φ z = 0.11) of P exists.
Under weakly stratified conditions, P has little variation in the vertical direction because of strong vertical mixing. Since the vertical mean of the density-driven flow u d is zero, it thus hardly contributes to net longitudinal advective transport of phytoplankton. As a result, the values of φ z and φ x for u = u r (i.e., only driven by river flow) are close to those for DF-W (u = u r + u d ).
During strong stratification, P has a marked vertical structure, which implies that contribution of the alongestuary advection of P is sensitive to the vertical distribution of u. In the upper layer, u attains high velocities (maximum 0.62 m s −1 ), which advects phytoplankton quickly seaward. Moreover, the amplitude of u d (∼ 0.51 m s −1 ) is much larger than that of u r . Thus, u d is more efficient than u r in advecting high P to the estuary mouth, and the values of φ z and φ x for u = u d are close to those for DF-S. Besides, in the lower layer u d advects low P landward, which, together with weak vertical mixing, contributes to maintaining vertical gradients in P .

Vertical Shape of A z
Under both weakly and strongly stratified conditions, when the vertical shape of vertical eddy viscosity A v and eddy diffusivity κ v is changed from the default shape A z | DF to either a constant shape A z | const or a parabolic shape A z | para , the values of A v and κ v in the upper layer significantly increase. As a result, in the latter two cases, the negative contribution of the vertical turbulent diffusion (as shown in Fig. 7) to the accumulation rate of P is amplified and φ z therefore decreases. During weak stratification, the values of κ v in the upper layer for A z | DF are sufficiently high such that P is almost vertically uniform and φ z is very low, as shown in Figs. 6a and 7c. Thus, a further increase of κ v in the upper layer causes little decrease in the values of φ z under weakly stratified conditions.
The amplitudes of subtidal current u for the cases A z | const and A z | para are substantially smaller than those for A z | DF . As is discussed in the previous subsection, a reduced magnitude of u leads to an increased along-estuary gradient of P . Moreover, under highly stratified conditions, P largely decreases in the upper layer when A z | DF is changed to A z | const or A z | para . Hence, in the latter two cases, phytoplankton arriving at the estuary mouth are further reduced in density and φ x therefore noticeably increases.
Note that the values of φ z and φ x for A z | const under weakly stratified conditions are largely distinct from those for strongly stratified conditions. This further demonstrates that the differences in the P patterns between weakly and strongly stratified conditions are mainly caused by the difference in the intensity of turbulence (measured by depthaveraged values of A v and κ v ) rather than by the difference in the vertical distributions of A v and κ v .

Osmotic Stress
Under weakly stratified conditions, with a higher loss rate m 0 of phytoplankton in salt water, the net specific growth rate (μ−m) is lower in the salt water area (0 ≤ x < 27 km). Hence, P decreases faster towards the estuary mouth and φ x therefore increases with increasing m 0 , as is shown in Fig. 7.
Under strongly stratified conditions, an increase in m 0 leads to faster loss of P in the salt water area (3.5 < z ≤ 16 m), where positive growth is not possible due to light limitation. The subsequent increase of the vertical gradient of P further amplifies the contribution of vertical mixing. Hence, P in the upper layer decreases faster towards the estuary mouth and φ x therefore increases. Employing a vertically constant m const implies a decrease of (μ − m) in the upper layer and a increase of (μ − m) in the lower layer. Consequently, the difference in P between the upper and the lower layer decreases. Furthermore, as P in the upper layer decreases faster towards the estuary mouth for lower (μ − m), the value of φ x increases. Note that in the model it is assumed that the mortality rate instantaneously changes once the phytoplankton encounter a different salinity. In reality, time lags might be involved, which is a potential interesting point for future studies. On the other hand, this study showed that sensitivity of P patterns to the formulation of m is rather weak.

Insights from Additional Sensitivity Studies in the Electronic Supplement
Parameters κ h , u * and z 0 in the Turbulence Closure Both field data and theoretical studies have shown that the values of longitudinal turbulent diffusivity κ h significantly vary with along-estuary position (Helder and Ruardij 1982;MacCready 2004). Interestingly, when an along-estuary varying κ h was used in the model, the values of φ z and φ x (Fig. S2a) hardly differ from those of the default cases (in which a spatially constant κ h was used). This is because the along-estuary diffusive transport is much weaker than the longitudinal advective transport induced by subtidal current (see "Factors Contributing to Local Change of P " section).
In the experiments shown earlier, the same value of the tidally-averaged friction velocity u * has been applied to both weakly and strongly stratified conditions instead that it is computed from the bottom shear stress. Godin (1991) pointed out that the bottom friction increases with both tidal and subtidal current. When the magnitude of subtidal current is sufficiently smaller than that of tidal current, the value of u * has often been assumed proportional to the amplitude of tidal current (Burchard and Hetland 2010;Chen and de Swart 2016). However, for stratified estuaries, the magnitude of subtidal flow is often so large that it affects bottom shear stress. The amplitude of tidal current is higher during spring tides, whereas the amplitude of subtidal current is higher during neap tides. As the bottom stress due to combined tidal and subtidal currents is not output of the present model, no conclusive relation between the value of u * and the phase of spring-neap tidal cycle can be established. Additional experiments show that when the value of u * is increased, the value of φ z decreases and that of φ x increases (see Fig. S3a) under both weakly and strongly stratified conditions. This is because a higher value of u * results in larger intensity of turbulence, which amplifies the negative contribution of the vertical turbulent diffusion to the accumulation rate of P . However, choosing different values for u * does not affect the earlier analysis on the mechanisms governing P patterns under either weakly or strongly stratified conditions.
Regarding the bottom roughness length z 0 , it is known to vary during spring-neap cycle (Cheng et al. 1999). However, when the value of z 0 is varied in the model, the values of φ z and φ x (Fig. S4a) hardly change during both weak and strong stratification.

Parameters v, μ max , H N and H I in the Biological Module
If the value of sinking velocity v of phytoplankton is increased, under strongly stratified conditions the value of φ z decreases and that of φ x increases (see Fig. S8). The reason is that an increased v leads to a larger decrease of P in the upper layer, as is illustrated by Fig. 8f and is discussed in "Factors Contributing to Local Change of P " section. During weak stratification, the range over which φ z decreases and that over which φ x increases are small. This is because sinking of phytoplankton has little influence on P patterns during weak stratification, as is revealed by Fig. 8e.
For a larger value of the maximum specific growth rate μ max the value of φ z increases and that of φ x decreases (Fig. S9a) under both weakly and strongly stratified conditions. These changes in P patterns are caused by the increase of the specific growth rate μ of phytoplankton (see Eq. 15). Specifically, P values in the upper layer increases with increasing the net growth (μ − m)P , as is illustrated by Fig. 8a and b. Moreover, an increase of (μ − m)P leads to a slower decrease of P towards the estuary mouth. If the value of half-saturation constant of nutrient-(or light-)-limited growth H N (or H I ) is increased, μ decreases (see Eq. 15). Consequently, when increasing H N (or H I ), P patterns are similar to those for decreasing μ max .
Beside the above-mentioned parameters, results (details in the Electronic Supplement) of sensitivity studies regarding other model parameters and boundary conditions reveal that these parameters and conditions only have some quantitative influence on P patterns. Finally, if the net growth of phytoplankton is completely switched off, that is (μ − m) = 0, during both weak and strong stratification P values are high in the entire domain and they increases towards the bottom (see Fig. S12). Thus, to reproduce the observed P patterns with the present model, it is necessary to explicitly model the net growth of phytoplankton instead of treating phytoplankton as a tracer.

Subtidal Flow Components
Only river-driven flow u r and density-driven flow u d have been accounted for in calculating subtidal current in this model. Other subtidal flow components, such as winddriven flow and flow induced by asymmetric tidal mixing, have been also found to contribute to both amplitude and vertical structure of subtidal current (Jay 1991;Cheng et al. 2013). Considerable variation occurs in both the direction and the magnitude of wind stress in the prototype estuary, that is the Columbia River estuary, as observed by Roegner et al. (2011). Furthermore, modelling subtidal flow induced by asymmetric tidal mixing requires explicitly resolving the tidal variations of the vertical eddy viscosity and the vertical shear stress, which is beyond the scope of the hydrodynamic module used in this study.

Domain Geometry
A rectangular cross-section shape and a flat bottom have been assumed for the whole model domain in this study. One important advantage of these assumptions is that it facilitates the employment of analytical solutions for hydrodynamic processes, which allows for easier control and interpretation of the model results. However, previous studies have found that bottom topography influences physical conditions in estuaries. For instance, the intensity of turbulence is dependent on water depth, which further influences subtidal flow. Furthermore, bottom topography, and also lateral constrictions impact salinity distribution directly (Jay et al. 2015), and indirectly by generating internal waves (Geyer and Smith 1987). The latter could have caused the oscillation in the observed salinity pattern as is shown in Fig. 5b. Moreover, the water depth and width variation also directly influence phytoplankton dynamics (Liu and de Swart 2015). Thus, when studying the detailed distribution of phytoplankton, a more realistic geometry should be considered.

Seasonality in Salinity Distribution and Phytoplankton Species
The default experiments in this study are designed for typical physical and biological conditions in spring season in the Columbia River estuary. However, settings in other seasons are different, for example, salinity distribution highly depends on river discharge (Jay and Smith 1990). Furthermore, in summer, phytoplankton observed in the Columbia River estuary mainly originate from the ocean instead from river (Roegner et al. 2011). Moreover, in late summer, Myrionecta rubra, rather than freshwater diatoms, dominate and they form blooms in the lower reach of the Columbia River estuary (Herfort et al. 2012). Thus, when investigating P patterns on seasonal time scales, the above mentioned aspects should be accounted for.

Conclusions
To gain fundamental insight into the mechanisms underlying phytoplankton dynamics in stratified estuaries, an idealised model was applied to quantify the effect of vertical turbulent exchange processes, subtidal currents and salinitydependent phytoplankton loss on spatial distributions of density P under weakly and strongly stratified conditions. Specifically, this was done by fulfilling four research aims: (1) to examine the effect of vertical stratification (which largely affects the magnitudes and vertical distributions of eddy viscosity A v and eddy diffusivity κ v ) and (2) to quantify the importance of the advection by subtidal currents in determining the vertical and along-estuary gradient of P (which are quantified by parameters φ z and φ x , respectively); to assess the sensitivity of the values of φ z and φ x to (3) the vertical distributions of A v and κ v and (4) the salinity-dependent loss rate of phytoplankton.
The model in this study was able to capture the main features of the observed P patterns in the Columbia River estuary. Specifically, under weakly stratified conditions, the values of P are almost vertically uniform and they decrease towards the estuary mouth, whereas under strongly stratified conditions, high values of P occur only in the upper layer and extend seawards. Analysis of model results showed that during weak stratification, the estuary domain is unproductive because of light limitation and strong vertical turbulent mixing. The riverine freshwater species that are advected into the estuary by subtidal current are rapidly mixed downward to the aphotic zone where phytoplankton experience loss. As a result, P gradually decreases seaward. During strong stratification, sinking of phytoplankton, together with light limitation and weak vertical turbulent mixing in the upper layer, results in the unproductiveness of the estuarine domain. After freshwater phytoplankton enter the estuary, high values of P are allowed in the upper layer since the weak vertical turbulent mixing isolates phytoplankton in the surface euphotic zone. Furthermore, phytoplankton accumulated in the upper layer are quickly advected seawards by the subtidal current before being sufficiently subject to loss and sinking processes.
Sensitivity studies revealed that the the advection of P by subtidal currents is important in obtaining a lower value of φ x (viz., a smaller along-estuary gradient of P ) under weakly stratified conditions and in obtaining a smaller φ x and a larger φ z (viz., a larger vertical gradient of P ) under strongly stratified conditions. During weak stratification, the values of φ z and φ x are insensitive to the vertical distributions of A v and κ v . However, the opposite is true during strong stratification. Thus, under the latter conditions, in order to closely reproduce the features of observed P patterns, it is necessary to account for stratification controlled vertical distributions of A v and κ v . A higher loss rate in salt water, which results from is a higher osmotic stress of freshwater phytoplankton, leads to a larger φ x during weak stratification and in a larger φ z during strong stratification. In all these sensitivity experiments, strong vertical gradients of P (i.e. large values of φ z ) only occur if the intensity of turbulence, which is measured by depth-averaged values of A v and κ v , is weak.
Additional sensitivity studies showed that the choices of the values of input parameters only quantitatively influence the values of φ z and φ x , but do not affect the analysis and conclusions previously presented. Nevertheless, since the model in this study is highly simplified, in the future application of the model to study other problems, several aspects should be taken into consideration, for instance, accounting for additional components in subtidal currents, spatially varying geometry and seasonal variations of salinity distribution and phytoplankton species.