Lagrangian Residence Time in the Bay of Gdańsk, Baltic Sea

The impact of synoptic scale and mesoscale variability on the Lagrangian residence time (LRT) of the surface water in the Bay of Gdansk was investigated using the results from an eddy-resolving model. The computed LRT of 53–60 days was up to four times longer than the estimated flushing time reported by Witek et al. (2003). The highest residence times were those of Puck Bay and near the coast, below 50 m water depth, especially during the winter. These sites also had the highest annual mean in LRT. During the summer, when the level of biological activity is high, the LRT distribution was very heterogeneous and patchy, possibly due to the eddy field and to variable wind forcing. Long-term run tracking of the inflowing water from the Vistula River showed a broad spectrum of tracer distribution. The potential impact of a much higher LRT on the near-coastal nitrogen cycle, coastal filter function and genetic differentiation is discussed, and the consequences for coastal zone management are considered. The incorporation of residence time into the Marine Strategy Framework Directive descriptors would result in an improved, unbiased evaluation of good environmental status.


INTRODUCTION
In the 1970s, the relevant time scales of mixing and flushing were referred to in the literature as "flushing time", "turn-over time", "transit time", "mean transit time", "residence time", "mean residence time, " "age", and "mean age." Definitions were provided by Bolin and Rohde (1973) and generalized by Zimmerman (1976), with some modifications, especially in the distinction between local and integral time scales. The differences in the various time scales were described by Jouon et al. (2006). According to Zimmerman (1976), Figure 1 and the legend illustrate the definition of time scales. Figure 1 demonstrates, that the presence of loop currents or eddies would results in longer time scales. Nonetheless, all of these time scales were constructed based on the assumption of long periodic transport times and well-mixed water. However, this does not apply to near-coastal frontal systems, especially in the presence of river plumes draining into the coastal sea. For example, Otto (1978) showed that the variable water circulation in the Southern Bight of the North Sea can result in significant variations from mean values, as also demonstrated for climatologic time scales by Hainbucher et al. (1986). Thus, there is a clear need for investigations of time scales with respect to synoptic and mesoscale variability.  Zimmerman (1976). The time difference t 1 -t 0 from entering the control volume to the actual position x is the "age" of a single water body. The time difference t 2 -t 1 from the actual position x to the edge of the control volume is the "residence time" of the water body. The sum of both is the "transit time" covering the time span from entering to leaving the control volume. If many water bodies are considered, the different time scales can be computed for each water body. Averaging over all water bodies results in "mean age", "mean residence time" or "mean transit time". Properties such as "turn-over time" (Bolin and Rohde, 1973) or "flushing time" (Zimmerman, 1976) cannot be displayed in this graph because they are defined by total volume divided by total flux. The "influence time" (Delhez et al., 2014) although defined differently, has a similarity to "age". Consequently, the presence of eddies or loop currents would cause a longer pathway through the control volume and would result in longer time scales. A plot of the residence time as function of the starting point results in the LRT as function of space and time. The area averaged LRT is equivalent to the mean residence time.
According to Zimmerman (1976), the residence time of a water body in a control volume is defined as the time required for the water body to leave that volume (see also Figure 1 and the legend). This concept of residence time has been widely used in different geographic areas (Dippner, 1993;Cucco and Umgiesser, 2006;de Brauwere et al., 2011;Daewel and Schrum, 2013). An alternative approach was that of Delhez et al. (2014), who proposed the concept of influence time, defined as the time required to replace the water in the area of interest by renewing water.
A fundamental problem of these definitions of time scales is that all are based on a fixed reservoir, such that mixing and flushing time scales are geometrically influenced and the dynamics of synoptic or mesoscale variability are not considered. To overcome the problem of geometric fixed volumes, Dippner (1994) developed the concept of dynamic flushing time, which assumes a variable volume that depends on the position of a frontal system. In the present work, we apply a similar approach to compute the local Lagrangian residence time (LRT) for the Bay of Gdańsk (BG) and Vistula River (VR).
Bay of Gdańsk is one of the larger bays in the southern Baltic Sea. The orography is characterized by the Peninsula Hel, which protects the western part of the Bay from the direct influence of the Baltic Sea. Dynamically the BG is a Bay, which is driven by variable wind forcing and the river plume front of VR. The influence of tides is negligible and also the impact of wave induced mass transport is of minor importance as in many other bays (Dippner, 1987).
The VR is one of the largest rivers draining into the southern Baltic Sea, discharging 15% of the total amount of nitrogen and 19% of the total amount of phosphorus from land-based sources into the Baltic Sea (HELCOM, 1998). A determination of the climatologic winter mean (January andFebruary 1928-2005) molar N:P ratio of inorganic nutrients in the upper 30 m of the BG water-column based on monitoring data (Feistel et al., 2008) yielded values close to Redfield ratio (Figure 2). By contrast, N:P values in the surface water of the BG, measured on a cruise of the "FS Elisabeth Mann-Borgese" (29.2.-9.3. 2016), determined very high ratios, up to 7.8 times higher than the Redfield ratio of 16 (Figure 3). Further data obtained during the cruise are reported in Bartl et al. (2018Bartl et al. ( , 2019. The high and unbalanced Redfield ratios in the coastal zone may alter nitrate cycling (Lunau et al., 2012). To understand the efficiency of nitrate removal, nutrient retention capacity and the complexity of the coastal filter function in the BG (Thoms et al., 2018;Bartl et al., 2019) requires knowledge of the water residence time. Based on the LOICZ (Land Ocean Interactions in the Coastal Zone) model of Gordon et al. (1996); Witek et al. (2003) computed a flushing time of 15 days in the BG, which they named "water exchange time" or "water residence time".
The disadvantage of the method used by many authors is that a fixed volume and a constant transport rate are assumed. This approach does not consider the dynamics of eddies and river plume fronts. While under normal conditions the VR plume propagates eastward, parallel to the coast, during easterly or northeasterly winds it detaches from the coast and meanders in the center of the BG. The resulting eddies have a seasonal varying radius of deformation of 850 m during winter, 4500 m during spring and 6500 m during summer (Voss et al., 2005;Bartl et al., 2018) and a spin-down time of approximately 4.5 days (Voss et al., 2005). These non-linear processes can be taken into account if the residence time is computed using an eddy-resolving model within a Lagrangian framework. Dippner (1994) showed that the dynamic flushing rates in a Lagrangian coordinate system are much lower than geometric flushing rates. Therefore, we hypothesized that the LRT in the BG is much longer than the estimated 15 days (Witek et al., 2003).

The Circulation Model
The circulation model (Holtermann et al., 2014) is an eddy-resolving version of the General Estuarine Transport Model (Burchard and Bolding, 2002). It computes the threedimensional, non-linear, baroclinic primitive equations in a staggered C-grid (Arakawa and Lamb, 1977), assuming hydrostatic and Boussinesq approximations. The model domain, which covers the Central Baltic Sea with very high resolution, is one-way nested in a coarser resolution model that simulates the entire Baltic Sea. Figure 4 shows the outer (Baltic Sea) and inner (Central Baltic Sea) model domains along with the topography and location of Gdańsk Bay in the model. The adaptive scheme used in the vertical direction follows the evolution of shear and stratification in the interior of the water column (Hofmeister et al., 2010). Advection was discretized using a frontal-conserving total-variation-diminishing scheme, following the directional split approach of Pietrzak (1998). Temporal discretization is based on a mode splitting of barotropic and baroclinic pressure gradients (Backhaus, 1985). Vertical mixing is computed with a second-order k-ε closure scheme that is extended by two prognostic equations for the turbulent kinetic energy k and the dissipation rate ε (Umlauf and Burchard, 2005), and horizontal mixing using the scheme of Smagorinsky (1963). Further details are given in Holtermann et al. (2014).
To account for riverine freshwater fluxes, the daily discharges of the Baltic Sea rivers were extracted from the Global Runoff Data Centre database. For the Baltic Sea, only 60% could be covered with direct observations (including the VR). For the remaining rivers, individual climatologic annual cycles were constructed. Afterward, the total annual discharge of all Baltic Sea rivers was computed and compared with the Baltic Marine Environment Protection Commission-Helsinki Commission (HELCOM) estimates (HELCOM, 1986). Then, the climatologically rivers were rescaled such that the sum of observed discharge and scaled climatologically river discharge matched the HELCOM values. The freshwater input of VR, which is the main river in BG, is displayed in Figure 5.
The boundary conditions were adopted from a larger-scale model that covers the entire Baltic Sea (Figure 4) with a setup similar to the one used by Hofmeister et al. (2011). The initial conditions were constructed from multiple Conductivity-Temperature-Depth (CTD) profiles and interpolated to the model grid by applying a distance-weighted algorithm. Data on atmospheric fluxes were obtained from the German Weather Service (DWD). Those data are available with 3 h temporal resolution and were used in order to calculate the air-sea fluxes of momentum, heat and humidity based on the bulk formulas of Kondo (1975).
The topography was linearly interpolated from a 1 NM grid (Büchmann et al., 2011). The model has spherical coordinates with a resolution of 0.00914 • (501-597 m) and 0.0054 • (596 m) for zonal and meridional directions, respectively. The time step is 3 s for the barotropic and 60 s for the baroclinic mode. The results were stored every 3 h. Although the complete model run covers the period 2007-2017 a detailed validation was performed for the years 2007-2009 in order to coincide with the Baltic Sea Tracer Release Experiment (BATRE) (Holtermann et al., 2014).
The model was validated against observed sea level elevations, temperature and salinity profiles at different locations in the Central Baltic Sea and fit well to the observations FIGURE 5 | Runoff of VR for 2009 from the Global Runoff Data Center. (Holtermann et al., 2014). For the presented analysis, we use the computed surface velocities of 2009.

The Lagrangian Model
According to Maier-Reimer (1979), the local LRT is computed within a Lagrangian framework. If a water body is interpreted not as a continuum, but as a finite set of water particles transported by the Eulerian velocity field, then each particle is characterized by a co-ordinate x and an age a. Accordingly, the numerical problem of solving an advection-diffusion equation is reduced to the simple problem of solving a set of ordinary differential equations describing the change in coordinates and properties for each particle. Hence, transport, diffusion and age are calculated as described in Eq. 1: where i is the i-th particle, t is the time, t is the time step, P is a "band-width" of turbulent fluctuation and µ is a random number uniformly distributed between 0 and 1. The advection velocity u L is obtained from the Eulerian velocity field of the circulation model by transformation onto Lagrangian coordinates (Maier-Reimer, 1988;Dippner, 1990). A reflecting boundary condition for the random walk is used at the closed boundary. From a numerical perspective, the main disadvantage of the Eulerian system is the advective part of the transport equation, which yields the well-known problem of numerical diffusivity, especially in the presence of strong gradients (Dippner, 2005). However, by transformation onto Lagrangian coordinates, which is the core of the technique, this part simply vanishes. For the diffusion part of the equation, which a priori is assumed to be of minor relevance, a simple Monte Carlo technique is used. Eq. 1 is solved using a 4thorder Runge-Kutta method. In this simulation, the particles are advected at the horizontal surface velocities, excluding a vertical component. Thus, the particles can be interpreted as buoyant particles at the surface. Two different complementary applications were computed to analyze the spatial as well the temporal variability. In the first application, we analyzed the horizontal variability and initialized at each wet grid point of the model a tracer with an age (a = 0d), excluding the Vistula Lagoon and compute the time for all 19,458 tracers required to leave the model domain. This was done with the velocity fields for 2009, beginning January 1st and diffusion (P = 0 m/s) was not considered. This procedure was repeated for all other months starting on the 1st of each month (Figures 6, 7). To estimate the influence of diffusion, the runs were repeated using different diffusion velocities. Figure 8 displays both, the time required for all tracers to leave the model domain (T 100 ) and the time at which 95% of the tracers have left the area (T 95 ) for all month and three different diffusion velocities, 0, 0.1, and 0.2 m/s. Averaging over all months yielded an annual mean LRT (Figure 9), with the standard deviation providing information on the spatial variability of the LRT in the BG (Figure 10).
To estimate the LRT of the Vistula plume only and the temporal variability of LRT, in a second application 105,120 tracers were initialized in the vicinity of the VR's entrance.
For 1 year, a tracer was released every 5 min and advected until all tracers had left the domain. The entire computation covered the period from January 1, 2009 to March 28, 2010. With this approach, the LRT is identical to the Lagrangian transit time as shown in Figure 1 (Zimmerman, 1976). The spectra of the LRTs of the tracers and their cumulative spectra were shown in Figure 11. Again to estimate the potential influence of diffusion, the run was repeated using different diffusion velocities. The spectra for diffusion velocities of 0.1 and 0.2 m/s were also shown in Figure 11.

Nutrient Export
According to Nixon et al. (1996), the LRT (in months) can be used to compute specific, residence-time-dependent estuarine quantities using the bulk formulas for the export of total nitrogen (TN) and total phosphorus ( These bulk formulas were derived from regression analyses of relative TN or TP exports (% of total N or P inputs) as a function of mean water residence time (month) (Nixon et al., 1996, Figures 1, 2 therein). Since the regression analyses are based on data from semi-enclosed estuaries, lagoons and freshwater lakes, their application to the open BG is limited and regarded solely as qualitative estimations of the relative amounts of exported and retained nutrients.

RESULTS
The monthly spatial distributions of the LRT showed rather heterogeneous patterns (Figures 6, 7) and a pronounced intra-annual variability. Details of monthly circulation patterns and wind fields were documented in the Supplementary  Information (SI). The shortest LRTs during the considered year were in January and October, which were caused by the dominant southerly to south-westerly winds resulting in a strong flushing. During November, December and February, the longest LRT was in the Puck Bay. These long LRTs were caused by mainly southerly winds from southwest to southeast. The consequence was a stagnant circulation in the Puck Bay affecting long residence times, which can last more than 3 month in November and December (Figures 7, 8). For 95% of the tracers, 147 days beginning in November and 117 days beginning in December were required to leave the bay. The combination of these wind fields with the semi-enclosed geometry of the Puck Bay resulted in the highest mean and standard deviation of the LRT (Figures 9, 10). During these months, the stagnant circulation for longer than 3 months is also influenced by the wind history and the accumulative consequences of the wind forcing (see Supplementary Information). By contrast, the LRTs during spring and summer were patchy and highly heterogeneous. The annual mean residual circulation in the BG computed from all velocity fields was 1.7 cm/s. Wind fields in opposite direction to the mean circulation resulted in a stronger perturbation of the circulation field, e.g., a detachment of the river plume front from the coast, the formation of baroclinic instabilities and the formation of eddies (Voss et al., 2005). The formation of meandering patterns and filaments clearly reflected the influence of the eddy field and the variability of wind forcing on the LRT (see also Supplementary Information).
Assuming that the T 95 value of each month can be considered as a monthly mean residence time (Figure 8), the annual mean residence time, obtained by averaging over all months, was 60 days. If months with extremely low and high values (January and November) were excluded, the averaged residence time was 56 days. Figure 8 also displays the simulations with different diffusion velocities. The colored impulses were the T 95 values and the gray bars the T 100 values for each model run. The T 95 values were identical in January, October, November and December, During the other month of the year the differences were small. The T 100 values showed a lager variability. With respect to most of the T 95 values, the time at which 95% of the tracers have left the area, diffusion seems to play a minor role. A possible explanation is given below.
The spatial pattern of the mean LRT (Figure 9), averaged over all months, could be split into three different areas. (i) north of 54.6 • N, where the LRT was <24 days; (ii) in Puck Bay, where the LRT was longest, >48 days, and (iii) over a small belt in front of the Polish coast, the near-coastal area, where the LRT was between 24 and 36 days. The water depth in the latter two areas was shallower than 50 m. In the graph of the standard deviation of the LRT (Figure 10), two different areas were identified: in the Puck Bay, where the variability in the LRT was highest, and the  western near-coastal area, where the LRT exceeded 36 days. In the rest of the bay the variability in the LRT was <24 days.
Investigating the temporal variability and the fate of VR water, long-term simulations were performed with and without considering diffusion (Figure 8). The spectra as well as the cumulative spectra were identical, indicating that turbulent diffusion might be of minor relevance on this scales. Two larger peaks seen at 7 and 16 days have amplitudes of 4.6% and of 3.7%, respectively, while the minor peak at 29 days had an amplitude of 1.7%. Along with the peaks, the spectra were broad, with a smooth e-folding in its tail and an amplitude of 0.5% until 55 days. A view of the cumulative spectra showed that after 16 days, i.e., close to the estimated flushing time reported by Witek et al. (2003), only 47% of the tracers had left the BG. The LRT of 95% of the tracers was 53 days, which was more than three-fold higher than the value estimated by Witek et al. (2003). No significant changes in the tracer spectra were identified by considering diffusion in the long-term simulation.  The maximum speed in the surface velocity was 1.4 m/s, which is large compared to the maximum diffusion velocity of 0.2 m/s. The narrow grid of the eddy-resolving model allows to resolve small scale processes such as stretching and shear deformation, which contribute generally to diffusion. Therefore, an additional diffusion velocity seems to be not necessary in the model.

DISCUSSION
Our results indicate water residence times up to four times longer than the 15 days flushing time estimated by Witek et al. (2003). This difference was remarkable given that we used only the surface velocities. Using the model of Delhez et al. (2014), the computed influence time for surface water in the BG was in the same range as our results (Bartl, pers. comm.), which suggests that the computations of residence time in a Lagrangian model and influence time in a Eulerian model seem to be equivalent. This similarity will be the subject of a future investigation. The potential impact and consequences of the LRT in the surface water of the BG are discussed in the following.

Spatial Differences of LRT
The residence times in Puck Bay were significantly longer than in the rest of the BG, especially in winter (November to March), although it has no allochthonous nutrient source. Puck Bay is characterized by a shallow water depth, effective vertical mixing and a high population density of different benthic communities (Janas et al., 2004). Its sediments provide a large remineralization potential and play a key role in nutrient turnover (Andersen and Kristensen, 1992). Thus, together with the long surface-water LRT, recycled nutrients from the sediments may strongly impact the productivity in Puck Bay.
The 50-m depth contour in the BG (Figure 3) separates near-coastal sandy sediments from offshore muddy sediments (Thoms et al., 2018). This differentiation also appears in the LRT, which was longer in the coastal area shallower than 50 m water depth than in the northern offshore area. The longer surface-water LRT in the near-coastal area may be conducive to the sedimentation of organic matter deriving from primary production on the near-coastal seafloor. This would fuel benthic nutrient turnover, especially nutrient retention via fluxes (Thoms et al., 2018) and nitrification (Bartl et al., 2019), in turn facilitating permanent nutrient-removal processes such as denitrification and burial. Thus, the surface-water LRT may have an indirect impact on the permanent removal of nutrients from this coastal system, which functions on time scales much longer than the mean water residence time (Asmala et al., 2017). By contrast, the bottom water residence time may have a more direct influence than the surface LRT on benthic nutrient retention and permanent nutrient removal. Overall, the water residence time strongly influences how often a nutrient cycles within a coastal system before it is removed or exported. This complex behavior is currently the subject of a detailed model-based study computing the influence time at the surface and bottom of the BG (Bartl, pers. comm.).

Nutrient Cycle and Filter Function
Inputs to the Baltic Sea of TN and TP from rivers, atmospheric deposition and nitrogen fixation have been slowly decreasing since the 1980s and are currently estimated to be 638 kt/year for TN and 28 kt/year for TP. These loads include dissolved inorganic, organic and particulate sources (HELCOM, 2011). Inorganic nutrients in the form of nitrate and phosphate contribute most to the TN (45% dissolved inorganic nitrogen, DIN) and TP (51% dissolved inorganic phosphorus, DIP) loads of the VR (Pastuszak and Witek, 2012).
However, from 1989 to 2006 the DIN concentrations in coastal waters increased (HELCOM, 2009) but since then seem to have stagnated, given the recent decreases in the loads from major rivers (HELCOM, 2011). DIN and DIP loads to the Baltic Proper are roughly 200 kt N and 7 kt P (taken from Figure 3.8, HELCOM, 2009), with a molar N:P ratio of roughly 13, which is in agreement with the long-term climatic mean (Figure 2). Whether this ratio systematically changes in near-coastal areas due to microbial or chemical processes is unknown, but the significantly lower N:P ratio in the Baltic Proper is well documented (Voss et al., 2006). The residence time is the most important factor regulating nutrient cycling (Nixon et al., 1996). Nutrient retention, i.e., the biological cycling of nutrients via their incorporation into biomass and remineralization, slows nutrient transport from the land to the sea whereas nutrient removal via denitrification or burial permanently reduces the inventory (Asmala et al., 2017). Nutrient retention facilitates a nutrient residence time that is longer than the mean water residence time (Asmala et al., 2017), which in turn allows permanent nutrient removal to mitigate riverine nutrient loads in the coastal zone.
In the surface water, the active processes supporting nitrogen retention are the assimilation of inorganic nutrients into biomass and regeneration. Estuarine conditions similar to those of the BG are found in Pomeranian Bay, off the Oder River Lagoon. In a study of that bay by Korth et al. (2013), nitrate assimilation rates of up to 235 nmol N L −1 h −1 were determined in March, suggesting daily rates of 2.82 µmol L −1 d −1 in the surface waters (1 m water depth, constant rates for 12 daylight h). Assuming similar rates in the near-coastal area of the BG, nitrate assimilation would consume an appreciable amount of 17-152 µmol L −1 during a LRT of 6-54 days. The nitrate concentration of the VR in March was ∼186 µmol L −1 (see Figure 9.6 in Pastuszak and Witek, 2012), resulting in a nitrate input of 31-281 µmol L −1 within 6-54 days, assuming a discharge rate of 1600 m 3 s −1 (see Figure 9.5 in Pastuszak and Witek, 2012) into an estuarine surface water volume of 4.94 km 3 (BG surface area × 1 m water depth). Thus ∼54% of riverine nitrate may be incorporated into biomass and thereby retained. A similar quantity (57%) can be estimated using the equation of Nixon et al. (1996), according to which approximately half of the VR's input is retained and half is exported. The large temporal variability in the LRT and in nitrate uptake rates (Korth et al., 2013) in spring points to a large potential variability in nitrogen retention. Since no nitrate uptake rates are available from other seasons, we can only speculate here. In summer, when primary production is highest (Witek et al., 1999) concurrent with long LRT (Figures 6, 7), all riverine DIN is likely incorporated into biomass. In autumn and winter, however, LRT is short except in the Puck Bay (Figure 7) and primary production (Witek et al., 1999) is low while riverine N loads can be high resulting in pronounced export of N. Consequently, the timing of the surfacewater residence time together with nutrient discharges from the VR and nutrient uptake rates determines the efficiency of nutrient retention in the BG.
Half of the plankton biomass resulting from nutrient assimilation was degraded in the surface waters and remineralized to inorganic nutrients within 30 days (Garber, 1984). Since the LRT is similar throughout the spring (March-May) and riverine nutrients are continuously delivered, a bloom could theoretically develop and persist during the spring months. In addition, assimilation and degradation, including ammonification and nitrification, in the coastal surface waters for the entire spring period can be assumed. Indeed, mean nitrification rates of 59 nmol L −1 d −1 were detected during a spring bloom in the surface waters of the BG, and further offshore they were much higher 258 nmol L −1 d −1 (Bartl et al., 2018). This observation is consistent with the high N:P ratio in the nutrient load, which is ∼140 in winter and spring but only ∼10 in summer (Voss et al., 2006). The highest Redfield ratios measured in the surface waters of the BG during winter were in Puck Bay and the near-coastal area shallower than 50 m water depth (Figure 3). These ratios are far above the ratio of ∼16:1 needed to meet the metabolic requirements of plankton. An experiment with estuarine water in which a gradient in N:P ratios was applied demonstrated that nitrification rates are significantly enhanced under a surplus of nitrogen (Lunau et al., 2012). A field study along an estuarine gradient also reported a significant increase in nitrification rates at sites where nitrate levels were highest (Dähnke et al., 2008). These two studies independently support enhanced nutrient retention under high N-loads and a sufficient residence time.

Impact of LRT on Genetic Differentiation
The results of this study also help to explain the strong genetic differentiation among phytoplankton populations in the BG. Previous studies of the genetic differentiation of two bloomforming species in the Baltic Sea, the diatom Skeletonema marinoi and the toxic dinoflagellate Alexandrium ostenfeldii, showed that they were genetically distinct from other local Baltic populations (Tahvanainen et al., 2012;Sjöqvist et al., 2015). Our study showed that the bloom periods for S. marinoi and A. ostenfeldii (February-April and August-September, respectively) were characterized by long residence times, which would prevent transport of the blooms into the main basins and, for LRTs >30 days, facilitate anchoring of the blooming populations in the bay. At the end of their bloom periods, the two species enter a resting stage in which the cells settle to the sea floor, where they form benthic seed banks (Kremp et al., 2009;Godhe and Härnstöm, 2010). Resting stages that accumulate in the inner parts of Puck Bay not only inoculate new local blooms but also facilitate the retention and adaptation of local populations (Sundqvist et al., 2018), which in the case of toxic, bloom-forming A. ostenfeldii will adversely affect both co-occurring biota and public health (Setälä et al., 2014).

Consequences for the Integrated Management of Coastal Zones
Inputs of nitrogen and phosphorus from land into marine coastal systems are responsible for enhanced primary production and therefore eutrophication (Nixon, 1995). According to the Marine Strategy Framework Directive (MSFD, Anonymous, 2008): "Eutrophication is a process driven by enrichment of water by nutrients, especially compounds of nitrogen and/or phosphorus, leading to: increased growth, primary production and biomass of algae; changes in balance of organisms, and water quality degradation. The consequences of eutrophication are undesirable if they appreciably degrade ecosystem health and/or the sustainable provision of goods and services" (Ferreira et al., 2010). The consequences of large-scale eutrophication in the Baltic Sea include reduced water transparency, oxygen depletion, and changes in the marine food web (Fonselius, 1972;Bonsdorff et al., 1997;Karlson et al., 2002;Conley et al., 2011;Voss et al., 2011;Gustafsson et al., 2012).
In addition to being the most important factor regulating nutrient cycling, residence time determines the nutrient retention time, total nutrient export (nitrogen and phosphorus), and the total removal (buried and denitrification) of nitrogen (Nixon et al., 1996). In addition, the LRT is of major relevance for the persistence of harmful contaminants or pollutants as well as toxic algal blooms and influences the efficiency of coastal filter function. All these aspects must be considered in assessments of good environmental status using the MSFD descriptors "biological diversity, " "food web, " "eutrophication" and "seafloor integrity." However, the Task Group Reports on "Biological diversity" (Cochrane et al., 2010), "Food webs" (Rogers et al., 2010), "Eutrophication" (Ferreira et al., 2010) and "Seafloor integrity" (Rice et al., 2010) do not explicitly consider residence time as an important variable influencing and regulating processes in near-coastal marine ecosystems. MSFD descriptor 5 (eutrophication) refers to the adverse effects of eutrophication as including "losses in biodiversity, ecosystem degradation, harmful algae blooms and oxygen deficiency in bottom waters" (Ferreira et al., 2010). Nonetheless, because eutrophication is strongly linked with the other MSFD descriptors mentioned above, via strong bentho-pelagic coupling, residence time is an important driver for all of them. The incorporation of residence time into these MSFD descriptors would result in an improved, unbiased evaluation of good environmental status.
These considerations are also of major relevance in the context of nutrient reduction efforts. Conley (2000) suggested that a reduction in phosphorus loads in spring will contribute to reducing both oxygen deficiency and fast-growing macrophytes, whereas a reduction in nitrogen loads would reduce summer chlorophyll concentrations and improve the conditions for benthic vegetation. A much larger residence time during summer would result in a prolongation of high-levels algal biomass production. To avoid undesirable algal blooms, reductions of both nitrogen and phosphorus are urgently required for the conservation of near-coastal areas and for the achievement of a good ecological status.

DATA AVAILABILITY STATEMENT
The model results supporting the conclusion of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
JD, EC, and PH made the computations. IB, AK, FT, and MV contributed to the text.

FUNDING
EC was the recipient of support for the project T2 (Energy Budget of the Ocean Surface Mixed Layer) of the Collaborative Research Centre TRR 181 on Energy Transfer in Atmosphere and Ocean, funded by the German Research Foundationproject number 274762653. IB and FT were funded through the BONUS COCOA project, by BMBF grant number FKZ 03F063A. PH was financed by grant HO5891/1-1, from the German Research Foundation. Additional support was received by the Open Access Fund of the Leibniz Association.

ACKNOWLEDGMENTS
The runoff data from Vistula River has been provided by the Global Runoff Data Centre, Koblenz, which is greatly acknowledged.