BALTIC SEA- ANALYSIS OF TWELVE BASINS USING THREE-DIMENSIONAL FLOW DYNAMICS
Guillaume Vigouroux
December 2014
TRITA-LWR Degree Project 14:17 ISSN 1651-064X
LWR-EX-2014:17
c Guillaume Vigouroux 2014
Master Degree Project Environmental Engineering and Sustainable Infrastructure Division Land and Water Resources Engineering
Done in association with the Water Resources Engineering Research group Royal Institute of Technology (KTH)
SE-100 44 STOCKHOLM, Sweden
Reference to this thesis should be as follows: Vigouroux, G. (2014) "Carbon-based nutrient cycling in the Baltic Sea - Analysis of twelve basins using three-dimensional flow dynamics".
TRITA-LWR Degree Project 14:17. 39 pages.
SUMMARY IN SWEDISH
Östersjön, som är belägen emellan Centraleuropa och norra Europa, är en av världens största bräckta vatten- drag och drabbad av övergödning. Detta beror främst på grund av en ökning av kväve och fosforbelastningar, skapade av människan. Under det senaste århundradet har algblomningen ökat och även blomningen av den blå- gröna algen, vilket orsakar olika problem för Östersjön, såsom anoxisk havsbotten och minskning av fiskpop- ulationen. Detta har i sin tur en negativ effekt på fritidsvärden, som då kommer att uppleva en förlust. De problemen har negativa ekologiska och socioekonomiska konsekvenser och därför har åtgärder tagits av HEL- COM, som är konventionen för att värna miljön i Östersjön.
Övergödning i Östersjön är beroende av många parametrar såsom belastningar av växtnäring, temperaturen och fluiddynamiksegenskaper. De parametrarna utvecklas naturligt men förandras också av människors handlin- gar, som kan vara direkt, till exempel policy, eller indirekt såsom klimatförändringen. Modellering är därför ett sätt att förutsäga effekten i politiska åtgärder och mer allmänt, konsekvenser av parametrarförändringar.
Kiirikki et al. modellen för alger, vattenkvalitet och sediment, har utvecklats för att beskriva övergödning och koncentration av växtnäring i Finska viken. Den har gett bra resultat angående blågröna algblomningar och den inre belastning som finns i sedimenten utan att införa många ovissa parametrar, och som just därför har valts för den föreliggande studien. Studiemålet är att bedöma modellen tillämpligheten för hela Östersjön.
För att implementera modellen har Östersjön delats i tolv bassänger separerade med lagergångar. De bassängerna har sammankopplats alltefter flöden och delats upp i två vertikala lager, som skiljs där densitets- gradient är maximal eftersom utbytet av växtnäring minskas. Kiirikki et al. modellen tillämpas av varje bassäng.
Modellen har kalibrerats i 2001 och 2002 och simuleringar har gjorts från 2003 till slutet av 2009.
Resultaten efter kalibrering visar att modellen överensstämmer med kontrolldata för den södra delen, även om den modellerade lösta oorganiska kväve koncentrationen är för hög under vinterperioder och några fas- förskjutningar konstateras. För den norra delen är primärproduktionen väl modellerad men den lösta oorganiska kvävekoncentrationen är för låg. Dessa skillnader med kontrolldata och andra modeller visar att modelltilläpnin- gen och kalibreringen kan förbättras. Till exempel behövs en global känslighetsanalys genomföras för att förd- jupa kunskapen av parametrarinfluens och på så sätt kan kalibrering förbättras. Ytterligare kan bassängspecifika parametrar och solstrålningar implementeras och bättre flödesdata kan användas för att få mer realistiska resul- tat. Avslutningsvis är Kiirkki et al. modellen lämplig för att beskriva övergödning och vattenkvalitet i Östersjön även om den behöver förbättras för att användas av beslutsfattarna.
ACKNOWLEDGEMENT
I would firstly like to thank the people that helped me during this work. I am grateful towards Vladimir Cvetkovic for his unconditional support and pertinent advices. I am thankful to Anders Jönsson and Benoît Dessirier for their help and great ideas that made this work truly interesting. I also have to thank Sofie Soltani for her help and support in the office. Finally, I am appreciative to Prabin Paul and Bijan Dargharli for their help.
I also want to thank the people helping indirectly to the realisation of that work. I am greatly thankful to my parents, who always support me. Thank you as well to Mengni, Liwen, Xavi, Yuanying, Tushar, 9gag for making the days in the office interesting and fun. Thanks to Yuki, Civeny and Sanna and to all my friends for those pleasant more than two years at KTH. Last but not least, big up to Bingcheng.
TABLE OFCONTENTS
SUMMARY INSWEDISH. . . . III
ACKNOWLEDGEMENT . . . . V
ABBREVIATIONS ANDGLOSSARY. . . . IX
ABSTRACT. . . . 1
1 INTRODUCTION. . . . 1
1.1 Baltic Sea. . . . 1
1.2 Eutrophication . . . . 2
1.3 Eutrophication models. . . . 2
1.4 Overall plan . . . . 4
2 METHODS. . . . 4
2.1 Global description . . . . 4
2.2 Kiirikki et al. model . . . . 4
2.3 Hydrodynamics. . . . 5
2.3.1 Blocks definition . . . . 5
2.3.2 Layer separation . . . . 7
2.4 Implementation of the model . . . . 7
2.4.1 Resolution method. . . . 7
2.4.2 Data . . . . 10
2.4.3 Assumptions and uncertainties . . . . 11
2.4.4 Model set-up . . . . 11
2.4.5 Calibration and validation . . . . 11
3 RESULTS AND DISCUSSION . . . . 12
3.1 Comparison between basins. . . . 13
3.2 Comparison with monitoring data . . . . 13
3.3 Comparison with other model results . . . . 19
3.4 Limitations and extensions of the model. . . . 19
4 CONCLUSION. . . . 21
REFERENCES . . . . 22
Data . . . . 23
Other references . . . . 23
APPENDIXA: ODESYSTEM FROM THEKIIRIKKI ET AL.MODEL. . . . 24
A.1 Equations . . . . 24
A.2 Functions . . . . 24
A.2.1Rates . . . . 24
A.2.2Limiting factors . . . . 25
A.3 Model variables . . . . 25
A.4 Parameters . . . . 25
APPENDIXB: SUB-BASINS PROPERTIES AND MONITORING STATIONS. . . . 28
LIST OF FIGURES
1 Map of the Baltic Sea . . . . 1
2 Conceptual model of the eutrophication . . . . 3
3 Conceptual model of the Kiirikki et al. eutrophication model. . . . 6
4 Map of the sub-basins of the Baltic Sea. . . . 8
5 Depth of the pycnoclyne . . . . 9
6 Location of the results and monitoring station. . . . 14
7.A Simulation results for the Eastern Gotland Basin for 2002. . . . 15
7.B Simulation results for the Eastern Gotland Basin for 2002. . . . 16
8 Simulation results for the Bothnian Bay . . . . 17
9 Simulation results and monitoring data for the Eastern Gotland Basin . . . . 18
LIST OF TABLES 1 MASE before and after calibration. . . . 12
2 MASE for the validation period (2003-2009). . . . 20
3 Comparison of the primary productivity. . . . 20
ABBREVIATIONS AND GLOSSARY
3-D Three-dimensional
C Carbon
CE-QUAL-ICM Three-dimensional eutrophication model for the Chesapeake Bay and major tributaries COHERENS Coupled Hydrodynamical Ecological model for Regional Shelf seas
COMBINE Cooperative Monitoring in the Baltic Marine Environment
det detritus
DIN Dissolved Inorganic Nitrogen
DIP Dissolved Inorganic Phosphorus
diapycnal mixing Mixing across surface of constant density EFDC Environmental Fluid Dynamics Code
FVCOM Finite-Volume, primitive equation Community Ocean Model GEMSS Generalized Environmental Modeling System for Surfacewaters GFDL Geophysical Fluid Dynamics Laboratory
GIS Geographic information system
halocline Vertical zone in which salinity changes rapidly with depth
HELCOM Baltic Marine Environment Protection Commission (Helsinki Commission) ICES International Council for the Exploration of the Sea
MASE Mean Absolute Scaled Error
MOM Modular Ocean Model
N Nitrogen
P Phosphorus
POM Princeton Ocean Model
pycnocline Vertical zone in which density changes rapidly with depth
Ri Richardson number
SMHI Swedish Meteorological and Hydrological Institute SHARK Swedish maritime archives (Svenskt HavsARKiv) STRÅNG Mesoscale model for solar radiation
thermocline Vertical zone in which temperature changes rapidly with depth UMASSD University of Massachusetts Dartmouth
WHOI Woods Hole Oceanographic Institution
ABSTRACT
Eutrophication is a major problem in the Baltic Sea and is the result mainly of the increase of the anthropogenic nutrient loading. Thus, the links among water quality, sediments, and eutrophication have to be understood in order to predict the consequences of our actions and of the climate change on the Baltic Sea. Therefore, water quality models that take into account the hydrodynamics have to be developed to help policy makers.
In that perspective, the Kiirikki model, an ecosystem and sediment model, coupled to a box approach has been used to describe the water quality of the Baltic Sea. The latter has been divided into twelve sub-basins according to the topography, each of them separated into two vertical layers. The Kiirikki model has been implemented on each sub-basin and the hydrodynamics are used to link sub-basins between them.
After calibration, it can be seen that the model results are consistent with the monitoring data for the southern part, even if the dissolved inorganic nitrogen levels are too high during the winter and some phase shifts are observed. For the northern part, the primary production is well modelled but there is an offset concerning the dissolved inorganic nutrient. Thus, it can be concluded that the implementation of the Kiirikki model is a realistic tool to describe the water quality and eutrophication of the Baltic Sea. However, the differences indicate that the Baltic Sea model cannot be use for policy making yet and more work is needed to improve the model such as a global sensitivity analysis as well as the use of site specific parameters.
Keywords: Eutrophication; Baltic Sea; Hydrodynamics; Nutrients; Modelling; Sediments; Ecosystem.
1. INTRODUCTION
1.1. Baltic Sea
The Baltic Sea is situated between Central and Northern Europe. With an area of 377, 000k m2, it is a relatively small sea but one of the largest brack- ish water bodies in the world. It is enclosed by nine countries: Denmark, Germany, Poland, Lithuania, Latvia, Estonia, Russia, Finland, and Sweden (Fig. 1).
Its catchment area is 1, 641, 650k m2 and inhabited by over 85 million people (SMHI, 2004; HELCOM, 2010a). The Baltic Sea is connected to the North Sea only through the Danish Straits, which, with its vast catchment area, makes it easily vulnerable to pollu- tion.
The Baltic Sea is heterogeneous and has a complex geometry; it is constituted of straits, sills, archipelago and open sea areas and therefore can be separated in natural sub-basins, each of them having different physico-chemical and biological properties (HEL- COM, 2010a; Omstedt and Axell, 2003). Moreover, the inflow of fresh water through river discharge is superior to the inflow of saline water from the North Sea, which makes most of the part of the Baltic Sea strongly stratified, with a permanent halocline at about 60m, for the central part (Reissmann et al., 2009). Therefore, the hydrodynamics properties of the Baltic Sea are complex and difficult to model.
The Baltic Sea plays an important socio-economic role for the surrounding countries. A lot of indus- tries depend on the Baltic Sea, for example for the transportation of goods and passengers to and be- tween the different countries or through the provi- sion of fishes by fisheries and fish farms. The Baltic
Sea supports also some tourism activities, in particu- lar in the archipelagos and is also important for the population, as it provides multiple recreational activ- ities as fishing or swimming. Most of those activities depend on ecosystem services, and therefore on the environmental health of the Baltic Sea (HELCOM, 2010a). Thus preserving or improving the quality of this ecosystem is socially and economically impor- tant.
Fig. 1: Map of the Baltic Sea (source: Wikipedia).
Due to its particularities, the Baltic Sea possesses a unique ecosystem with a low species richness but genetically deviant from other populations. This is due to difficult physical conditions, for marine
and freshwater species, with high salinity variations (from 2 to 25P SU ) and low winter temperatures (Johanesson and André, 2006). This ecosystem is threatened by multiple anthropogenic pressures such as the release of pollutant, not only by the large coastal cities but due to the large and densely pop- ulated catchment area, the overexploitation of some maritime resources or the maritime transport (HEL- COM, 2010a).
One of the most important anthropogenic pres- sures is the increase of the amount of nutrients (phos- phorus and nitrogen) in the Baltic Sea, which causes eutrophication and changes the whole dynamic of the Baltic ecosystem. Those nutrients come from different sources, such as agriculture, transported by the rivers or leaching directly in the sea, industries, wastewater treatment plant, maritime transport or airborne emissions (HELCOM, 2010a). The nu- trients are after transported by the water current and taken up by planktons. Due to the complex- ity of the hydrodynamics, the diversity of the nu- trients sources and its potentially important conse- quences, eutrophication constitutes one of the main challenges of the Baltic Sea.
1.2. Eutrophication
The main trophic states characterizing a water body are the oligotrophic state, which has low pri- mary production (algal biomass), low nutrient con- tent and deep photic zone, theeutrophic state, which has high primary production, high nutrient content and important variation in the dissolved oxygen con- centration, and themesotrophic state, with interme- diate properties (Dodds et al., 1998). Eutrophication can be defined as "an increase in the rate of supply of organic carbon to an ecosystem" (Nixon, 1995).
Here, it is considered to be the ecosystem response due to an increase of the nutrient resources, leading to a change of the trophic status.
Even though eutrophication is a natural process, its causes are most of the time of anthropogenic na- ture, as human activities have globally increased the release of nutrients. In the Baltic Sea, the main exter- nal sources are the agriculture, with the use of fertil- izer that are leached from the soil and transported by the rivers or directly to the sea, the wastewater treat- ment plants, releasing nutrients from the population and industry, some other industries, such as the fish farms and the atmospheric deposition, enhanced by the shipping traffic (HELCOM, 2013b). Nutrient loads have increased importantly during the second half of the 20th century, about four times for the ni- trogen and eight for the phosphorus (Larsson et al., 1985). Moreover, the increase of organic matter due to eutrophication can deplete the bottom water of oxygen during its decomposition and thus creating anoxia events, which can also cause an internal load-
ing, by releasing phosphate from the nutrient by a re- duction of oxidized iron compounds (Kiirikki et al., 2001).
By changing the ecosystem and the relative dis- tribution of the species, eutrophication creates im- portant ecologic as well as socio-economic problems.
The main ecological problems are the increase of the biomass of phytoplankton, the change in the species composition, with the increase of cyanobacteria al- gae, of benthic and epiphytic algae and of the biomass of consumer species, and the loss of biodiversity, with the reduction of the fish population. These problems have economic and social consequences with the diminution of the harvestable fish popu- lations and the loss of recreational value, through toxic algae blooming, diminution of the water trans- parency and of the aesthetic value of the water body (Smith and Schindler, 2009). However, this increase of primary productivity can be considered as an ad- vantage for the production of fuel, for example.
The eutrophic state is a stable and resilient state of the ecosystem and it is therefore quite difficult to come back to an oligotrophic state. Thus, in the Baltic Sea case, it is important to take measures to avoid a shift to the eutrophic state. HELCOM, which is the Baltic Marine Environment Protection Commission, works on the protection of the marine environment of the Baltic Sea, defined different ob- jectives for reducing the eutrophication (HELCOM, 2008; HELCOM, 2010a):
• clear water,
• concentrations of nutrients close to natural lev- els,
• natural level of algal blooms,
• natural distribution and occurrence of plants and animals,
• natural oxygen levels.
1.3. Eutrophication models
Eutrophication models are a type of ecosystem models, therefore, they have to consider the influ- ence of the relevant physical parameters, such as the nutrient concentrations, the temperature, as well as the relations between algae and their ecosystem and thus their place in the food web (Fig. 2) (HELCOM, 2006).
However, considering all the relations with the food web complicates the model as it requires cal- ibrating more parameters and thus more data is needed and asks for more computational power.
Therefore, most of the models do not consider the fish influences, such as the CE-QUAL-ICM by Cerco and Cole (1993) or the model developed by Fennel (1995) and others do not consider the zooplankton influences either (Tyrrell, 1999). Such models sim- plify the interaction with the food web in a growth
Fig. 2: Conceptual model of the eutrophication (source: HELCOM, 2006).
and decay coefficient and can still realistically de- scribe the eutrophication.
Meteorological conditions and the hydrodynamic response of the sea, such as currents and vertical mix- ing, play a crucial role in the modelling of eutroph- ication, especially in systems that cannot be consid- ered as fully mixed, as they influence the availability of nutrients and the transport of algae (Neumann, 2000). Therefore, hydrodynamic modelling has to be integrated for modelling the eutrophication. This can be done either by using external hydrodynamic simulations that are then incorporated in the ecosys- tem model, or by integrating the ecosystem model with hydrodynamic simulations.
There are multiple hydrodynamic models that can be used to simulate the physical properties of the water body of interest. Many of them derive from Princeton Ocean Model (POM) developed by Blum- berg and Mellor (1987), such as FVCOM, a finite vol- ume 3-D primitive equation coastal ocean circulation
model developed by UMASSD and WHOI joint ef- forts. Other models that build on POM, are also used such as GEMSS, a 3-D hydrodynamic, transport and water quality model, developed at J. E Edinger and Associates. (Hubertz, 2014), COHERENS, a three-dimensional multi-purpose numerical model, designed for application in coastal and shelf seas, es- tuaries, lakes, reservoirs (Luyten, 2013) or MOM, the modular ocean model developed by the GFDL.
Some hydrodynamic models integrate an ecosys- tem or eutrophication model. This is the case of the Environmental Fluid Dynamics Code (EFDC), based on the POM that includes an eutrophication module, or the Delft3D, which is a modelling suite simulating flows but also water quality and ecology (Hubertz, 2014).
Modelling studies on the eutrophication of the Baltic Sea have been carried out by coupling hydro- dynamics and eutrophication models. The Swedish Coastal and Ocean Biogeochemical model has been
coupled to the Rossby Centre Ocean model (RCO- SCOBI) and the Neumann ecosystem model coupled to MOM have been used to both model and predict the climate change effect on the Baltic Sea (Meier et al., 2011; Neumann, 2010). The Kiirikki et al. model has also been used to model the eutrophication in the Gulf of Finland, built on top of a 3-D water quality model (Kiirikki et al., 2001).
1.4. Overall plan
Eutrophication has multiple consequences on the ecosystem, most of them being negative from an ecological and socio-economical point of view. Eu- trophication of the Baltic Sea depends on a lot of pa- rameters such as the nutrient loads, the temperature, and the hydrodynamic properties. Those parame- ters evolve naturally, but are also modified by anthro- pogenic actions, that can be direct such as policy or indirect, through climate change. It is therefore im- portant to have tools allowing us to predict the ef- fectivity of policy measures but also, more generally, the consequences of change in the parameters.
The Kiirikki et al. model has been developed for describing the eutrophication, water quality and sed- iment conditions of the Gulf of Finland, without in- troducing a multitude of uncertain parameters, and has given realistic results (Kiirikki et al., 2001; Ki- irikki et al., 2006). Moreover, it has also been shown to give the best results in regard of the modelling of the internal loading process, and is thus chosen for the present study (Dessirier and Soltani, 2011). The goal of this study is to use this model on the whole Baltic Sea, partitioned in different blocks, by cou- pling hydrodynamic data to the water quality model for a large scale domain, and to assess if this approach can be used to simulate the eutrophication of the Baltic Sea, and to which extent such simulations are comparable to measurements.
In order to assess the performance of the model, a simulation will be carried from April 2001 to the end of 2009 and will be compared with monitoring data as well as other model results in Section 3. In Section 2 the methods and approach used to obtain the results will be described.
2. METHODS
2.1. Global description
For this study, the modelling approach has been to implement the ecosystem and sediment models de- veloped by Kiirikki et al. (Kiirikki et al., 2001; Ki- irikki et al., 2006) on the different sub-basins of the Baltic Sea.
The eutrophication and water quality are de- scribed by the Kiirikki et al. model, which is a sim- plified version of an ecosystem model, aiming to de- scribe growth of two groups of algae, cyanobacteria
and other phytoplankton and their relation with the nutrients under different forms, and taking into ac- count only the necessary physical parameters. This model is described in detail in Section 2.2.
As the Baltic Sea is not homogeneous, but rather constituted by different sub-basins separated by sills, this separation has been taken into account for the implementation of the Kiirikki et al. model. The separation used has been made according to the COMBINE program of HELCOM and is described in Section 2.3.1. Each sub-basin is then divided into two layers, one representing the surface layer and the other the deeper part and a volatile sediment layer is also taken into account. The method for the sep- aration is explained in Section 2.3.2. Therefore the hydrodynamics described by Dargahi and Cvetkovic (2014) is essential in this study, for the transport among the sub-basins and the delimitation of the ver- tical layers as well as for governing the biogeochemi- cal processes.
The implementation of the model, the main as- sumptions and the data used are presented in Section 2.4.
2.2. Kiirikki et al. model
The ecosystem part of the model calculates the concentration of two competing groups of algae:
the cyanobacteria (principally the genera Nodularia and Aphanizomenon) and the other phytoplankton (constituted of the diatoms and the flagellates), as well as the concentration of different variables driv- ing the dynamics of those algae, such as the dissolved inorganic nutrients (DIN and DIP) and the nutri- ents, as well as the carbon in detritus (Ndet, Pdet and Cdet). This model has been developed for the Gulf of Finland. It does not include detailed food web re- lations or oxygen concentrations and therefore does not aim to be a complete ecosystem model; however it is a realistic cyanobacteria model (Kiirikki et al., 2001).
This ecosystem model is based on the oceanic phytoplankton model of Tyrrell (1999) and follows different principles, as described by Kiirikki et al.
(2001):
• The phytoplankton is constituted of two com- peting groups of algae: the nitrogen fixing cyanobacteria and the other phytoplankton.
The other phytoplankton is more efficient when DIN and DIP are available. When DIN is low and DIP still available, the nitrogen fix- ers are more competitive. However, the nitro- gen fixing process requires energy and is thus not efficient at low temperatures. Nutrients are consumed by both algal groups according to the Redfield ratio (Redfield, 1958).
• The losses of biomass are described by a tem- perature dependent loss term, lower for the
cyanobacteria, as they are avoided by the grazers due to their toxicity (Sellner et al., 1996). The growths of both algae are dependent on the solar radiation, depending also on the presence of ice, and nutrient availability according to Michaelis- Menten kinetics and on the temperature accord- ing to Frisk (1982). The growth is also limited by a self-shading factor, representing the carry- ing capacity of the system. To avoid extinction during the winter period, a minimal biomass, for which the loss rate approaches zero, is de- fined.
• The dead biomass of algae is injected in the detri- tus according to the Redfield ratio. The detritus are sinking and also mineralized into dissolved inorganic nutrients. When they reach the bot- tom, they are converted into sediments.
The sediment part of the model is implemented according to the description made by Kiirikki et al.
(2006). The sinking detritus coming from the decom- position of the algae settle as volatile sediments when they reach the bottom. Under aerobic conditions, denitrification can take place and some of the phos- phorus in sediment is bound to ferric iron, and thus only some of the mineralized nutrients are released from the sediments. This is not the case under anaer- obic conditions and all of the mineralized nutrients are released (Kiirikki et al., 2006).
Phosphorus bound to ferric iron (F e(I I I )) is also taken into account. Ferric iron occurs in sediment and has phosphorus binding abilities, which are lost when it is reduced to ferrous iron (F e(I I )) (Gunnars et al., 2002). Under anoxic conditions, microbial ac- tivities are driven by the reduction of ferric iron or of sulfate (SO42−), leading also to the reduction of ferric iron due to the formation of hydrogen sulfide (H2S), and producing also large flux ofC O2(Kiirikki et al., 2006).
Under aerobic conditions, the ratioD I ND I P should re- main close to the Redfield ratio. This implies that the denitrification rate and the binding of phosphorus to ferric iron should be equal. On the contrary, un- der anaerobic conditions this ratio is lower than the Redfield ratio, thus the release of iron bound phos- phorus is higher than the mineralization of nitrogen (Kiirikki et al., 2006).
However anoxic conditions are hard to describe with the oxygen concentration as the gradient of oxy- gen can be important near the bottom due to the mineralization of organic matter. Settled organic matter is mainly oxidized toC O2by organisms and therefore, the flux ofC O2 is used by the model to determine anoxic events. When the flux ofC O2 is low, oxygen consumption is also low and sediments are able to retain part of the mineralized nitrogen and phosphorus, which is not the case when the flux is high. The limit between aerobic and anaer-
obic conditions is for aC O2 flux of approximately 200− 470 m g .C .m−2.d ay−1(Kiirikki et al., 2006).
This model has given realistic results for the wa- ter quality simulations in the Gulf of Finland, but is subject to some limitations. Some parameters are site specific for the Gulf of Finland, for example, con- cerning the sediment part or the growth rates, and their value may differ for the Baltic Sea. Moreover, the limitation due to trace elements for the growth of cyanobacteria is not taken into account for the Gulf of Finland but might apply in the Baltic Sea (Ki- irikki et al., 2001). More importantly, by not con- sidering the complete food web relations, this model could lose its validity in the case of drastic food web changes. However, considering an important part of the food web would introduce a multitude of uncer- tain parameters and thus increase the model’s uncer- tainty.
The eutrophication model is conceptually repre- sented in Figure 3.
2.3. Hydrodynamics 2.3.1. Blocks definition
The Baltic Sea is relatively shallow, with a mean depth of 54m and presents important variation of its depth with the presence of sills and deeper parts.
Decomposition into sub-basins can be made accord- ing to those sills and is justified as the flows between them are largely governed by the sills (Leppäranta and Myrberg, 2009).
The current inside and between the blocks are of two types. The surface currents are, on a short time scale, mainly caused by winds, but the influence of winds is weak on the long term, due to their high variability. On a longer time-scale (some months), a density driven (baroclinic) circulation appears due to the positive fresh water budget. The water is mixed by mesoscale eddies and deep water circulation (Lep- päranta and Myrberg, 2009). Deep water circulation is created by inflow events of saline water from the North Sea. They can be barotropic (sea level differ- ences) or baroclinic (difference of salinity) and occur irregularly. Major Baltic Inflows occur at a decadal time, ventilate the deep water layer with oxygen rich water and affect the whole stratification (Reissmann et al., 2009).
The chosen block definition (Fig. 4) was made by HELCOM for the COMBINE program, which define a common framework for the monitoring of the Baltic marine environment (HELCOM, 2013a).
This block definition is interesting as it follows the topography of the Baltic Sea: the sub-basins are sep- arated by the sills and therefore the mixing pro- cesses are more important inside each sub-basin than among the sub-basins. The only open boundary condition is situated in the south west of the Baltic
Fig. 3: Conceptual model of the Kiirikki et al. eutrophication model.
Sea and take into consideration the influence of the North Sea. Moreover, as this definition is based on the COMBINE program, it is part of a framework used for the management of the Baltic Sea.
2.3.2. Layer separation
Vertical mixing plays a crucial role in the dynam- ics of marine ecosystem, by participating in the trans- port of nutrients and making them available for algal growth. This vertical mixing is carried out by dif- ferent processes. Inflowing water from the North Sea creates strong bottom current, participates to the mixing by entrainment and creates an uplift of the water column. Inertial waves and internal wave breaking enhanced the vertical turbulent mix- ing. Coastal upwelling can transport sub-halocline water in the surface water and surface cooling dur- ing the winter induce unstable stratification leading to convection of the water (Reissmann et al., 2009).
These diapycnal mixing processes, by allowing the deep phosphorus to reach the surface layer, play a major role in the eutrophication. Therefore, the py- cnocline is an interesting parameter for the layer sep- aration in the Baltic Sea.
Density is an important factor for mixing as the mixing decreases when the density difference in- creases (Leppäranta and Myrberg, 2009). The den- sity of water depends mainly on the temperature, the salinity and the pressure. The gradient of salinity in the Baltic Sea is quite strong as the inflow of fresh wa- ter from rivers and precipitation is bigger than the in- flow of salty water from the North Sea (Leppäranta and Myrberg, 2009). The halocline is usually found at approximately 60 meters. During summer, a shal- low thermocline is formed and protects the halocline against erosion, therefore making it relatively weak.
On the contrary, there is no thermocline during the winter and the halocline is stronger (Reissmann et al., 2009). The pycnocline, which is defined as the zone where the gradient of density is the highest, depends therefore on the variations of temperature and salin- ity and varies seasonally (Fig. 5) (Reissmann et al., 2009). This pycnocline can be used for the definition of a surface and a deep layer as the vertical mixing and the transport of nutrients through the pycnocline is reduced.
The Richardson number (Ri) is the ratio between the stabilizing buoyancy forces and the destabilizing shear production of turbulence:
N2
∂ U
∂ z
2
whereN is the buyoancy frequency and U is the ve- locity (Howard et al., 2004).
Values of Ri < 0.25 indicate that shear instabili- ties and therefore vertical mixing are likely to occur
and values ofRi< 1 that advective instabilities could occur (Howard et al., 2004). High values ofRi indi- cate a strong stratification and little dispersive mix- ing. Therefore it is an important parameter for the determination of vertical mixing.
The Secchi depth is also useful for the determina- tion of vertical layers as it provides an estimation of the depth of the euphotic zone: it is approximatively twice the Secchi depth (Aarup, 2002).
Layer delimitation method
Each block is separated in two fully mixed, homo- geneous layers. The separation is done at the pycno- cline, which is calculated for each column constitut- ing the block and then averaged over the area where there is no mixing. Mixing is considered when the maximal gradient is less than a critical value (taken at 0.1) or the Richardson number is less than 0.25 (winter and coastal uppwelling). If there is mixing for more than half of the area, the block is consid- ered to be constituted of only one layer, fully mixed.
The phytoplankton is supposed to grow only in the surface layer, at a maximal depth of 10 me- ter, which correspond to an approximation of the depth of the euphotic layer for the Baltic Sea (Müller- Karulis et al., 2012). The exchange of nutrients be- tween the surface and the deep layers are governed by the vertical velocity determined by the hydrody- namics. Another method could be based on the salin- ity gradients to estimate vertical turbulent exchanges (Howard et al., 2004).
2.4. Implementation of the model 2.4.1. Resolution method
The eutrophication and water quality model has been implemented with Python as it is a high-level, object oriented and interpreted programming lan- guage (Python Software Foundation, 2014). It is composed of two main parts: a preprocessing part that prepares the data used in the simulation part, which runs the model developed by Kiirikki et al.
There is also a post processing part used for the treat- ment and the graphical representation of the results.
Preprocessing
The preprocessing is done in two steps. In the first step, the information used for the definition of the problem is retrieved from the hydrodynamic data, such as the time step used for the hydrodynamic sim- ulation, the cells constituting the grid and the depth of the different layers. The different blocks are also defined as a collection of cells and the interfaces, which are the area between two blocks in contact, is defined as a collection of points. The nutrient loads for the different blocks for each day of the simulation period are also extracted from the load files.
Fig. 4: Map of the sub-basins of the Baltic Sea.
Fig. 5: Depth of the pycnoclyne, at the 21/06/2001 (left) and 21/12/2001 (right).
Then, the hydrodynamic data is extracted and pro- cessed for each simulation time step. The hydrody- namic data is described in details in Section 2.4.2.
Each time step, the velocity vector, the tempera- ture, the salinity, and the Richardson number are ex- tracted, for each point of the grid and vertical layer.
They are then converted into cell based values by av- eraging over the four points constituting a cell. The water density of each vertical layer of the cell is calcu- lated using the formula given by Millero and Poisson (1981)ρ = ρ0+AS+BS1.5+C S2, whereρ0is the den- sity of pure water,S is the salinity and the coefficients A, B and C depend on the temperature. The depth of the pycnocline is calculated for each cell and the verti- cal layer definition for each block is made according to the method described in Section 2.3.2. For each cell, the volume, the area and the volume averaged temperature and salinity for each layer, the upward and downward flow between the layers and the tem- perature and velocity at the bottom, are determined.
To get the physical properties of the block, those val- ues are aggregated over the cells constituting it. The flows at the interfaces are also calculated as the nor- mal part of the velocity vector for each point multi- plied by the area corresponding to each point. The flows going in different directions of the interfaces are separated and each of these flows is constituted of four components representing the different possibil- ities for the flow to go from the surface or deep layer or the surface or deep layer. As the pycnocline sepa- ration can be different for each cell, all four compo-
nents can be non-null. Finally a maximal acceptable time step is calculated for the interface, which corre- sponds to the minimal time needed for all the water of one layer of the blocks to be changed.
The results are then saved in an hdf5 file (hierar- chical data format) in order to have a good organiza- tion and portability of the data. The solar radiation and boundary conditions are also incorporated to the file.
Simulation
The simulation part solves the Kiirikki et al.
model for each block for the time of the simulation.
In order do so, the time is discretized into time step of one day and the model variables are separated in three categories:
• The ecosystem variables, which are constituted of the cyanobacteria and the other phytoplank- ton. They are only present in the surface layer of the block.
• The nutrient variable, which are constituted of the DIN, DIP, as well as the nitrogen, phospho- rus and carbon in detritus for both the surface and the deep layer.
• The sediment variables, constituted of the nitro- gen, phosphorus and carbon in the volatile sed- iments and the iron bound phosphorus.
The simulation part starts by reading the file cre- ated during the preprocessing in order to initialize
the blocks and interfaces with their physical proper- ties and initial condition. A restart is also possible to set the initial conditions as the last value of the restart file. A vector of simulation dates with a time step of one day is also created as the simulation time step can be different of the hydrodynamic one.
After the initialization, the solving of the system is constituted of three steps:
• For each hydrodynamic time step, the depth of the separation between the vertical layers is changed and so the concentration of the nutri- ent variables are changed, as well as the vol- ume of the layers, to take into consideration the changes in the pycnocline depth.
• Then for each simulation time step, the flow be- tween the blocks is considered. The interfaces are first sorted by increasing maximal acceptable time step and this time step is considered for the exchange of water between the blocks of the interface so that no more than all the water of one layer of the block is exchanged during this step. The surface variables are exchanged be- tween the blocks according to the surface flows, while the nutrient variables are exchanged be- tween the layers of the blocks according to the four components of the flow in both directions.
The sediment variables are not changed during this step.
• Finally the equations representing the Kiirikki et al. model are implicitly solved for each block for a time step of one day. In a first time the parameters and functions used in the equations are calculated and then the equations are solved using a Python solver. The equations are given in the Appendix A.1.
The separation of those two steps is based on the CE-QUAL-ICM eutrophication model as it uses also a block approach (Cerco and Cole, 1994).
During the simulation, the variables are also saved in a hdf5 file for each time step so they can be repre- sented by the post-processing part.
2.4.2. Data
In this part, the different datasets used and their processing in order to extract the useful information are presented.
Hydrodynamic data
The hydrodynamic data used by the model corre- spond of the results of a GEMSS simulation carried out B. Dargahi for the SEABED project, for the pe- riod between 04/04/2001 and 31/12/2009 (Dargahi and Cvetkovic, 2014). This dataset was constituted of three types of file:
• HDM files that represent the time in Julian day
of the model outputs.
• GRD files that represent the depth of the verti- cal layers and the coordinate and indexes of the grid points used in the model output.
• CTM files containing the hydrodynamic re- sults for each date and each point of the grid.
Those results are instantaneous and considered as punctual. They are constituted of the ve- locity vector, the temperature, the salinity, the Richardson number, as well as other parameters not used such as the shear stress or the vorticity.
These files were then used by the preprocessing part to extract the useful information. Section 2.4.3 dis- cusses the approximations that are implied by the used of this instantaneous and spatially discrete data.
Block delimitation data
The map showing the block definition is given in Figure 4. The dataset used for the delimitation of the block has been downloaded from the HEL- COM map service and is a GIS shapefile representing the blocks used by the COMBINE program (HEL- COM, 2013c). The block definition has been made by using a GIS software from this GIS shapefile and the hydrodynamic file representing the grid. From this operation, a text file attributing each cell of the grid is given and used by the preprocessing part for the definition of the blocks.
Monitoring data
The monitoring data has been obtained from the SHARK database of SMHI and the ICES database (SMHI, 2013; ICES, 2014). Two or three outer sea monitoring stations have been considered for each block, if available. The list of the stations used and their corresponding blocks is given in Appendix B.
The monitoring data was constituted of the posi- tion of the measurement, the depth, the date and the nitrite, nitrate, ammonium, phosphate and chloro- phyll concentrations. From this information, the DIN, DIP and algal density were computed for each layer of the blocks.
This monitoring data has also been used to create a boundary condition vector with a station situated at the South-west part of the Baltic proper.
Loads data
The nutrient loads have been calculated using the Fifth Baltic Sea pollution load compilation from HELCOM (PLC-5.5), which gives us the total ni- trogen and phosphorus input per year between 1994 and 2010 (HELCOM, 2013b). The loads values were given by country and by sub-basin of the Baltic Sea, each sub-basin corresponding to one or more sub- basin defined by the COMBINE program. There are three types of loads: the atmospheric deposition, the point source and the riverine loads.
The atmospheric deposition is supposed to be con- stant throughout the year (Rolff et al., 2008). For the aggregated sub-basins, the separation into the sub- basin used by the model was made proportionally to the area.
As the point sources are releases of nutrients by industries and waste water treatment plants, they de- pend mainly on the human activity and were there- fore also considered as constant throughout the year.
A dataset from HELCOM with the position of the point source and their loading value was found for the year 2006 (HELCOM, 2010b). Those loads were separated by blocks and the separation into the de- fined sub-basin was made proportionally to the load values calculated.
The riverine loads were considered to be propor- tional to the river flow because they correspond mainly to soil leaching that takes place with the pre- cipitation. The daily river flow data was obtained by hydrologic simulation for the SEABED project for most of the river of the Baltic Sea. The flows were aggregated by country and by blocks and the concen- trations were calculated by dividing the load values by the yearly discharges. Those concentrations of ni- trogen and phosphorus were then multiplied by the daily discharges to get the loads.
The total loads have been converted into DIN and DIP assuming that 90% of the total nitrogen is DIN and 75% of the total phosphorus is DIP.
Solar radiation
The solar radiation data was obtained from the SMHI STRÅNG database (SMHI, 2014). Those re- sults were obtain with the STRÅNG model, which is a mesoscale model calculating the solar radia- tion. The daily global irradiance was extracted from the database inW .m−2.d ay−1and converted in M J .m−2.d ay−1. This data was extracted for the sim- ulation time (between 04/04/2001 and 31/12/2009) at the location of Stockholm.
2.4.3. Assumptions and uncertainties In addition to the assumptions made by the con- ceptual model of Kiirikki et al., the modelling ap- proach used and the data are also subject to some as- sumptions and uncertainties.
A first assumption of the modelling approach is made by the separation of the Baltic Sea into sub- basins and layers. Those are considered to be fully mixed and to have homogeneous physical and chem- ical conditions. This assumption is false at a small scale as the deeper zones of the Baltic Sea are strongly stratified and the nutrient concentration are different in the coastal and open sea areas. However, the goal of the model is not to give an exact value of the algal density for each point but to give an overview of the eutrophication situation over time and for the differ- ent part of the Baltic Sea.
The hydrodynamic is a crucial part off this mod- elling work, for the transport between sub-basins and vertical layers and for driving the ecosystem and sediment model, and has been thoroughly verified (Dargahi and Cvetkovic, 2014). However, the use of instantaneous and spatially discrete hydrodynamic data is a major source of uncertainty for the model as this data does not respect the continuity equation.
Temporally averaged data cannot be used as it is im- portant to take into consideration that there is water coming in and water leaving the block, that corre- spond to velocity values with opposite signs. Tempo- ral or spatial average would then give a unique value, therefore neglecting a part of the flow. A discussion about the possibilities of amelioration will be made in Section 3.4.
The other sets of data are also subject to different assumptions or uncertainties. The monitoring data is supposed to be representative of the DIN, DIP and algal concentrations for a sub-basin. There are also some uncertainties concerning the conversion of the total loads into DIN and DIP, and further work is needed in that area. Moreover, for simplicity reasons, the solar radiation is considered to be the same for all the Baltic Sea, which can be a source of uncertainties for the algal growth.
For the calculation of the density, the effect of the pressure has also been neglected. However, the value of interest is the gradient of density for the determi- nation of the pycnocline and as the pressure varies linearly with the depth; its effect can be neglected.
2.4.4. Model set-up
Initially, the preprocessing part has been used to aggregate the hydrodynamic data, as well as the loads, the solar radiation and the boundary conditions.
The time period of the simulation is between 04/04/2001 and 31/12/2009. Monitoring data could not be found for all the blocks for the time of the start of the simulation and, for some variables, there is no monitoring data. Therefore, a monitoring sta- tion was used to set the initial condition of the avail- able variables and 0 were used for the others. With those first initial conditions, a winding up simulation of 20 years, repeating the two first years ten times was carried out. This method has the advantage to make sure that the model converge towards a stable state and to use those values as initial conditions for all the variables.
Then, the values at the end of the 20 years simu- lation were used as initial conditions and the simula- tion was carried out for the entire time period. The results obtained for the simulation are presented in Section 3.
2.4.5. Calibration and validation
The comparison between monitoring and mod- elled data has been done using the mean absolute
Table 1: MASE before and after calibration.
Sub-basin Surface
DIN Deep DIN Surface
DIP Deep DIP Algal
density
Bothnian Bay (before calibration) 10.4 15 59.0 91.4
After calibration 6.0 7.0 0.7 0.8
Southern Baltic Proper (before
calibration) 13.7 14.1 5.7 4.6 8.2
After calibration 3.1 2.8 1.9 1.2 2.3
scaled error (MASE), which is the mean of the error scaled by the error of the naive forecast method. This method has the advantages of being scaled and there- fore can compare time series with different ranges and can be used for values close to zero (Hyndman and Koehler, 2006).
The calibration has been carried out by trying to reduce the MASE for the DIN and DIP concentra- tions between April 2001 and April 2003. When dif- ferent values were available for the same day, the aver- age has been taken. More emphasis has been put on the surface concentrations as when there are differ- ent values for the same day, their dispersion are less important than for the deeper concentrations. The chlorophyll monitoring data has not been used as it is scarce for most of the sub-basins. The sensitive parameters have first been identified by a local ap- proach, changing one parameter at the time. Then, a range of variation has been identified for the sen- sitive parameters, according to the literature and the parameters were changed together. The parameters are given in Appendix A.3. Table 1 gives the MASE for the northern and southern most blocks before and after calibration, for the calibration period (2001- 2002).
The validation consists of the comparison between monitoring and modelled data using the MASE, for the whole simulation period and is presented in Sec- tion 3.2.
3. RESULTS AND DISCUSSION In this part, the behaviour of the different vari- ables and their links with each other are explained.
The Figures 7.A and 7.B represent the different vari- ables for the Eastern Gotland Basin (Fig. 6).
The bloom of phytoplankton starts in the begin- ning of April, is maximal in the start of May, with an areal density of 165g .m−2and ends in the begin- ning of June (Fig. 7.a). This bloom is due to an in- crease of the solar radiation and of the water temper- ature, coupled with the availability of both DIN and DIP. As the concentration of DIN and DIP decrease and the self-shading factor increases, the growth rate decreases, to be inferior than the death rate and the phytoplankton density decreases during May. There
is a smaller bloom in July, due to an increase of the DIN and DIP concentration. The bloom in July is stopped by low DIN concentrations and by the rapid augmentation of the cyanobacteria density.
The cyanobacteria density starts increasing slowly at the end of the first phytoplankton bloom but the cyanobacteria bloom starts only in the middle of July and the maximal density of 38g .m−2is reached from August to the middle of September (Fig. 7.b). Un- like phytoplankton, cyanobacteria are able to fix at- mospheric nitrogen and thus do not depend on the DIN availability. This capacity, in addition to easily available DIP and increasing water temperature make it possible for them to outgrow the phytoplankton from the middle of July. Their density is then con- stant due to a decrease in their growth rate, as the DIP concentration is low in August.
The detritus correspond to the dead algae, and the nutrients (carbon, nitrogen and phosphorus) in detritus have similar trends. The surface detritus concentrations have a first peak in the middle of May, corresponding to the end of the phytoplankton bloom and therefore when the algal death decreases due to a decrease of the algal density (Fig. 7.g, 7.i and 7.k). The surface detritus concentrations then de- crease as they are either mineralized and transformed in DIN and DIP, or they sink until reaching the deeper part. The detritus increase again due to the smaller algal blooms. The deeper detritus undergo the same mineralization processes or sink until they reach the volatile sediment part. The deeper detritus comes from the sinking and the flow of surface detri- tus and its concentration is smoother as it is not di- rectly linked to the algae death (Fig. 7.h, 7.j and 7.l).
Their concentrations start increasing at the middle of May and its maximum is reached at the beginning of June, stay high until October and then decrease again until the next algal bloom.
The surface DIN is easily available at the begin- ning of the year (320,m g .m−3) and decreases with the first algae bloom, in the beginning of May (Fig.
7.c). As the DIN availability reduces, the algal density decreases and DIN concentration increases both through mineralization of the detritus and loads at the end of June (between 0 and 50m g .m−3).
This creates the second smaller algal bloom and the
