Comment on bg-2021-156 Anonymous Referee # 3 Referee comment on " Modeling cyanobacteria life cycle dynamics and historical nitrogen fixation in the Baltic Sea

This study incorporated a cyanobacterial life cycle model with phosphorus dependency, which improved the prediction of diazotrophic cyanobacterial blooms in the Baltic Sea. The research is quite interesting and challenging; however, I found the whole manuscript lacks a clear hypothesis, clear clarification of why phosphorus is important, and the interpretation of results is not deep enough. I could see that the authors were trying to explain the methodology as it is a complicated study, however I got lost easily as there is not a clear approach or conceptual diagram to lead the readers. I have also got a few major concerns as listed below.

The Baltic Sea is exposed to significant impacts from eutrophication (HELCOM 2010) because of the combination of large increases of nutrient supplies since World War II (Gustafsson et al. 2012), permanent stratification (e.g. Leppäranta 2009), and long water residence times (Meier 2007) that reduce the deep water ventilation and causes wide spread oxygen deficiency. Several large-scale geoengineering interventions have therefore been proposed as solutions to the eutrophication problems (e.g. Conley 2012). Nations around the Baltic Sea have decided on a Baltic Sea Action Plan to reduce external loads of nutrients to the area (HELCOM 2007). Despite reduced nutrient inputs (Gustafsson et al. 2012), there is still an increase in abundance of filamentous cyanobacteria during recent decades (Finni et al. 2001, Kahru and Elmgren 2014, Reusch et al. 2018. The growth of the filamentous cyanobacteria is sensitive to the availability of phosphate (Moisander et al. 2007, and phosphorus loads are therefore of extra importance to decrease to allowable levels suggested by the Baltic Marine Environment Protection Commission HELCOM (2018). With some cyanobacterial taxa being able to utilize both phosphate and organic phosphorus (Schoffeleen et al. 2018) complicates modeling efforts and knowledge on their relative importance is still limited.
The aims of the current study were to gain understanding in phosphorus dynamics in the Baltic proper as well as demonstrate the workings and boundaries of the CLC model in order to use it for continuous monitoring and estimates of nitrogen fixation for management purposes. This will be done by I), run sensitivity experiments addressing phosphorus limitation to determine the optimum settings for the Baltic proper in relation to cyanobacteria blooms, II) include the CLC model in a high-resolution 3D coupled physical and biogeochemical model of the Baltic Sea, and III), validate it to observations of cyanobacteria carbon biomass and estimated nitrogen fixation measurements based on previous in situ measurements.

Method
The Baltic Sea is a semi-enclosed estuary which has limited water exchange with the adjacent North Sea (Fig. 1). In order to study bloom formations of filamentous cyanobacteria, we included a modified version of the cyanobacteria life cycle (CLC) model in a high-resolution three dimensional (3D) coupled physical-biogeochemical model of the Baltic Sea (Meier et al. 2003, Eilola et al. 2009, Almroth-Rosell et al. 2011) spanning 1850-2008. The modified CLC model is described in detail below, together with modifications of the biogeochemical model setup (schematically shown in Fig 2).

