Shipborne nutrient dynamics and impact on the eutrophication in the Baltic Sea

• Relative impact of shipborne nutrients does not exceed 10% locally. • The North-Eastern and South-Western Baltic Sea are the most impacted. • Nitrogen fixing cyanobacteria biomass decreases due to the shipborne nitrogen


Introduction
In natural waters, phytoplankton growth is usually limited by physical and biogeochemical factors (Hecky and Kilham, 1988;Chislock et al., 2013). Human-induced change in these limiting factors may have a considerable effect on the entire ecosystem where phytoplankton usually plays the key role (Schelske, 1975;Hecky and Kilham, 1988;Chislock et al., 2013). In coastal waters the limiting element is most often a nutrient; either nitrogen or phosphorus (Hecky and Kilham, 1988;Granéli et al., 1990). Changes in the input of the limiting nutrient into the particular coastal ecosystem may, in turn, alter phytoplankton biomass and cause other ecosystem alterations (Schelske, 1975;Schindler, 2006;Chislock et al., 2013). Since ecosystems are complex and their responses non-linear, secondary changes, which occur as a consequence, are more difficult to assess in advance (Duarte et al., 2009;Savchuk, 2018).
Maritime transport, among other sources, increases nutrient input into the marine ecosystem via emissions to the air and the concurrent deposition, as well as direct discharge Bartnicki et al., 2011;Herz and Davis, 2002). Emissions from marine transport contribute significantly to air pollution globally, with commercial shipping estimably responsible for around 4-14% of the global nitrogen emissions (Wang et al., 2008;Corbett and Fischbeck, 1997;Corbett and Köhler, 2003;Eyring et al., 2005). Marine pollution from direct sanitary wastewater discharge is an especially pronounced problem with large cruise ships, which represent b1% of the global merchant fleet, while estimably being responsible for 25% of all the waste generated by merchant vessels (Perić et al., 2016;Herz and Davis, 2002). Global maritime transport is experiencing an increasing trend (Eyring et al., 2005;Marmer et al., 2009).
Worldwide, the coastal areas of Western and Southern Europe (North Sea and the Mediterranean), the United States and South-East China have been more thoroughly studied regarding the amount and effect of shiprelated emissions on the environment (Viana et al., 2014;Marmer et al., 2009;Aksoyoglu et al., 2016;Chen et al., 2019;Djambazov and Pericleous, 2015).
In Northern Germany and Denmark that are surrounded by numerous shipping lanes, the contribution of shipping emissions to the atmospheric nitrogen dioxide is around 15% in winter and 25% in summer, whereas in the western entrance of the English Channel the shipderived nitrogen dioxide contribution can be up to 90% . This demonstrates that in the case of lacking additional sources, shipping can be the primary contributor to the atmospheric input of nutrients, providing extra fuel for the already intense eutrophication in the coastal waters of Western Europe. Quantitatively, the contributions from shipping emissions to nitrogen dioxide pollution show a large spatial variability, with maximal contributions in the Mediterranean basin and the North Sea, and the average being 7-24% across Europe (Viana et al., 2014).
A highly industrialized and rapidly transitioning area with dense marine traffic lies in South-East Asia. Modelled annual ship emissions on the most active routes reach up to 10 4 kg N/km 2 /yr (Chen et al., 2019). The ship-induced increase of nitrogen deposition is not only found along the shipping routes but also spread to wide land regions. The highest simulated total N deposition from ship emissions appeared in the coastal areas, which reached 1000 kg N/km 2 /yr. In the sea area, the N deposition flux was ranged from 200 kg N/km 2 /yr to 900 kg N/km 2 /yr. The Arctic Ocean is an area of growing global interest, as increasing navigability enables more intense/frequent shipping activity that will bring about higher input of nutrients and pollutants into a formerly pristine environment (Gong et al., 2018).
Multiple studies have implied that additional input of nutrients has an effect on (marine) ecosystems, especially by increasing primary production and decreasing oxygen concentrations due to the increased production of dead organic matter (Herz and Davis, 2002;Perić et al., 2016;Hagy et al., 2004;Spokes and Jickells, 2005;Troost et al., 2013;Djambazov and Pericleous, 2015). In a study by Troost et al. (2013) additional nitrogen from atmospheric deposition lead to increased primary production in the North Sea. There are temporal and spatial differences in the effect of additional nutrient input, which determines ships as effective proxies to areas lacking other nutrient sources. The offshore areas (N30 km off the coast) are, in general, more affected due to absence of land input and the system being nitrogen limited during long periods in summer (Troost et al., 2013). It is thus mainly the nitrogen-limited phytoplankton species that benefit from the atmospheric deposition, that has been estimated to be accountable for 13.8-15% of primary production in several studies (Troost et al., 2013;Spokes and Jickells, 2005). Atmospheric nitrogen input can, hence, support a significant increase in phytoplankton biomass and will, along with other nitrogen sources, enhance long-term eutrophication effects in the intensely shipped coastal waters of Western Europe and England (Spokes and Jickells, 2005;Djambazov and Pericleous, 2015). Similar results have been obtained or are expected to occur in the actively industrializing South-East Asia region and the increasingly shipped Arctic Sea (Zhao et al., 2015;Gong et al., 2018). Through its effects on primary production, addition of nutrients may also influence the rest of the ecosystem via changes in the carbon equilibrium or shifts in the species composition on higher trophic levels of the food chain (Voss et al., 2001;Van de Waal et al., 2010).
The Baltic Sea has been, similarly to other intense shipping areas, under human pressure for a long time, as well as facing natural challenges from being a northerly inland sea with long residual time, limited water exchange and slow degradation processes in a temperate region (Zillén et al., 2008;Reusch et al., 2018). The anthropogenic nutrient input to the Baltic Sea has been the major cause of eutrophication and the consequent extensive cyanobacteria blooms (Elmgren, 1989;Reusch et al., 2018;Jonson et al., 2015;Aksoyoglu et al., 2016;HELCOM, 2014). Contribution of the shipping sector to the total atmospheric deposition of oxidized nitrogen into the Baltic Sea is driven by the source strength as well as by the meteorological conditions which means that the annual contribution does not only vary with changing shipping emissions, but also with inter-annual variability of the meteorology (Bartnicki et al., 2011). Bartnicki et al. (2011) andJonson et al. (2015) studied atmospheric deposition of nitrogen to the Baltic Sea for the time period between 1995 and 2012 with atmospheric chemistry-transport model EMEP and found contribution from the Baltic Sea and North Sea shipping to the total deposition of oxidized nitrogen to be on average 18 kt N/year which makes a relative contribution of 13-20%. A similar result was obtained by Raudsepp et al. (2013) in the Gulf of Finland, where NOx deposition to the sea from ship exhaust gases was estimated to be about 12% of the annual atmospheric NOx deposition. The total atmospheric nitrogen deposition is estimated to be about 220 kt N/year (HELCOM, 2013a(HELCOM, , 2013bSimpson, 2011) and the waterborne nitrogen input ca. 760 kt N/ year, which makes the atmospheric contribution of nitrogen around 22%. Therefore, the shipping contribution is about 1.25-3.3% of the total nitrogen to the Baltic Sea.
Phosphorus enters the Baltic Sea mainly as waterborne while the atmospheric contribution is calculated to be approximately 5.5% of the total phosphorus input (HELCOM, 2013b). Shipborne phosphorus enters the sea via waste generated onboard the ships such as sewage (black water), food waste and grey water (Wilewska-Bien et al., 2018). Other possible phosphorus sources of significance are bilge water and atmospheric deposition (Wilewska-Bien et al., 2018;Neumann et al., 2018bNeumann et al., , 2018c. The atmospheric contribution of phosphorus from ships is small because marine fuels contain very little (b15 ppm) phosphorus (Wilewska-Bien et al., 2018;ISO 8217, 2012). Compared to the overall phosphorus input to the Baltic Sea, the phosphorus from shipping is estimated to comprise around 0.3% of the total annual input (HELCOM, 2013b;Wilewska-Bien et al., 2018).
Similarly to the North Sea, any additional nutrient deposition in the nitrogen limited offshore areas of the Baltic Sea influences the phytoplankton biomass during spring by changing the interspecies dynamics of different functional groups (Tilman, 1982;Stockner and Shortreed, 1988). Furthermore, increased production of organic matter also increases oxygen demand and may expand the hypoxic areas below the upper mixed layer of the Baltic Sea water column (Zillén et al., 2008;Conley et al., 2002). The share of nutrient input from shipping to the Baltic Sea is small compared to the total input, but its relative importance may become important due to spatio-temporal variance of the different sources and the natural spatio-temporal dynamics. Further, despite being a small share, the shipborne nutrient input contributes to the cumulative nutrient input as one of many small sources.
The objective of this study is to evaluate the response of the Baltic Sea ecosystem to excess nutrients from shipping (SHIP scenario minus NOSHIP scenario). While other studies have either focused on a certain area of the Baltic Sea (Raudsepp et al., 2013;Neumann et al., 2018aNeumann et al., , 2018bNeumann et al., , 2018c, or a selected discharge or nutrient type (Wilewska-Bien et al., 2016, the current study takes into account all shipping-related nutrient sources: atmospheric emissions and the concurrent deposition, as well as direct discharges of different categories. We evaluate the impact of shipborne nutrients on the overall nutrient-phytoplankton-oxygen dynamics across the entire Baltic Sea and determine which processes are responsible for a transformation of the nutrients. We compare the situation under current regulation of ship emissions (SHIP) to a zero shipborne nutrient emission scenario (NOSHIP).

The Baltic Sea
The Baltic Sea is a brackish semi-enclosed sea in northeastern Europe with a surface area of 422,000 km 2 and a volume of 21,205 km 3 (Leppäranta and Myrberg, 2009). The Baltic Sea is characterized by relative shallowness with an average depth of 54 m, while N1/3 of the sea is shallower than 30 m (Fig. 1), giving it a small total water volume relative to the surface area. The maximum depth of 459 m is located in the deep trench of the Landsort Deep. The Gotland Deep with a maximum depth of 250-m in the central Baltic Proper is considered a dynamic deep area with a high significance in shaping hydrographic and biogeochemical fields of the Baltic Sea. The long-term salinity is determined by net precipitation (e.g. Jaagus et al., 2018) and river discharge across the Baltic Sea coast (Hansson et al., 2011) and by the saline water inflows from the North Sea through very narrow and shallow channels in the Danish Straits (BACCII Author Team, 2015). The saline and oxygenated water inflows to the Baltic Sea, especially the Major Baltic Inflows, occur only intermittently (e.g. Mohrholz, 2018).
The surface salinity varies horizontally from~10 near the Danish Straits down to~2 at the northernmost and easternmost subbasins of the Baltic Sea. The halocline, a vertical layer with rapid depthdependent changes of salinity that separates the well-mixed surface layer from the weakly stratified layer below, is located at the depth range of 60-80 m (Matthäus, 1984). The bottom layer salinity below the halocline depth varies from 15 in the south down to 3 in the northern Baltic Sea (Väli et al., 2013). Long-lasting periods of oxygen depletion in the deep layers of the central Baltic accompanied by salinity decline and overall weakening of vertical stratification are referred to as stagnation periods. Extensive stagnation periods occurred in the 1920s/1930s, in the 1950s/1960s and in the 1980s/ beginning of 1990s (BACCII Author Team, 2015).
The upper mixed layer temperature of the Baltic Sea is characterized by a strong seasonal cycle driven by the annual course of solar radiation (Leppäranta and Myrberg, 2008). Maximum water temperatures are reached in July and August and minimums during February, when the Baltic Sea becomes partially frozen. The strongly seasonal sea ice coverage has a vital role in the annual course of physical and ecological conditions. In general, sea ice starts to form in October and may last until June. Depending on the year maximum ice extent could vary in the range of 30,000 km 2 (e.g. in 2015) to 260,000 km 2 (e.g. in 2011). In case of a fully ice-covered Baltic, the maximum ice extent is 422,000 km 2 , which was last observed during the 1940s (Vihma and Haapala, 2009). Temporal trends of sea ice extent could be a valuable indicator of the climate change signal in the Baltic Sea region. It has been estimated that a 1°C increase in the average air temperature results in the decline of ice-covered area in the Baltic Sea by about 45,000 km 2 (Granskog et al., 2006). Seasonal thermocline, developing at the depth range of 10-30 m in spring, is strongest in summer and is eroded in autumn. In autumn and winter the Baltic Sea is thermally mixed down to permanent halocline at the depth of 60-80 m (Matthäus, 1984). The 20-50 m thick cold intermediate layer forms below the upper mixed layer in March and is observed until October within 15-65 m depth (Chubarenko and Stepanova, 2018;Liblik and Lips, 2011). Decrease in maximum ice extent may influence vertical stratification of the Baltic Sea (Hordoir and Meier, 2012) The deep layers of the Baltic Sea are disconnected from the ventilated upper ocean layers, and temperature variations are predominantly driven by mixing processes and horizontal advection. The large-scale mean horizontal circulation is dominantly cyclonic in the Baltic (Meier, 2007;Jędrasik and Kowalewski, 2019). Deep water circulation consists of dense bottom currents of the inflowing saline water in the southern Baltic Sea, while convection, mixing, entrainment and vertical advection of water masses leads to interactions between upper and lower layers (Leppäranta and Myrberg, 2008). Instantaneous surface circulation pattern is driven by wind in the Baltic Sea.
The Baltic Sea has been suffering from eutrophication for at least half a century already. Out of external sources, rivers contribute most to the Baltic Sea nutrient input (HELCOM, 2018a). Seven biggest rivers -Daugava, Gota, Nemunas, Neva, Oder, Tornio and Vistula -cover 50% of the Baltic drainage area of 1.74 km 2 (HELCOM, 2018a). The average flow of the largest rivers Neva and Vistula is 2310 m 3 /s and 1112 m 3 /s, respectively, whereas smallest contributing rivers have a flow b600 m 3 /s (HELCOM, 2018a). Rivers with catchments in densely populated agricultural areas, such as Vistula, Nemunas and Oder in the southern part of the Baltic Sea have a bigger nutrient input compared to northern areas, where large parts of river drainage areas are under forest (HELCOM, 2018a). Although, nutrient input from the rivers has decreased by 12% in case of nitrogen and 25% of phosphorus over couple of decades, no improvement of the Baltic Sea state has been detected (HELCOM, 2018b). A warming trend of 0.08°C/year in the upper 50 m layer and 0.04°C/year in the deep layers (N60 m) reinforce a cascade of biogeochemical processes which is called the "vicious circle" (Savchuk, 2018;Meier et al., 2012;Vahtera et al., 2007). Due to the "vicious circle" there is enhancement of cyanobacteria blooms, which hampers nitrogen reduction attempts and sustains an elevated trophic state via accelerated oxygen consumption during organic matter oxidation, an increase of denitrification and nitrification rate in sediments and enhanced release of phosphates from the accumulated sediments due to hypoxia and anoxia (Seitzinger, 1988;Vahtera et al., 2007;Savchuk, 2018). On the other hand, the cyanobacteria blooms stimulate summer production in the entire food web, from zooplankton and benthos to fish (Karlson et al., 2015).

Methods
The effects of additional shipborne nutrients on the marine primary production is estimated using the coupled physical and biochemical model system GETM-ERGOM (Burchard and Bolding, 2002;Bruggeman and Bolding, 2014) for the Baltic Sea. Combining the simulation results from Automatic Identification System (AIS) based emission modelling using The Ship Traffic Emission Assessment (STEAM) and atmospheric deposition fields from the Community Multiscale Air Quality (COSMO-CLM/CMAQ) model system (Rockel et al., 2008;Matthias, 2008;Byun and Schere, 2006;Karl et al., 2018) using consistent STEAM shipping emissions to the atmosphere, the input of nitrate, ammonium, phosphate and organic matter are applied as mass fluxes to the surface layer of the sea. The year 2012 is considered as the reference year (SHIP model simulation). A NOSHIP model simulation is performed excluding the above-mentioned external input of nutrients from shipping activity. In general, annual mean temperatures were 0.5°C to 0.7°C above normal and it was wetter in the Nordic and Baltic regions in 2012 (Achberger et al., 2013). The year 2012 can be considered as typical for hydrographic and biogeochemical conditions relative to the climatological period 1993-2014. Horizontally averaged annual temperature and salinity profiles, as well as sea ice extent and volume were close to the mean (Von Schuckmann, 2019). The year 2012 shows decrease of salinity below the halocline in the Gotland Deep (Von Schuckmann, 2019) and relatively high spatial extent of oxygen depleted water accountic for the hypoxic area of 60,000 km 2 in the Baltic Sea (Savchuk, 2018). Previous saline water inflow which signal of the salinity increase was detected in the Gotland Deep took place at the end of 2006. The mean total freshwater discharge into the Baltic Sea for the year 2012 was 14% higher than long-term average in year 2012 (Kronsell and Andersson, 2013).
Spring bloom started relatively early, in March already, with the start date varying little across the entire Baltic Sea area (Raudsepp et al., 2019a). Peak bloom day stretched from the end of March to the end of April. The spring bloom ended at the end of May, similarly to the other years since 2007. The spring bloom's spatiotemporal coverage in 2012 was close to the mean, but the phytoplankton summer bloom was among the smallest (Raudsepp et al., 2018). The latter is consistent with a minimum of the interannual oscillations of cyanobacterial blooms with a period of about 3 years in the Baltic Sea (Kahru et al., 2018).
TN and TP pools were in a stable level of about 6000 ktons and 680 ktons, respectively since 2005, as well as DIN and DIP pools of about 1000 ktons and 480 ktons, respectively (Savchuk, 2018). In general this is consistent with the period of stable nutrient input, within otherwise decreasing trend since 1980 (Savchuk, 2018).
The modelling system consists of the ship emission model, The Ship Traffic Emission Assessment Model (STEAM) (Jalkanen et al., 2009;Jalkanen et al., 2012;Johansson et al., 2013;Johansson et al., 2017), the atmospheric chemistry modelling system, Climate Limited-area Modelling Community (COSMO-CLM) and The Community Multiscale Air Quality (CMAQ) model system (Rockel et al., 2008;Matthias, 2008;Byun and Schere, 2006;Karl et al., 2018), and the coupled marine physical model, General Estuarine Transport Model (GETM) (Burchard and Bolding, 2002;Bruggeman and Bolding, 2014), and the biogeochemical model, the Ecological Regional Ocean Model (ERGOM) (www.ergom. net; Neumann and Schernewski, 2008). The STEAM is ship Automatic Identifiation System (AIS)-based emission model providing shipping emissions for the CMAQ model system and direct ship discharges to the water. The atmospheric deposition fields from the CMAQ model provide the input of nitrate, ammonium, phosphate and organic matter mass fluxes to the surface layer of the sea. Further on, the GETM-ERGOM system uses the atmospheric deposition fields and direct ship discharges for the estimation of the effects of shipborne nutrients on the Baltic Sea ecosystem.
The model system simulations were carried out for the year 2012. In general, annual mean temperatures were 0.5°C to 0.7°C above normal and it was wetter in the Nordic and Baltic regions in 2012 (Achberger et al., 2013). The year 2012 can be considered as typical for hydrographic and biogeochemical conditions relative to the climatological period 1993-2014. Horizontally averaged annual temperature and salinity profiles, as well as sea ice extent and volume were close to the mean (Von Schuckmann, 2019). The year 2012 shows decrease of salinity below the halocline in the Gotland Deep (Von Schuckmann, 2019) and relatively high spatial extent of oxygen depleted water accountic for the hypoxic area of 60,000 km 2 in the Baltic Sea (Savchuk, 2018). Previous saline water inflow which signal of the salinity increase was detected in the Gotland Deep took place at the end of 2006. Freshwater discharge. Spring bloom started relatively early, in March already, with the start date varying little across the entire Baltic Sea area (OSR 3). Peak bloom day stretched from the end of March to the end of April. The spring bloom ended at the end of May, similarly to the other years since 2007. The spring bloom's spatiotemporal coverage in 2012 was close to the mean, but the phytoplankton summer bloom was among the smallest (Raudsepp et al., 2018). The latter is consistent with a minimum of the interannual oscillations of cyanobacterial blooms with a period of about 3 years in the Baltic Sea (Kahru et al., 2018). TN and TP pools were in a stable level of about 6000 ktons and 680 ktons, respectively since 2005, as well as DIN and DIP pools of about 1000 ktons and 480 ktons, respectively (Savchuk, 2018). In general this is consistent with the period of stable nutrient input, within otherwise decreasing trend since 1980 (Savchuk, 2018).

STEAM model description
The STEAM (Jalkanen et al., 2009(Jalkanen et al., , 2012Johansson et al., 2013Johansson et al., , 2017 uses Automatic Identification System (AIS) data to describe ship traffic activity and the detailed technical knowledge of the ships for the calculation of the atmospheric emissions and discharges directly from the ships to the sea. In this work, new capabilities were built into STEAM, which enabled the description of various discharges, like Black Water (BW), Grey Water (GW), Ballast Water (BLW), Bilge Water (BLG), Food Waste (FW), Scrubber water from open and closed loop systems (SWO, SWC), Stern Tube Oil (STO) and several species of antifouling paints (AFP) from ships. These developments of STEAM are described in a separate manuscript (Jalkanen et al., in preparation). BW is defined as sewage, wastewater that originates from toilets, medical facilities, premises for living animals or other wastewaters when mixed with those drainages (MARPOL Annex IV). GW is collected from dishwater, shower, laundry, bath and wash-basin drains, and its discharges are not limited by the international law (MARPOL Annex V). BLG is a mixture of different substances from machinery, spills and overflow tanks and gets accumulated in the lowest part of the ship (Klein Wolterink et al., 2004;IMO, 2006). FW generated on-board can be any 'spoiled or unspoiled' foods and food scraps (MARPOL Annex V). GW, BW and FW are functions of people carried onboard, both crew passengers.
Vessel activity for year 2012 was described with AIS data sent by the Baltic Sea fleet and provided to us by the Helsinki Commission (HELCOM), which consists of all Baltic Sea countries. This data consists of over 320 million automatic position reports. In addition, technical description of the fleet was based on data from IHS (IHS Global, 2016). Passenger capacity utilization was modelled based on the quarterly reports of major passenger vessel operators, who together carry over 20 million passengers each year. Based on these reports, passenger capacity utilization was estimated as 50% throughout the year, except for cruise vessels, which have been reported to use close to 90% of their capacity (HELCOM, 2015;Wilewska-Bien et al., 2018). The time spent on-board was separately estimated for the crew and passengers based on AIS data. This resulted to estimates of the discharge amounts. The spatiotemporal releases of accumulated quantities were modelled as defined by the MARPOL convention Annexes I, IV and V, which define what, where and how discharges can be released to the sea. Some of the releases occur randomly, like GW, BW and BLG, which were modelled as continuous discharges in areas defined by current MARPOL rules and national environmental legislation if they are stricter than MARPOL (like oily releases in Finnish waters and SWO discharges in German waters).
The atmospheric emissions were delivered as input for the chemical transport model CMAQ, which was used to calculate the atmospheric transformation and transport of pollutants. The STEAM outputs consisting of gridded daily inventories of nutrients in BW, GW and BLG as well as total volumes of discharges were used as input data for the GETM-ERGOM.

CMAQ model description
The CMAQ model used for calculation of the deposition fields (Byun and Schere, 2006;Appel et al., 2013Appel et al., , 2017 is an atmospheric Eulerian chemistry transport model. It computes atmospheric concentrations and deposition of numerous trace gas species and aerosol components, depending on emissions and the physical state of the atmosphere. CMAQ was developed by the US EPA about 20 years ago (Byun and Ching, 1999). It is freely available through the CMAS center and permanently updated. In the SHEBA project, the model was run in version 5.0.1.
The model was set up in three nested domains, a grid with 64 × 64 km 2 for entire Europe, 16 × 16 km 2 for northwestern Europe and 4 × 4 km 2 over the Baltic Sea. In the vertical, the model extends up to 100 hPa in a sigma hybrid pressure coordinate system with 30 layers. Twenty of these layers are below approximately 2 km; the lowest layer extends to ca. 36 m above ground. The entire year 2012 was run with a spin-up period of one month for the initialization of the model runs. Meteorological fields were calculated with the COSMO-CLM model (Rockel et al., 2008) and interpolated to the CMAQ grid with an adapted version of the Meteorology-Chemistry Interface Processor (MCIP) (Otte and Pleim, 2010).

GETM description
The GETM is a numerical hydrodynamics model which is solving sea state by means of salinity, temperature, currents and water level (Burchard and Bolding, 2002). The modular concept of GETM makes it possible to use various parameterizations and numerical techniques to solve numerically three-dimensional primitive ocean equations. The time split technique allows using shorter timestep for solving free-surface evolution in barotropic mode and longer timestep for internal baroclinic mode solving transports. Spatially the equations are discretized on a staggered Arakawa-C grid using spherical horizontal coordinates and bottom-following adaptive layers discretization in vertical direction. The vertical resolution has been enhanced by reducing layer thicknesses near the boundaries and in the vertical ranges where stratification is strong. Such adaptive coordinate system produces less numerical mixing compare to general sigma-coordinate discretization (Gräwe et al., 2015). The horizontal advection has been solved using third-order TVD scheme with P2-PDM limiter. Directional splitting technique has been applied according to Pietrzak (1998). To minimize known pressure gradient errors internal pressure has been solved using zinterpolation method according to (Shchepetkin and McWilliams, 2003). The vertical subgrid turbulence is solved using k-ε model using algebraic turbulence closure for momentum equations has been applied via General Ocean Turbulence Model (GOTM, Umlauf and Burchard, 2005). Background vertical diffusivity has been set to 10 −6 m/s. Horizontal viscosity coefficient has been set to 10 m 2 s −1 according to Wallcraft et al. (2005) considerations. Air-ocean momentum and heat fluxes were calculated using Kondo (1975) bulk parametrizations.
The model domain covers the whole Baltic Sea with closed boundary in the Kattegat. The bathymetry has been derived from the Baltic Sea Digital Database (BSBD 0.9.3) and interpolated on a 1 nautical mile grid, which is the horizontal resolution of the model. In vertical direction 40 bottom-following and adaptive layers have been defined, ensuring b5 m vertical resolution in halocline and near the surface. The timestep for the biogeochemical processes is set to baroclinic timestep, which was 500 s.

ERGOM description
ERGOM (Neumann, 2000(Neumann, , 2002 is an extended version of the Nbased NPZD model taking into account processes like N-fixation, phosphorus limitation, denitrification and phosphorus binding into iron compounds. 12 state variables are used which describe the N cycle in molar N units. The inorganic nutrients, which are consumed by primary producers, are defined as dissolved nitrate, ammonium and phosphate. Primary producers are divided into three functional phytoplankton groups: diatoms, flagellates and N-fixing cyanobacteria. Chlorophyll-a (chl-a): is the sum of all three functional groups of algae and detritus. Nitrogen in phytoplankton and detritus (mol N/kg) is converted into molar carbon-content according to the Redfield ratio (Redfield, 1934). Grazing of phytoplankton is described as the growth of zooplankton. Phyto-and zooplankton are transformed into dead organic matter, which sinks and contributes to the sediment pool as detritus. Under oxic conditions part of the detritus is remineralized back to dissolved nutrients; this process uses oxygen and has a temperature-dependent rate. Under anoxic conditions denitrification reduces nitrate to molecular nitrogen, which leaves the system. If nitrate is depleted under anoxic conditions, detritus is oxidized with sulfate and generates dihydrogen sulfur gas, which is considered as negative oxygen. Under oxic conditions reactive phosphates are bound into iron-phosphates, which sink out of the water-column and accumulate in the sediment layer. In case of anoxia with the presence of sulfuric acid, iron-oxide gets reduced and phosphates are released back to the system as nutrients available to the primary producers. A fraction of these iron-phosphate complexes is also assumed to be buried permanently, depending on the sediment thickness and sedimentation rate (Neumann and Schernewski, 2008). More detailed descriptions of the ERGOM model are available in Neumann (2000); Neumann et al. (2002); Neumann and Schernewski (2008); Radtke et al. (2012) and Lessin et al. (2014a).
The initial conditions for salinity and temperature have been taken from a hindcast simulation of Maljutenko and Raudsepp (2014) to avoid a long spin-up period and instabilities, which occur during cold starts. Stable density fields from December 2000 were chosen based on similar measured salinity values and stratification conditions in the Gotland Deep in 2012. A short spin-up period of one month was applied to adjust the model with new atmospheric forcing and river inflows. Initial nutrient pools have been adopted from climatic average fields for January from the 40 year long hindcast ERGOM simulation for the Baltic Sea (Kõuts et al., 2019). The oxygen profile measured in 2011 December at the Gotland Deep has been assimilated to the initial field of oxygen for the whole Baltic Proper domain to adjust the extent of the anoxic area for the year 2012. The open boundary conditions have been closed to neglect the nutrient fluxes from the North Sea and therefore to study the Baltic Sea as a closed system. This assumption affects the biogeochemistry in the Kattegat and southwestern Baltic but as we deal with one year and one year simulation repeated for five years, we assume the impact to be minor in the context of the entire Baltic Sea.

Coupled GETM-ERGOM justification
Placke et al. (2018) have compared main hydrographic and circulation fields of the long-term model simulation results performed by GETM, MOM and RCO circulation models in the Baltic Sea. Their assessment remains without general conclusion that all considered fields are reproduced better by one of the model than by the others. Therefore, we do not have solid justification of using one particular circulation model in our study. Although, Gräwe et al. (2015) have argued that the adaptive vertical coordinate system (Burchard and Beckers, 2004;Hofmeister et al., 2010) implemented in GETM has an advantage in comparison to fixed vertical coordinate system implemented in MOM and RCO by having less numerical mixing due to a better resolution of a strong vertical stratification. Our argument of using GETM is that we have successfully validated the model for long-term simulation period of 1996-2006(Maljutenko and Raudsepp, 2014. Comprehensive comparison of the Ecological Regional Ocean Model (ERGOM) (www.ergom.net; Neumann and Schernewski, 2008) coupled to the physical Modular Ocean Model (MOM 3.1) (e.g. Pacanowski and Griffies, 2000), the Swedish Coastal and Ocean Biogeochemical (SCOBI) model (Eilola et al., 2009;Almroth-Rosell et al., 2011) coupled to the Rossby Centre Ocean (RCO) circulation model (RCO-SCOBI) and the BAltic sea Long-Term large-Scale Eutrophication Model (BALTSEM) (Gustafsson, 2003) long-term simulation results of nutrient and oxygen dynamics in the Baltic Sea show no advantage of any model . The ERGOM and SCOBI models are relatively similar in terms of the set of state variables used and key biogeochemical processes implemented in the models (e.g. Eilola et al., 2011).
The SCOBI model has been mainly used in the coupling with RCO model for climate change studies (e.g. Meier et al., 2018) and recently coupled with NEMO (Raudsepp et al., 2019b). The MOM-ERGOM coupled model system has been widely used for the biogeochemical studies of the Baltic Sea (e.g. Neumann et al., 2017). A coupled GETM-ERGOM model has been used for different studies of the biogeochemistry of the Baltic Sea in multi-year simulations (Lessin et al., 2014a(Lessin et al., , 2014b and in particularly for investigating the ship impact to marine biogeochemistry (Raudsepp et al., 2013). Besides, the model system of the 3D ocean circulation Hiromb-BOOS model (HBM) coupled with the biogeochemical ERGOM model (Maar et al., 2011(Maar et al., , 2014Wan et al., 2012) has been used for studying nutrient loads impact to primary production in the western Baltic Sea (Maar et al., 2016) and for evaluation of atmospheric nitrogen inputs including ship-borne nitrogen to marine ecosystems (Neumann et al., 2018a(Neumann et al., , 2018b(Neumann et al., , 2018c.

Meteorological forcing
Meteorological fields to drive the CMAQ and GETM-ERGOM system were calculated with COSMO CLM version 5.0 (Geyer, 2014). COSMO was driven with ERA-Interim reanalysis fields using the spectral nudging technique to force the model to stay close to the reanalysis. The runs were performed on a 0.11°× 0.11°rotated lat/lon grid with 40 vertical layers up to approx. 20 km altitude. A nested model run on a 0.025°× 0.025°grid with 50 vertical layers was applied to the Baltic Sea region and then interpolated on a 4 × 4 km 2 grid with a modified version of the MCIP (Otte and Pleim, 2010). We have applied hourly atmospheric forcing for the GETM-ERGOM model system.

Input of nutrients
Nutrient input from the 36 largest rivers has been derived from the European hydrology model (E-HYPE) hindcast simulation of Donnelly et al. (2015). The diffuse inputs of adjacent regions along the coast have been added to the nearest rivers. Since only total nitrogen (TN) was available from hindcast simulations, we used fractions of 0.70 and 0.05 for nitrate and ammonium, respectively. The remaining fraction of 0.25 from TN contributes to the organic nitrogen pool as detritus input. The input of phosphates and organic phosphorus was assumed 0.25 and 0.75 fraction of the total phosphorus, respectively. The used fractions were inferred from numerous studies on nutrient loads to the Baltic Sea (e.g. Stålnacke et al., 1999;Vahtera et al., 2007;Savchuk et al., 2012).
Nutrient input from shipping to the Baltic Sea has been considered from two main sources: atmospheric depositions and discharge from shipping. The emission of NOx to the atmosphere from shipping has been taken from hourly emission dataset of STEAM simulation covering the full calendar year of 2012. In STEAM, each vessel is considered as unique case considering vessel specific ship activity and technical description. Daily emissions of NOx, SOx, CO, Elementary Carbon (EC), Organic Carbon (OC), SO 4 and Ash were generated as emission grids, which were used as input for atmospheric modelling.
Emissions for all other anthropogenic sources except shipping are based on EMEP emission data for the year 2012. They include NOx, NH 3 , SO 2 , CO, NMVOCS and PM10. Biogenic emissions were calculated with BEIS 3.4 (Schwede et al., 2005;Vukovich and Pierce, 2002). Sea salt emissions followed the parametrization of Gong (2003), but excluding surf zone emissions because of overestimation of the emission flux in some regions (Neumann et al., 2016). All emissions were temporally and spatially disaggregated and distributed to the CMAQ model's grids with the SMOKE for Europe emission model (Bieser et al., 2011a). The vertical distribution of the emissions from certain sectors followed the standard profiles given in Bieser et al. (2011b).
Atmospheric deposition of NOx species into the sea have been considered from the following wet and dry depositions of: particulate NO 3 (PM_NO 3 ), nitrogen trioxide (NO 3 ), nitrogen dioxide (NO 2 ), nitric oxide (NO), nitric acid (HNO 3 ), nitrous acid (HONO), nitrogen pentoxide (N 2 O 5 ), peroxyacetyl nitrate (PAN), oxidized peroxyacetyl nitrate (OPAN) and peroxynitric acid (HNO 4 ). Particulate ammonium (PM_NH 4 ) and ammonia (NH 3 ) contribute to reduced nitrogen deposition. The atmospheric phosphorus deposition has been considered as 3.5% of mineral ash deposition calculated in CMAQ. The hourly atmospheric depositions from a 4 × 4 km 2 atmospheric deposition grid have been interpolated on the 2 × 2 km 2 shipping discharge grid used for the direct nutrient inputs.
Discharge of nutrients from BW, GW, FW, SWO and SWC have been calculated from the discharge volumes and concentrations in wastewaters (Wilewska-Bien et al., 2018;Jalkanen et al., in preparation., Wilewska-Bien et al., 2016, Jönsson et al., 2005, McLaughlin et al., 2014, Furstenberg et al., 2009, Hufnagl et al., 2005. The organic nitrogen is considered as dead organic matter input and is recycled to dissolved nutrients in the water by bacteria. The actual numbers of different inputs are given in Table 1.

GETM-ERGOM validation
GETM-ERGOM system validation is performed for the year 2012. We have calculate model and data bias and root mean square error (RMSE) for basic state variables as water temperature, salinity, nitrate, phosphate, oxygen and chlorophyll a (Chl-a) at the station BY15 of the Gotland Basin. Surface layer temperature is significant abiotic variable that controls the seasonal cycle of the biogeochemical processes and is well reproduce by the model (Fig. 2a, Bias = −0.5°C, RMSE = 0.8°C). Bottom temperature (Fig. 2a, Bias = −0.3°C, RMSE = 0.4°C) and bottom salinity (Fig. 2b, Bias = 0.7 g/kg, RMSE = 0.7 g/kg) do not vary significantly in time, as their dynamics is mainly related to the Major Baltic Inflows which are absent in 2012 (Raudsepp et al., 2018). Surface salinity (Fig. 2b, Bias = 0.5 g/kg, RMSE = 0.5 g/kg) as well as entire vertical temperature (Fig. 2d) and salinity (Fig. 2e) profiles are well reproduced except that halocline is about 10-m deeper in the model than in the measurements. The latter affects vertical distribution of nitrate (Fig. 2f), phosphate (Fig. 2g) and oxygen (Fig. 2h). Seasonal variations of nitrate (Bias = 1.1 mmol N/m 3 , RMSE = 1.5 mmol N/m 3 ), phosphate (Bias = −0.1 mmol P/m 3 , RMSE = 0.2 mmol P/m 3 ) and chl-a (Bias = −0.9 mg/m 3 , RMSE = 1.0 mg/m 3 ) in the surface layer are depicted in Fig. 2c. Their basic seasonal cycle is well reproduced by the model -nutrients concentrations are high in spring before the spring bloom of diatoms, depleted in summer and start to increase in autumn; spring diatom bloom and summer cyanobacteria bloom are present in the model results. Deficiencies of the model are that phosphate concentration is underestimated in spring and spring bloom is delayed in time (Fig. 2c). Large model errors in the nitrate concentration between the depth of 60 and 150 m (Fig. 2f) and in the phosphate concentration from 60 m depth to the bottom (Fig. 2g) originate from the errors in the model initial fields (not shown) indicated by the variability of the model profiles of corresponding variables.
The GETM-ERGOM system consists a number of nonlinear related equations and calibration parameters. To perform mathematically correct uncertainty estimations and model sensitivity analyses is complex task on its own. Besides of the uncertainties of the GETM-ERGOM system, we have uncertainties related to the STEAM model and their treatment of emission and discharge of each particular ship, uncertainties related to COSMO CLM and CMAQ model; uncertainties due to river load of nutrients and uncertainties related to the GETM-ERGOM setup. The provided list of potential sources of uncertainties is far from being complete. In the present application of the GETM-ERGOM we rely on the wide use of this model system for similar research (Schernewski and Neumann, 2005;Neumann and Schernewski, 2008;Radtke et al., 2012;Lessin et al., 2014a;Lessin et al., 2014b;Neumann et al., 2018aNeumann et al., , 2018bNeumann et al., , 2018cKõuts et al., 2019). Therefore, we have kept the set of calibration parameters similar to the previous studies and rely on the validation results of the model. The strong argument for the reliability of the model results and conclusion drawn from the study is that we have performed two identical model simulation, i.e. with and without input of the ship-borne nutrients. Mainly, differences of the spatiotemporal distribution of the biogeochemical variables and their fluxes are analysed, not the absolute values. The calculated RMSE values are considered as proxies for the uncertainties of state variables.
The choice of biogeochemical model parameters and systems dynamics is sensitive to the physical model parameters and external forcing (Burchard et al., 2006;Miladinova and Stips, 2010).

Sources of nutrients from shipping
Annual input of nitrate is two orders of magnitude larger than the input of other nutrient compounds (Fig. 3a). Nitrate deposition consists of atmospheric deposition and direct discharges to the water. Atmospheric deposition of nitrate prevails over the direct discharge, which results in continuous distribution over the Baltic Sea. There is a southwest gradient of the nitrate input distribution due to increase of shipping intensity towards the southern Baltic Sea and Danish Sounds (Sime, 2014) and due to prevailing atmospheric circulation. Prevalent winds are from W/SW, so the area receives shipping emissions from the North Sea as well.
N-NH 4 input also combines atmospheric deposition and ship discharge (Fig. 3b). High input from discharges is observed on the shipping lanes, while more continuous spatial map of the atmospheric input is negative. The negative values increase towards the southwestern Baltic according to the increase of shipping activity. Discharge of N-NH 4 from ships, which originates mainly from BW, are compensated by reduced atmospheric deposition of N-NH 4 in the SHIP case compared to the NOSHIP case along the main shipping routes. Input remains positive on the routes of passenger ferries in the Gulf of Finland and northern Baltic Proper. Same holds for the southwestern Baltic and the Danish Sounds, although reduced deposition of ammonium in the SHIP case compared to the NOSHIP case is high there. P-PO 4 , which enters the sea via direct discharge, has a distribution pattern with elevated concentrations directly on the shipping lanes, where the input takes place (Fig. 3c).
Discharge of organic compounds of nitrogen and phosphorus takes place in the shipping lanes (Fig. 3c,d). BW, GW and FW are discharged to the sea in shipping lanes, but only allowed outside of 12 Nm from the coast. Only not comminuted or disinfected discharges were taken into account because at the time of modelling regulations did not require nutrient reduction in wastewater. The untreated wastewater was modelled beyond 12 Nm from coastline and when ships were en route moving at 4 knots. The wastewater generated within the 12 Nm boundary is collected onboard the ships and released outside the 12 Nm zone with a determined rate. Parameters like kinetics, mixing or transport (advection) were not included in the modelling. In the case of cargo ships, it results in point-discharge at the start of the shipping lanes outside of the 12 Nm coastal zone. The much larger volumes released from passenger ferries are assumed to be continuous, which results in a line discharge. The discharge rates for different type of ships are included in the modelling runs.

Time series of marine environment response to ship deposition compared to natural processes
The Gotland Basin was selected to study the immediate effect of shipborne nutrient input on the biogeochemistry of the Baltic Sea. Time series of main biogeochemical parameters were selected from the BY15 monitoring station (Fig. 1), which is a HELCOM monitoring Table 1 Annual total excess nutrient input (SHIP minus NOSHIP), nutrient input in the form of direct discharge from the ships and percentage of contribution from different sources to direct discharge.

Annual nutrient input
Share of shipping emissions  (Fig. 1). BY15 as a deep station (250 m) represents deep-layer biogeochemical processes, which are detrimental in steering the eutrophication process in the Baltic Sea. A major part of the research on the influence of MBIs on hydrophysical and biogeochemical conditions in the Baltic Sea are based on the measurements there (e.g. Savchuk, 2018;Mohrholz, 2018). It is also considered as a representative area of spring and summer blooms -with medium to high chl-a concentrations (Kahru et al., 2007;Zhang et al., 2018).
Nitrate, phosphate, phytoplankton and zooplankton variations undergo the commonly known annual cycle in the Baltic Sea without visible differences between the SHIP and NOSHIP simulations. In order to estimate the impact of the SHIP scenario on the Baltic Sea ecosystem we calculated the relative changes of the main biogeochemical variables as (SHIP-NOSHIP)/SHIP.
Nitrate and phosphate concentration increases on the surface until April-May and go into steep decline mid-May (Fig. 4a,b). This coincides with the diatom bloom, which peaks at the end of June when nitrate becomes depleted (Fig. 4c). Simultaneously to the steep decline of nitrate  and phosphate at the end of May, shipping-derived relative changes in nitrate and phosphate become apparent: more nitrate (max. +10%) and less phosphate (max. −3%) is available on the surface, in SHIP case compared to NOSHIP case (Fig. 4a,b). This lasts from June till the end of July/beginning of August. Relative changes of nitrate and phosphate are then reversed, starting at the end of August, with less than average surface nitrate being available from August to the second half of September and more phosphate from the beginning of August to the end of September. Relative increase (max. +5%) of diatoms is evident mostly after their bloom, whereas flagellates experience a relative increase (+12-13%) before their bloom peak (Fig. 4c,d). Cyanobacteria go through a relative decrease (−10%) from June to October and coincides with the relative shipping-related increase of diatoms and flagellates (Fig. 4e). Zooplankton experiences very little relative change due to shipping, with a slight increase (+1%) during their bloom in end of May/June and decrease (−4%) from August to October (Fig. 4f), slightly shifted from the cyanobacteria bloom peak (Fig. 4e).
Snapshots of vertical profile of the biogeochemical variables at station BY15 show that the surface layer is representative of the entire upper mixed layer and shipping-related impact does not go below the pycnocline, except for oxygen (Fig. 5). While nutrient concentrations are determined by the depth of the mixed layer (around 75 m) (Fig. 5a,b,c), phytoplankton and zooplankton depend on the euphotic zone in the upper 30 m (Fig. e,f,g,h). Changes in dissolved oxygen are most pronounced at the transition depth from hypoxic to anoxic conditions (concentration of dissolved oxygen between 1 and 2 ml l −1 ), which is at 150 m depth at the BY15 station (Fig. 5d). Relative changes are most severe where oxygen values are close to zero initially. Absolute decrease in oxygen results in an increase of the anoxic bottom area by 50 km 2 at the end of the year compared to the total anoxic area of 40,000 km 2 in NOSHIP scenario.

Spatial distribution to complement time series
The exact location of the shipping lanes does not have a substantial effect on the spatial distribution of the surface excess nitrate (Figs. 1, 6a) as NOx is mainly emitted to the air, where it undergoes atmospheric chemistry transformation and redistribution due to atmospheric circulation. This means that the spatial pattern of atmospheric deposition of nitrate is strongly dispersed compared to the concentrated discharge pattern of e.g. phosphate (Fig. 3c).
Deposited atmospheric nitrate is redistributed in the marine environment due to currents and mixing, so there is no resemblance of excess nitrate (Fig. 6a) to the deposition pattern (Fig. 3a). The Bothnian Bay has the lowest excess nitrate in the surface layer, which corresponds to the low nitrate input in the area. The rest of the Baltic coastal areas have higher surface concentrations of excess nitrate than the open sea areas. This pattern of redistribution is explained by shallower water column near the coast where the water depth is smaller than the mixing depth. For instance, in the case of equal deposition of nitrate to the deep and shallow area relative to mixing depth, the nitrate is equally mixed over the surface mixed layer in the deep area, but equally mixed over the entire water column in the shallow area. As the layer thickness is bigger in the first case than in the second case, the concentration of nitrate is lower in the first case than in the second case. The excess nitrate has highest concentrations on the eastern coast of the Baltic.
The mean distribution of excess phosphate in the water has a fragmented pattern (Fig. 6b). The area of high excess phosphate is clearly seen in the northeastern Baltic Proper due to the continuous line discharge from passenger ships. The other region of elevated excess phosphate concentration is in the southern Baltic around the Danish Sounds, also explained by very heavy ship traffic which is concentrated in a relatively small sea area.
The elevated excess diatom concentrations correspond to the excess nitrate, most notably in the Estonian Archipelago Sea (Fig. 6a,c). Elevated concentrations are also evident in the northern Baltic Proper and the northern part of the Gulf of Finland, in the central part of the Gulf of Bothnia, on the Polish coast and in the Kattegat. No excess diatoms are seen in the phosphorus limited areas, e.g. southern and eastern Gulf of Riga, and the eastern coastal area of the Baltic Proper.
The distribution of excess flagellates, which bloom at the end of summer and in autumn, mostly concentrates into Kattegat and less intensively in the Northern Baltic Proper -Gulf of Finland area (Fig. 6d), which slightly reflects the pattern of excess phosphorus on the surface. The spatial extent of cyanobacteria blooms is reduced in vast areas in the Baltic Proper, the Gulf of Finland and several coastal areas (Fig. 6e). The areas with little to no change in cyanobacteria blooms overlap with the areas with phosphorus limitation -e.g. Bothnian Bay and big river estuaries. Deposition of organic matter into the sediment was calculated during two periods: intense phytoplankton bloom period May-June and decomposition period in October-December (Fig. 6f,g). Maximum excess sedimentation in the period of May-June is heaviest in the southern Baltic (south of Gotland), the West-Estonian archipelago, Gulf of Finland and the Aland Sea. The Gulf of Bothnia, the southern part of the Gulf of Riga and areas on the eastern coast of the Baltic Proper are less affected by excess organic matter. The distribution of surplus sediments is spread uniformly over the Baltic Proper and the Gulf of Finland by water circulation. Organic matter accumulation in May-July is spatially more extensive and contains more organic matter due to primary production than in October-December when sediments are already being remineralized and the flux from the surface is declining. In autumn organic matter mostly accumulates in the medium deep and deep areas in the central part of the basins. The period from October to December reflects a more long-term distribution of sediments, shaped by transport. Another distinctive feature of the sediment distributions from both periods (Fig. 6f,  g) is the small-scale spatial heterogeneity. This is related to the smallscale topographic irregularities that hinder smooth horizontal transport and favor sedimentation on the slopes of topographic shoals.
The pattern of organic matter accumulation in the two periods is in most parts similar to oxygen distribution, with decreased oxygen areas overlapping with high organic material accumulation spots (Fig. 6h,i). Summer bottom oxygen reduction is mostly seen in shallow, coastal areas in the Archipelago Sea in Estonia, as well as the coasts of southern Finland and Sweden. Large areas are affected by shipinduced oxygen decline west of Bornholm island and in the Danish Straits, where marine traffic is also intensive. Central deep and mostly anoxic areas of the Baltic Proper are not affected by ship-induced oxygen decline. Winter oxygen conditions show a different pattern with minimum concentrations in the areas of medium depth, i.e. 60 m depth. The shallower sea areas are well ventilated by wind and thermal mixing and are not affected. The most severe impact is seen in the Kattegat, southern Baltic, Bornholm Basin, Stolpe Channel and Gdansk Basin (Fig. 6h, i). Winter oxygen conditions reflect the accumulation impact of shipborne nutrients and represent oxygen conditions that remain poor after dead organic matter has decomposed. Longitudinally, the reduction of bottom oxygen is stronger along the eastern bottom slope than along the western slope of the Baltic Proper due to the prevailing cyclonic circulation (Omstedt et al., 2014).  (Fig. 3). The time instance for the oxygen profile d) corresponds to cyanobacteria bloom peak. Relative difference between the SHIP and the NOSHIP scenario is shown as green line. Potential density (σ t ) is shown as yellow line.

Discussion
This study is based on one year of ship traffic data, the corresponding emissions and discharges and the atmospheric and marine conditions of the year 2012. The intensity of ship traffic can vary annually (HELCOM, 2010) and the same holds for atmospheric and marine conditions, which influence the actual deposition of nutrients to the sea and their spreading. A study by Raudsepp et al. (2013) for repeated ship deposition in changing hydrodynamic conditions in the Gulf of Finland shows interannual variations of different biogeochemical variables in the range of 50%. The ship routes do not change from year to year, hence the general pattern of direct nutrient discharge from ships remains unchanged, which means that main discharge areas are southwestern and northeastern Baltic Proper. Nitrogen discharge is two orders of magnitude lower than atmospheric deposition of nitrogen. Atmospheric nitrogen deposition could be inter-annually more variable in terms of the spatial distribution. However, as the results show a rather smooth and only slightly varying spatial nitrogen deposition pattern, Fig. 6. Spatial distributions of variable differences between the SHIP and the NOSHIP scenario. The difference in surface concentrations at the time of their maxima in the NOSHIP scenario during the period from February to June, for a) nitrate, and b) phosphorus. The difference of c) diatoms, and d) flagellates, during the time of the peak bloom in spring for the NOSHIP scenario. e) Difference of summer bloom (cyanobacteria) during the peak bloom for the NOSHIP scenario. Distribution of the maximum difference of the nitrogen in sediments f) after the spring bloom from May to July, and g) after the bioactive period from October to December. The distribution of maximum differences between the SHIP and the NOSHIP scenario for the near-bottom oxygen concentrations during the summer months from h) July to September, and i) for the period after the bioactive season from October to December.
we expect that atmospheric circulation will redistribute NOx emitted to air rather homogeneously in the other years, as well. The spatial distribution of excess nitrate, phosphate and organic matter from shipping discharge is more heterogeneous than atmospheric deposition of nitrate. The circulation can vary a lot in the Baltic Sea from year to year (Lehmann et al., 2002;Meier and Kauker, 2003;Väli et al., 2013), so the areas of high and low nitrate concentration can vary accordingly between years.
Diatoms and flagellates benefit from the ship-borne nutrients which reflects in the increase of their biomass. The spatial distribution of excess diatoms and flagellates might depend on the competition for the nitrogen between two groups. Nutrients play an important role in the structure of phytoplankton communities, which can be explained by varying nutrient-acquiring abilities of different species (Lomas and Glibert, 2000). However, it must be taken into consideration that it is difficult to determine exactly how important nutrients are, as there are also many other factors at play (which we do not look at in this study) regarding the competition between diatoms and dinoflagellates -e.g. other nutrients (silica), light conditions, water stratification, ice conditions and specific characteristics of the species, e.g. mixotrophy, which gives a competitive advantage in some nutrient enrichment situations (Lagus et al., 2004;Granéli et al., 1990;Lomas and Glibert, 2000;Spilling et al., 2014).
Cyanobacteria biomass is reduced the most in the areas with usually strong cyanobacteria blooms, which include southern and central Baltic and the western part of the Gulf of Finland (Kahru et al., 1994(Kahru et al., , 2007Kahru and Elmgren, 2014). This result is in accordance with previous studies on the effect of additional nitrogen on cyanobacteria blooms (Elmgren and Larsson, 2001). The areas with little to no change in cyanobacteria blooms overlap with the areas with phosphorus limitation -e.g. Bothnian Bay and big river estuaries (Kõuts et al., 2019). Excess detritus distribution and corresponding oxygen reduction reflect the pattern of excess diatom distribution and possibly water circulation in summer. The same is generally true for autumn except that shallow areas are well ventilated then. Deep areas in the central and northern Baltic Proper are unaffected as they are anoxic already (Fig. 6h).
The biogeochemical cycle in the Baltic Sea is defined by nitrogen being mostly the limiting nutrient. This is due to the system's ability to remove the added biologically available nitrogen through the microbial conversion back to N 2 gas via denitrification, which each year removes nitrogen corresponding to most of the input (Voss et al., 2005;Radtke et al., 2012). This takes place in the sediments as well as in oxygen-deficient deep waters and becomes more effective when oxygen-deficiency is widespread (Voss et al., 2005;Radtke et al., 2012). Therefore, any changes in nitrogen input affect the Baltic Sea ecosystem functioning.
To assess the impact of shipping-related excess nitrogen on the biogeochemical processes in the Baltic Seat we subtract fluxes of nitrogen from two model simulations, i.e. SHIP minus NOSHIP, like we did for the biogeochemical variables. The nitrogen balance in the marine system takes into account nitrogen content in the sea, both in water and sediments, nitrate and ammonium flux from the atmosphere and rivers, nitrogen brought to the system by nitrogen fixation and outflow of nitrogen to the air in the form of molecular nitrogen, i.e. the result of denitrification. River fluxes of nitrogen are the same in both cases giving zero contribution to excess nitrogen.
The results show that since the beginning of January until the end of March, shipborne nitrogen input exists in the form of nitrate in the water column (Fig. 7a). When the spring bloom starts, not all nitrogen load into the water remains there, i.e. total nitrogen in the ecosystems becomes lower than the shipborne nitrogen input. Total nitrogen in the marine system consists of nitrate in the water as the main contributor, nitrogen bound in phytoplankton, dissolved molecular nitrogen and nitrogen in detritus in the water and sediments. The difference of the input and the total nitrogen in the marine system increases slowly until mid-July.
The pool of dissolved inorganic nitrogen decreases due to the spring bloom of diatoms in May. Diatom biomass remains on an elevated level after bloom peak, compared to NOSHIP scenario, due to additional nitrates in the surface layer (Fig. 3a,c). With a short delay, the pool of nitrogen bound to detritus grows as phytoplankton starts to decay. Decomposition of organic matter leads to an increase of the ammonium pool, but as ammonium is rapidly oxidized to nitrate and molecular nitrogen, the concentration of ammonium-bound nitrogen remains low. Until mid-July when cyanobacteria start blooming, the content of inorganic nitrogen stays stable in the water, as well as nitrogen bound to diatoms, while nitrogen content in the sediments increases. In July the share of nitrogen in phytoplankton increases slowly due to the uptake of nitrogen by increasing flagellates. In principle, the nitrogen fluxes during that period could be described as concurrent uptake of inorganic nitrogen by phytoplankton, their death (formation of detritus) and decomposition of organic nitrogen into NH 4 as a short-lived compound, and finally transformation into molecular nitrogen and dissolved inorganic nitrogen again via denitrification.
The Baltic Sea as a nitrogen-limited ecosystem gives an advantage to the phosphorus-limited nitrogen fixing cyanobacteria (Granéli et al., 1990). Hence, additional nitrogen input due to shipping results in lower cyanobacteria biomass and an overall decrease of total excess nitrogen in the marine system. The decline in cyanobacteria biomass results in a decrease of the amount of dead organic matter and detritusbound nitrogen as well as a negative contribution to molecular nitrogen. The share of nitrogen in the other functional groups of phytoplankton and inorganic nitrogen in the water is also declining, but less strongly compared to cyanobacteria. Since the end of August until November, the amount of inorganic nitrogen in the water is stable, the negative contribution of cyanobacteria is decreasing, nitrogen in phytoplankton and in detritus varies slowly around stable content, mostly due to the flagellate bloom, but total nitrogen starts to increase slowly. Roughly since November primary production ceases, concentration of inorganic nitrogen increases due to shipborne input, nitrogen in phytoplankton decreases to zero, and nitrogen in detritus gains stable positive level. By the end of the year shipborne nitrogen consists mainly of inorganic nitrogen in the form of nitrate in the water and organic nitrogen in the form of detritus in the water and in the sediments. Smaller amount of dead organic matter results in lower denitrification in the NOSHIP case compared to the SHIP case and negative excess gaseous nitrogen.
Considering all the fluxes as explained above, the 5-year nitrogen balance is depicted in Fig. 5b. The difference in the atmospheric fluxes consists of shipping related deposition of nitrate and ammonium. Nitrogen fixing is lower with SHIP scenario due to smaller amount of cyanobacteria in the system and has negative excess flux. Also, denitrification is higher in the SHIP case due to the internal dynamics of the marine system and the excess flux is negative, which means that more nitrogen is removed from the system in the case of SHIP scenario. The difference in the nitrogen input to the sea and excess nitrogen is more than twofold after one year of shipping activity (Fig. 5a). When we extend our simulation to a five-year period, repeating the same annual ship input and hydrophysical conditions, we obtain that the total annual amount of nitrate, i.e. sum of nitrate in the water and nitrogen in detritus approaches a steady state. The ship input of excess nitrogen is largely compensated by decreasing nitrogen fixation and increasing nitrogen removal due to denitrification compared to the NOSHIP case.
The phosphorus content shows a decrease of the phosphate in the water column, i.e. excess phosphate is negative, but an increase of the phosphorus pool in the sediments (Fig. 5c). The phosphorus input to the Baltic is mainly eliminated by binding them in the sediments. In the Baltic Proper, where much of the sediment is oxygen-deficient, sediments are an inefficient sink, indicating that phosphorus, that is already in the system, is eliminated very slowly (Savchuk and Wulff, 2009). Continuous decomposition of organic material results in higher oxygen consumption and a negative trend of the oxygen content which decelerates in time (Fig. 5d).
In absolute values, an annual nitrogen input of 20 kt from shipping, is comparable with the nitrogen input of a big river of the Baltic, e.g. the Neva river (Stålnacke et al., 1999). Hence, it could also be assumed that ship-nutrient reductions have a similar response in the ecosystem as the reduction of nutrients in rivers, which has been more thoroughly studied in the Baltic Sea. River nutrient reduction has an effect on the ecosystem -nitrogen and phosphorus inventories decrease, but it usually takes more than two years to be significant according to simulations (Neumann et al., 2002) and even longer in natural systems, depending on their complexity (Duarte et al., 2009;Schindler, 2012). Also, a compensatory mechanism often contributes in relation to the N/P ratiocyanobacteria concentrations are known to increase after nutrient reductions (Higgins et al., 2017).
Coastal ecosystems are complex and respond differently when nutrient input is decreased, as the response trajectories are rather complex. Duarte et al. (2009) has stated that coastal ecosystems seldom return to a previous, oligotrophic level, but rather remain on a new, medium level. In general an oligotrophication process starts when nutrients are reduced -phytoplankton biomass will decrease and so will the oxygen deficient zones, until a new steady state is achieved. Well-known examples of ecosystems improving to a certain level are the Black Sea shelf area and Chesapeake bay (Kideys, 2002;Langmead et al., 2009;Lefcheck et al., 2018), where dissolved oxygen conditions have improved and phytoplankton blooms have slightly decreased with nutrient reduction and ecological communities have benefitted from that. Several studies suggest that decreasing inputs of nitrogen will not hasten recovery from eutrophication (in lakes) and may even hinder it by favoring nitrogen fixers or stimulating internal loading of phosphorus. Instead repeated adding of additional nitrogen could suppress cyanobacteria blooms and in the long-term reduce nitrogen import to the system (Schindler, 2012), which was also seen in our study. Similar relative contribution of shipborne dissolved inorganic and particulate organic nitrogen has been estimated for the southern Baltic Sea by Neumann et al. (2018c). Raudsepp et al. (2013) have shown annual decrease of nitrogen fixation due to ship nitrogen deposition by 2-6% in the Gulf of Finland. Other previous studies have estimated relative ship contribution to be around 10% for different nitrogen compounds into the Baltic Sea (Tsyro and Berge, 1998;Bartnicki and Fagerli, 2008;Bartnicki et al., 2011;HELCOM, 2005).

Conclusions
The shipping contributes about 0.3% of the total phosphorus and 1.25-3.3% of the total nitrogen input to the Baltic Sea. The amount of nitrogen directly discharged to the sea from the ships is about two orders of magnitude smaller than atmospheric input of excess nitrogen (SHIP minus NOSHIP). Excess ammonium deposition is negative.
Spatially, the Gulf of Bothnia has almost negligible deposition and discharge of nutrients due to shipping, while the North-Eastern and South-Western Baltic Sea and the Kattegat area are the most impacted. Relative impact of shipborne nutrients on the biogeochemical variables in the surface layer does not exceed 10% locally, but this is already four to eight times larger than the share of the ships in total nitrogen input to the Baltic Sea. Diatoms and flagellates have a marked increase in spatial distribution while the reduction of the spatial extent of cyanobacteria blooms is extensive and covers vast areas in the Baltic Proper as well as the Gulf of Finland and several coastal areas. The relative reduction of cyanobacteria concentration of 10% due to shipping is significant during their blooming period, whereas other phytoplankton functional groups experience notable relative changes during their low concentration period in summer. Some of the additional organic matter produced with the added nitrogen sinks to the bottom, where its decomposition consumes oxygen and increases the areas of oxygen-deficient bottoms by 50 km 2 in the deep area of the central Baltic Proper within one year of model simulation. This increase in anoxic bottom area decelerates in time.
In general, the nitrogen balance showed that shipborne nitrogen input does not result in equal increase of total excess nitrogen in the water system. In the mostly nitrogen-limited Baltic Sea ecosystem, the excess nitrate is consumed by diatoms and flagellates. The increase of diatoms and flagellates biomass leaves less phosphate in the water for the cyanobacteria that occur later in summer. Their biomass decreases, which results in less atmospheric nitrogen being fixed and brought to the marine environment. Besides, denitrification is higher due to the internal dynamics of the marine system, the excess flux is negative, which means that more nitrogen is removed from the system. Multi-year continuous input of ship-borne nitrogen does not accumulate in the marine environment, but total nitrogen reaches a stationary state. The excess phosphate in the water column is negative, but phosphorus pool in the sediments increases steadily.
In summary, nutrient input from shipping does not have a significant effect on the Baltic Sea ecosystem in terms of the ecosystem functioning as changes remain within 10%. Still, shipping accounts as an important nutrient source in the context of the Baltic Sea, and is comparable in size to a large river.