Ocean circulation model
The RCO (Rossby Centre Ocean) model is a Bryan-Cox-Semtner primitive equation circulation model with a free surface (Killworth et al., 1991). Its open boundary conditions are implemented in the northern Kattegat, based on prescribed sea level elevation at the lateral boundary (Stevens, 1990). An Orlanski radiation condition (Orlanski1976) is used to address the case of outflow, and the temperature and salinity variables are nudged toward climatologically annual mean profiles to deal with inflows (Meier, 2003). A Hibler-type dynamic-thermodynamic sea ice model (Hibler, 1979) with elastic-viscousplastic rheology (Hunke and Dukowicz ,1997) and a two-equation turbulence closure scheme of the k-ε type with flux boundary conditions (Meier et al., 2001) is embedded into RCO. The deep-water mixing is assumed inversely proportional to the Brunt-Väisälä frequency, with the proportionality factor based on dissipation measurements in the Eastern Gotland Basin (Lass et al., 2003). RCO is here used with a horizontal resolution of 2 nautical miles (3.7 km) and 83 vertical levels, with a layer thickness of 3 m. RCO allows direct communication between bottom boxes of the step-like topography (Beckmann and Döscher, 1997). A flux-corrected, monotonicity-preserving transport (FCT) scheme is applied in RCO (Gerdes, 1991). RCO has no explicit horizontal diffusion. For further details of the model setup, the reader is referred to (Meier, 2003) and (Meier, 2007).
The model performance of temperature and salinity was evaluated in Meier et al. (2018) and conforms well to observations but shows a higher position of the halocline and slightly lower bottom water salinity. The modelled temperature shows good agreement with the observations but some deviations with higher temperatures are found in the upper part of the halocline.

Biogeochemical model
The biogeochemical model SCOBI (Swedish Coastal and Ocean Biogeochemical model) has been developed to study the nutrient cycling in the Baltic Sea (Marmefelt et al., 1999, Eilola et al., 2009, Almroth-Rosell et al., 2011, Almroth-Rosell et al., 2015. SCOBI handles biological and ecological processes in the sea as well as sediment nutrient dynamics and is in this study coupled to RCO (e.g. Eilola et al., 2012, Eilola et al., 2013, Eilola et al., 2014. Resuspension of organic matter is calculated, with the help of a simplified wave model, from the wave and current-induced shear stresses (Almroth-Rosell et al., 2011). SCOBI has a constant carbon (C) to chlorophyll (Chl) ratio C:Chl = 50 (mg C (mg Chl)-1), and the production of phytoplankton assimilates carbon (C), nitrogen (N) and phosphorus(P) according to the Redfield molar ratio (C:N:P = 106:16:1) (Eilola et al., 2009). The molar ratio of a complete oxidation of the remineralized nutrients is O2:C = 138. Dead 4 90 95 100 organic material, represented by separate variables for nitrogen and phosphorus accumulates in detritus in the water column and in the sediments. For further details of the "standard" SCOBI model, the reader is referred to Eilola et al. (2009) and Almroth-Rosell et al. (2011.

Cyanobacteria life cycle model
The CLC model we used is a modified version of the detailed life cycle model developed by Hense and Beckmann (2006) including internal nutrient quotas, and the simplified version by Hense and Beckmann (2010) where life cycle transitions only depend on external factors. The CLC model equations as well as variables and parameters can be found in the Supplementary material. Similar to Beckmann (2006, 2010), we pooled the three main important nitrogen fixing taxa N. spumigena, Aphanizomenon sp. and Dolichospermum spp. into one functional cyanobacteria group. We are well aware that there are differences among the species (e.g. with respect to salinity or temperature dependence) and thus we may not expect to be able to reproduce specific local patterns. Nevertheless, our model will be able to reproduce the main seasonal and spatial patterns of biomass and nitrogen fixation.
Growth and life cycle transitions in our CLC model depend on external factors, which is similar to Hense and Beckmann (2010), but we kept the sinking and rising stage separated. We thus distinguish between three life cycle stages: the growing and nitrogen fixing stage (vegetative cells with heterocysts, called HET), the resting stage (akinetes, called AKI) and a stage (called REC) where we combine the recruiting (cells with gas vesicles) and the growing, non-nitrogen fixing stage For the transition between AKI (AKIB and AKIW) and REC we prescribe a fixed germination window -from April 20 to the end of April, instead of using a dynamic germination window as proposed by Hense and Beckmann (2010). This is 5 115 120 125 130 because the computational costs in a 3D framework are too high. Shifting the germination window, however, has only a small impact on the timing of maximum cyanobacteria abundance in summer and the magnitude of nitrogen fixation. The decadal mean nitrogen fixation is lower when germination is earlier by about 4-7%, while the maximum annual difference found for the entire period 1850-2008 is lower by 14%.
Growth of HET and REC are inhibited under anoxic conditions. For potential growth and transition of AKI to REC we assume a salinity dependent window between 3 and 10 PSU which is in fair agreement with the optimum growth of N. The AKIB is assumed to be rapidly immobilized in the sediment under salinity outside of the defined range. This effect is simulated by very large burial of AKIB when salinity is lower than 3 PSU or higher than 10 PSU. The temperature limitation for the growth of HET and REC (Supplementary material, Table S.3, Eqs. 8 and 26) reaches 10% of its maximum (where maximum equals 1 and indicates least limited) at a temperature of 11°C. The growth is unlimited by temperature at 28°C after which a temperature increase means a decline in growth.

Model forcing
The historical simulation uses reconstructed atmospheric, hydrological and nutrient load forcing and daily sea levels at the lateral boundary for the period 1850-2008 as described in detail in Meier et al. (2018) and Gustafsson et al. (2012) and references therein. The used High Resolution Atmospheric Forcing Fields for the period 1850-2008 were reconstructed using atmospheric model data for 1958-2007 together with historical station data of daily sea-level pressure and monthly air temperature observations. For the calculation of monthly mean river flows five different historical data sets were merged.
The basin integrated reconstructed nutrient loads from land and atmosphere to the present model are the same as used and described by Gustafsson et al. 2012. Nutrient loads contain both organic and inorganic phosphorus and nitrogen, respectively. In the present SCOBI version, the nitrogen and phosphorus detritus were separated and thus used both organic phosphorus and nitrogen from the forcing. This is the only difference in forcing from the present SCOBI model compared to the model used by Meier et al. 2018, where detritus consisted of one pool limited by the Redfield ratio. Daily mean sea level elevations at the boundary in the Northern Kattegat were calculated from the reconstructed, meridional sea level pressure gradient across the North Sea. In case of inflow, temperature, salinity, nutrients and detritus values were nudged towards observed climatological seasonal mean profiles for 1980-2005 at the monitoring station Å17 in the southern Skagerrak.
Nutrient concentrations before 1900 were assumed to be only 85% of present-day concentrations. A linear decrease of nutrient concentrations from 1950 and back in time to 1900 was assumed.

Observations
The Swedish National Marine Monitoring Program includes monthly tube sampling of phytoplankton abundance (including filamentous cyanobacteria) and water collection for chemical and physical parameters (e.g. inorganic nutrients, oxygen, salinity, temperature). This data is hosted by the Swedish National Oceanographic Data Centre at the Swedish Hydrological and Meteorological Institute and is freely accessible at www.smhi.se. For this work we have also used data of oxygen and nutrients from The Baltic Environmental Database (BED) which includes post processed monitoring station data from a number of institutes around the Baltic Sea. The data is freely available at http://nest.su.se/bed. The cyanobacteria biovolume (mm3 l-1) is calculated based on cell numbers and size of filaments (Olenina et al. 2006) and further to carbon concentrations (referred to as cyanobacteria biomass) based on Menden-Deuer and Lessard (2000). Concentrations of inorganic nutrients and oxygen were extracted from the database for station BY15 in the Eastern Gotland Basin, and cyanobacteria biomass for four stations in the Baltic proper for 1999-2008 (Fig. 1). The cyanobacteria biovolume was used to estimate nitrogen fixation rates (mmol N m-2 d-1) based on empirical species-specific measurements (Klawonn et al. 2016) according to Olofsson et al. (2021). The estimated nitrogen fixation rates were also calculated to annual nitrogen loads for the Baltic proper during 1999-2008, assuming a size of 200 000 km2, and provided as kton N yr-1.

Phosphorus dependence
In the original model by Hense and Beckmann (2006) that includes also the internal energy and nitrogen, the seasonal changes in cyanobacteria biomass are adequately modelled without taking phosphate into account. The rapid decrease of HET in autumn is then a result of an internal energy crisis caused by the high energy demand of nitrogen fixation together with decreasing temperatures and light. This is true also for the simplified model of Hense and Beckmann (2010) where the growth rate of HET is strongly limited by temperature. However, in the Baltic Sea, the phosphorus concentrations may limit the cyanobacteria biomass (Degerholm et al., 2006). We have therefore performed a series of sensitivity runs to evaluate the role of phosphorus uptake. In the experimental runs, we distinguish between uptake of inorganic and organic phosphorus, since both types are utilized by cyanobacteria (Schoffelen et al. 2018). A preferential uptake of dissolved inorganic phosphorus is, however, assumed in the model. REC and HET are assumed to assimilate C, N and P, and produce detritus, according to Redfield molar ratios.
In the first experiment (noP), we excluded phosphorus dependence from cyanobacteria in line with Hense and Beckmann (2010). In this case, the cyanobacteria can grow completely independent of phosphate availability in ambient water. In the next experiment (sPlim), we included strong limitations from both phosphate and organic phosphorus. In this case, the half saturation constants are large and the cyanobacteria growth depends strongly on the availability of phosphate.
In the third experiment, we assigned a very small value (10 -6 ) to the half saturation constants of both the phosphate and the organic phosphorus limitation terms, effectively removing the phosphorus limitation of cyanobacteria (wPlim). As long as phosphate exists in small amounts in ambient water, the growth is maintained independent of the concentration. However, in the absence of phosphate the growth is terminated.
In the fourth and final experiment, we kept the limitation by inorganic phosphate but removed the ability to utilize organic phosphorus in cyanobacteria (noOP). As can be deduced from Eq. 1 in Supplementary material, Table. S.3, the growth, in this case, gets no additional reinforcement from organic phosphorus.
The differences in parameter values between the phosphorus sensitivity runs are found in Table S.5 in the Supplementary material.

Results and discussion
We have used a novel combination of a 3D-model and a cyanobacteria life cycle model (CLC) for the Baltic Sea. In order to determine the optimum phosphorus settings for estimating cyanobacteria biomass and nitrogen fixation rates we performed four phosphorus limitation experiments where after the optimum were used for the estimates and validated to in situ observations.

Phosphorus sensitivity of cyanobacteria
The simulated biomass was generally larger than observations in all four phosphorus limitation experiments (Fig. 3). The experiment that generated the largest biomass was by far noP which completely excludes the impact of phosphorus in 8 195 200 205 210 ambient water. Since the cyanobacteria were not dependent on phosphate or nitrate, they grew extensively even in the first part of the 20th century (Fig. 4) when observations indicate that cyanobacteria blooms were seldom observed (Finni et al, 2001). Up until about 1980, the noOPlim (no additional contribution from organic phosphorus) generally generated the lowest annual mean. After noP, wPlim generated the highest annual mean. This is to be expected as the growth rate, given the same amount of organic and inorganic phosphate, in this case is largest (see Eq. 1 and Eq. 4 in Table S.3). Fig. 4 further displays a decline in cyanobacteria biomass from the mid 20th century to the 1980s for experiment noP. The reason for this decline is most likely the sharp increase in nutrient loads (Gustafsson et al., 2012) generating a competitive advantage of faster growing diatoms and flagellates leaving less DIN for the non-nitrogen fixing RECs . This is indicated also by the increase in diatoms and other phytoplankton biomass accompanying the cyanobacteria decline (Fig. 4, lower panel).
The experiments sPlim, wPlim, and noOP generated results closer to observations as compared to noP (Fig. 3). The lowest biomass was generated through wPlim for all stations except BY31. It is notable that wPlim also generated the best bloom timing as compared to observations. This experiment allows cyanobacteria to grow quickly even at low phosphate concentrations as long as the temperature is above approximately 11 degrees C (see Eq. 8 in Table S.3). Thereby the temperature sets the threshold of when the bloom will be initiated, and the phosphate dependence decides the end of it. The noP run generates a start of bloom that is close to wPlim but a later termination. The best seasonal timing validates our use of wPlim in our nitrogen fixation runs.

Including the cyanobacteria life cycle model
To estimate cyanobacteria biomass and nitrogen fixation rates we used the combined 3D and CLC model for the Baltic proper region using the phosphorus limitation setting wPlim based on the limitation experiments. Both, the model results and observations, captured an increase in filamentous cyanobacteria biomass in May-June and with a maximum abundance in July-August (Fig. 3). This is a typical seasonal cycle of cyanobacteria in the Baltic proper (Olofsson et al. 2020), and is an important improvement attained by using this model combination as compared to previous results (Hieronymus et al. 2018).
This improved seasonality is due to the inclusion of the cyanobacteria life cycle model (Hense and Beckmann, 2010). In earlier studies using models, the bloom of filamentous cyanobacteria was initiated too late in the season, resulting in a very low nitrogen fixation due to the temperature dependence in the model and decreasing water temperatures during fall (Hieronymus et al. 2018). By obtaining a bloom more constrained to the summer months, a larger nitrogen fixation due to higher temperatures was observed in the present study. The updated nitrogen fixation rates were also in the same range as estimates based on measurements for the same stations during the years 1999-2008, both in magnitude and timing (Fig. 5).
For nitrogen fixation, there was a slight difference where the model displayed a prolonged peak period in July -August while the observations showed a peak more contained to July. The strong coherence between model results and observed nitrogen fixation is somewhat surprising given the larger cyanobacteria biomass displayed by all model experiments compared to observations (Fig. 3).
We estimated the internal nitrogen load via nitrogen fixation to the Baltic proper based on monitoring and in situ measurements to a mean of 399 kton per year for 1999-2008, but with a large variation among years (SD ± 104). This is slightly below the external load from river runoff and atmospheric deposition of 430 kton per year (± 54), provided by HELCOM (2018). For the model combination we used herein, we got an estimated mean nitrogen load of 362 kton per year for experiment wPlim over the same years for the Baltic proper (calculated over an area of 216,600 km2). The estimated annual nitrogen load via nitrogen fixation to the Baltic proper has not changed over recent years ( -2017( in Olofsson et al. 2021, and is in the range of other studies for the Baltic proper (310 kton in Rolff et al. 2007;370 kton in Wasmund et al. 2001; 396 kton in Svedén et al. 2016) but below the estimated load of 613 kton in Wasmund et al. (2005), 514 kton in Gustafsson et al. (2013), and 511 kton in Schneider et al. (2009). For the Bothnian Sea, however, the estimated nitrogen load via nitrogen fixation has more than doubled from 1999-2004 until 2012-2017, with an increase in filamentous cyanobacteria along with decreased salinities (Olofsson et al. 2020).

Nutrients and oxygen
Diazotrophic cyanobacteria increase bioavailable nitrogen through release of ammonium from its newly fixed nitrogen (Ploug 2010;. They also impact surrounding organisms by competing for phosphate. We therefore demonstrate the mean seasonal cycle of phosphate and nitrate mean concentrations provided by the model for the periods 1999-2008 and 1960-1980 (Fig. 6). In the earlier period, wPlim consistently generates the highest biomass and noOP the lowest as expected from the lower growth rate obtained in this case (cf Eq 1., Table S.3 ). During the later period, dissolved inorganic nitrogen (DIN; nitrate and ammonium) is completely depleted after the spring bloom providing little opportunity for other phytoplankton than cyanobacteria to grow (Fig. 6). However, during the earlier part of the run, DIN is available even during summer allowing for higher biomass of the surrounding diatoms and other phytoplankton (Fig. 4)  The mean vertical profiles of phosphate, DIN and oxygen at stations BY5 and BY15 for experiment wPlim show an overall good representation by the model (Fig. 8). Below the mixed layer, the DIN concentrations are high compared to observations and the phosphate at BY5, a bit too low. The low surface DIN is a reflection of low nitrate concentrations as compared to observations which is also reported by Meier et al. (2012) and Saraiva et al. (2018). The low surface DIN is also seen in Fig.   7 where the noP experiment gives rise to higher DIN concentrations as nitrogen fixation due to strong cyanobacteria blooms in this case adds more DIN to the water column. Despite the shortcomings, the trends for both nutrients and oxygen are well captured by the model.

Summary and conclusions
Through a series of sensitivity experiments, we have shown that the inclusion of a weak phosphate limitation is essential for the CLC model in the Baltic Sea. Excluding this dependence generates too high concentrations of cyanobacteria, especially in the first part of the 20th century when cyanobacteria blooms were rarely observed. The large primary production in this case was also reflected in too high phosphate concentrations as eutrophication induced anoxia which gave rise to sedimentary phosphate release.
By including the CLC model into the 3D model for the Baltic proper we demonstrate a clear improvement in seasonality of blooms as compared to previous studies (Hieronymus et al. 2018 differences between the dominating taxa (Klawonn et al. 2016). Aphanizomenon sp. for example can perform equally high nitrogen fixation rates in 10 o C during spring as during the summer (Svedén et al. 2015), and they are responsible for the highest total nitrogen fixation in the region due to its long growth season (Klawonn et al. 2016). Aphanizomenon sp. may also use different sources of phosphorus, which may further discriminate the growth niches by the filamentous cyanobacterial species (Shoffelen et al. 2018). Phosphorus cycling is a complex topic, which also needs further studies in natural ecosystems, as high turnover rates of phosphorus of only about 2 h are hard to trace (Nausch et al. 2018). To include more species in the model might be of extra importance as climate change scenarios can change the community composition in the future (Wulff et al. 2018;Olofsson et al. 2020).
In this work, we have used a CLC model that includes benthic and pelagic akinetes from which the summer blooms originate. Research has shown that the life cycles of the different major bloom forming taxa are complex and there is no single answer on how they start growing after winter (Munkes et al. 2021). Experiments have suggested that all taxa form akinetes to some extent but the summer bloom of Nodularia spumigena and Aphanizomenon sp. originates mainly from small, overwintering water column populations while Dolichospermum spp. seems to originate from both akinetes and pelagic filaments (Wasmund et al. 2017, Suikkanen et al., 2010. The large improvement in seasonality when the lifecycle of cyanobacteria is modelled, as opposed to earlier modelling attempts that include only small winter populations, does however indicate that the separation into different lifecycle stages is of key importance for capturing the start and end of bloom.
Capturing the seasonality of cyanobacteria blooms is of great importance due to their impact on water quality as well as for obtaining better estimates of nitrogen fixation that contributes to eutrophication. This work constitutes a step forward for the modelling of cyanobacteria blooms in the Baltic Sea. The inclusion of CLC can with some further development be used to merge observations and modeling for obtaining better prognostic estimates of cyanobacteria blooms, which can be used for management purposes.

Code availability
The model code of the ocean model used for the simulations is publicly available from the Swedish Meteorological and Hydrological Institute, Norrköping, Sweden (https://www.smhi.se, E-mail: smhi@smhi.se).  . The inorganic nutrients nitrate, ammonia and phosphate are represented by NO3, NH4 and PO4, respectively. The phytoplankton groups A1 and A2 represent characteristics of diatoms and the flagellates and others. The bulk zooplankton ZOO grazes on phytoplankton A1 and A2 while the parameterized predation closes the system of equations. Nitrogen and phosphorus detritus are described by ND and PD, respectively. Oxygen dynamics are included and hydrogen sulfide concentrations are represented by "negative oxygen equivalents (1 ml H2S l-1 = -2 ml O2 l-1). The process descriptions of oxygen and hydrogen sulfide are simplified for clarity.      1960-1979 (upper) and 1999-2008 (lower). The data-points have been smoothed using a 1 month moving average.