Drivers of phytoplankton responses to summer storms in a stratified lake: a modelling study

Extreme wind events affect lake phytoplankton amongst others by deepening the mixed layer and increasing internal nutrient loading. Both increases and decreases of phytoplankton biomass after storms have been observed, but the precise mechanisms driving these responses remain poorly understood or quantified. In this study, we coupled a one-dimensional physical model to a biogeochemical model to investigate the factors regulating short-term phytoplankton responses to summer storms, now and under expected warmer future conditions. We simulated physical, chemical and biological dynamics in Lake Erken, Sweden, and found that wind storms could increase or decrease the phytoplankton 4 concentration one week after the storm, depending on antecedent lake physical and chemical conditions. Storms had little effect on phytoplankton biomass if the mixed layer was deep prior to storm exposure. Higher incoming shortwave radiation and hypolimnetic nutrient concentration boosted growth, whereas higher surface water temperatures decreased phytoplankton concentration after storms. Medium-intensity wind speeds resulted in more phytoplankton biomass after storms than high-intensity wind. Simulations under a future climate scenario did not show marked differences in the way wind affects phytoplankton growth following storms. Our study shows that storm impacts on lake phytoplankton are complex and likely to vary as a function of local environmental conditions.


Introduction
High wind speeds during storms reshape the lake physical and chemical environment in ways that alter phytoplankton biomass and growth. At the start of a chain of processes, wind stress at the lake surface induces internal mixing and this deepens the thermocline . These processes affect the vertical distributions of oxygen and nutrients, and in turn phytoplankton growth and vertical distribution. Upwelling of nutrientrich water occurs during mixing events and can alleviate nutrient limitation, potentially causing phytoplankton blooms (Soranno et al. 1997;Kasprzak et al. 2017;Whitt et al. 2019).
At the same time, surface temperature tends to decrease during storms (Kuha et al. 2016;Mesman et al. 2020), potentially reducing light-saturated phytoplankton growth rates (Trombetta et al. 2019). Additionally, a deeper mixed layer can increase light limitation for growth (Diehl et al. 2002) and deepening dilutes concentrations by mixing phytoplankton over a larger volume of water (Kuha et al. 2016). Sediment resuspension due to shear stress in shallower parts of the lake may simultaneously release nutrients and limit light availability (Ji et al. 2018). As such, there are conflicting effects of storms on nutrient and light availability (Stockwell et al. 2020), and the net effect of a storm on phytoplankton concentrations may depend on lake physiography and lake state prior to the event.
Changes in storm characteristics and lake thermal structure will affect phytoplankton responses to storms. As a result of climate change, extreme wind events will likely shift in frequency and intensity, with different parts of the globe experiencing increases or decreases (Mölter et al. 2016;Sainsbury et al. 2018). Concurrent with these changes in meteorological conditions, as surface water temperatures increase and stratification strengthens (O'Reilly et al. 2015;Pilla et al. 2020), more energy is needed to mix the water column (Schmidt 1928), so that in a warmer climate a wind event of a given magnitude and duration may cause less mixing. The depth of the mixed layer, at the time when a storm hits, also determines the degree of lake mixing; a deeper pre-storm mixed layer reduces entrainment of hypolimnetic water by surface waves (Imboden and Wüest 1995). Long-term climate effects on mixed layer depth (MLD) are still ambiguous; both shoaling and deepening mixed layers have been observed, in addition to non-significant trends (e.g. Kraemer et al. 2015;Pilla et al. 2020) and it is likely that local trends in transparency or wind speed are at least equally important as trends in warming to determine changes in MLD (Persson and Jones 2008;Woolway et al. 2019). Lastly, stratification is expected to occur earlier in the year as the climate warms (Woolway et al. 2021), and this will lead to a longer separation of epilimnion and hypolimnion. Therefore, there will be a greater build-up of nutrients in the hypolimnion during the stratified period Nowlin et al. 2005), which could be entrained into the mixed layer during a storm and become available for phytoplankton growth.
Understanding mechanistically how complex lake ecosystems are reshaped by storm events under present and projected future conditions is a challenging task. Process-based models have the advantage of allowing a quantitative comparison between different scenarios and identifying clear causal pathways even in complex systems and in conditions yet to be observed. Also, experiments that incorporate deep-water mixing and different scenarios regarding stratification, nutrient concentrations and climate warming are very demanding to set up (but see Giling et al. 2017), whereas models can relatively easily explore such a wide range of scenarios. Another issue involving the study of extreme events, is that such events are rare by definition and hard to predict. Moreover, storms act on short timescales (hours to days), while lake monitoring programs often include only weekly or monthly samples.
Therefore, biogeochemical data describing responses to storm events are scarce. In the present study, we use the General Ocean Turbulence Model (GOTM) model coupled to a biogeochemical model to study the impact of storms on phytoplankton dynamics. In Mesman et al. (2020), one-dimensional process-based models (including GOTM) were shown to simulate physical effects of storms with reasonable accuracy. This means that we can have some confidence that processes related to the transport and mixing of biogeochemical particles and solutes during storms are accurately simulated. Our approach allows us to draw conclusions about potential regulating factors for the response of phytoplankton to storms, and how climate warming may affect this response.
The main process under investigation here is the effect of storms on the redistribution of biogeochemical compounds between the epi-and hypolimnion. We use Lake Erken, a Swedish mesotrophic dimictic lake, as a case study, because of the available long-term time series of physical and biological variables that we used to calibrate the model. However, the findings in this study increase our understanding of the processes regulating phytoplankton responses to storms across stratifying lakes in general. Phytoplankton communities in these lakes are shaped by their need for both nutrients and light, which show opposing gradients in availability. Shallow polymictic lakes are likely to react differently to storms, with more emphasis on sediment resuspension and uprooting of macrophytes (Ji et al. 2018), and our model results are not applicable there. We assess scenarios covering a broad range of atmospheric and lake conditions, that reflect conditions present in many temperate, stratifying lakes.
Here, we investigate 1) how storm intensity, thermal structure, light availability, and nutrient availability control phytoplankton response to storms during summer stratification, and 2) how climate warming may influence the response of phytoplankton to storms. To answer these questions, we performed two numerical experiments. In the first experiment, we repeatedly simulate a storm event while changing storm intensity, incoming shortwave radiation, and pre-event mixed layer depth, surface water temperature, and hypolimnetic nutrient concentration in a full factorial design. In the second numerical experiment, we compare the response of phytoplankton to wind perturbations between present-day and future-climate air temperatures, at different times of the year and at different storm intensities. These simulations help us to disentangle and better understand the dynamic response of primary producers to storms in a changing world.

Site description
Lake Erken is a mesotrophic lake in Sweden (59°50′37′′ N, 18°35′38′′ E), with a mean depth of 9 m and a maximum depth of 21 m. The lake has a surface area of 24 km 2 and its retention time is 7 years (Blenckner et al. 2002). Lake Erken is dimictic, meaning that it experiences both winter ice cover and summer stratification, although short-term partial or complete mixing events are possible in summer in response to wind-induced mixing (Yang et al. 2016a). During summer stratification, both nitrogen (N) and phosphorus (P) can limit phytoplankton growth in the lake (Vrede et al. 1999), whereas during deep mixing or fully mixed conditions, light availability is the main limiting factor (Yang et al. 2016a). During summer, nutrient concentrations build up in the hypolimnion, and these nutrients are circulated through the complete water column after the autumn turnover ).
In most years, Lake Erken experiences a distinct spring bloom followed by a clear water phase, and then a second phytoplankton biomass peak in summer/autumn (Yang et al. 2016b), a pattern followed by many monomictic and dimictic lakes across the globe (e.g. Sommer et al. 2012). The spring bloom in Lake Erken is dominated by diatoms such as Aulacoseira spp., Stephanodiscus spp., and Asterionella formosa (Weyhenmeyer et al. 1999;Yang et al. 2016b). In summer, a major bloom-forming species is the cyanobacterium Gloeotrichia echinulata (Karlsson-Elfgren et al. 2003;Yang et al. 2016b).

Data collection
In this study, we used meteorological, water temperature, and biogeochemical data for the period 1999 to 2020. Meteorological data (wind speed, air temperature, air pressure, relative humidity, shortwave radiation, cloud cover, and precipitation) were collected using a weather station on a small island in the lake at hourly frequency. Moras et al. (2019) replaced missing meteorological data from nearby stations, selected by artificial neural network analysis, and we continued to use this dataset, supplemented by data until the end of 2020.
Discharge into Lake Erken was calculated by the HYPE model (Lindström et al. 2010), and validated using measured data from the main tributary of the lake. In this main tributary, discharge and temperature data were automatically monitored and summarised to daily values, while phosphate, total phosphorus, nitrate, and particulate organic matter concentrations in the inflow were measured once or twice per month. These data were collected only from 2004 onwards, and the first 5 years of the recorded nutrient loadings were recycled for 1999-2003, which was the spin-up period of the model (see below).
In the lake, hourly water temperature data were collected during the ice-free season with a thermocouple chain above the deepest point of the lake, every 0.5 m down to 15 m depth.
Starting in 2017, these data were collected year-round. Schmidt stability (Schmidt 1928;Idso 1973) and mixed layer depth were calculated from the water temperature profiles.
Schmidt stability was calculated using the "rLakeAnalyzer" R package (Winslow et al. 2019).
The mixed layer depth was defined as the depth where water density had increased by 0.15 kg/m 3 relative to the uppermost measurement (similar method as Wilson et al. 2020).
Water samples to determine nutrient concentrations (phosphate, total phosphorus, nitrate, and ammonium) were collected every two weeks during the ice-free season. During stratification, separate integrated samples of the epilimnion and hypolimnion were taken. In winter, if the ice was accessible, a single integrated nutrient sample was taken through a hole in the ice approximately every month. Nutrients were analysed using standard laboratory techniques, described in Ahlgren and Ahlgren (1976) and Goedkoop and Sonesten (1995). Chlorophyll-a data were collected at the same time and depth resolution as the nutrient data. Material concentrated by filtration on glass fibre filters were analysed using spectrophotometry as described in Ahlgren and Ahlgren (1976).

Model description and setup
The General Ocean Turbulence Model (GOTM) is a one-dimensional (1D) k-epsilon model that simulates vertical thermal and turbulence dynamics in freshwater and marine water bodies (Umlauf et al. 2005). GOTM is interfaced to the Framework for Aquatic Biogeochemical Models (FABM), which allows coupling of a physical model with a biogeochemical model (Table 1, Bruggeman and Bolding 2014). At every simulation time step in this coupled setup, the biogeochemical equations are applied to each layer in GOTM, including surface and sediment exchange, and GOTM regulates the transport of biogeochemical substances between the layers. Using FABM, GOTM was coupled to a modified version of the SELMA model, which itself is a modular version of the ERGOM model (Table 1, Neumann et al. 2002). This new version was named "Selmaprotbas", because apart from several code improvements, we implemented several features from the PROTBAS model described by Markensten and Pierson (2007).
The Selmaprotbas model describes oxygen, detritus, nitrogen, phosphorus, phytoplankton and zooplankton dynamics. Processes described in the model include (de-)nitrification, sediment resuspension, sediment solute release, mineralisation of detritus, phytoplankton growth regulated by nutrients and light, and grazing by zooplankton (Neumann et al. 2002).
Chlorophyll-a content and biomass are linked through a fixed chlorophyll-to-carbon ratio.
We modified the SELMA model code by 1) expressing biomass in carbon instead of nitrogen; 2) adding a silica cycle; 3) adding an option to use the phytoplankton light limitation and temperature growth dependence function described in Reynolds et al. (2001), 4) relating chlorophyll-a concentration directly to the carbon biomass of phytoplankton; 5) adding the possibility for buoyancy regulation of phytoplankton; and 6) allowing varying nutrient ratios over time in detritus and sediment. Advantages of these changes include a more comparable set-up to other biogeochemical models and a more complete description of potentially relevant processes. A more detailed common-language description of the model has been supplied in Supplement A and the model code is publicly available (see Software Availability). The meteorological conditions and inflow data collected at Lake Erken were used as inputs for the model, and the GOTM model was run with an integration time step of 1 hour and 0.5 m thick layers. A fourth order Runge-Kutta time integration scheme was chosen for the Selmaprotbas model. The Selmaprotbas model was run with two phytoplankton groups: diatoms and cyanobacteria, both of which had growth regulated by light, temperature, and the concentrations of phosphorus and nitrogen. The diatom group was calibrated specifically for the spring period, had high sinking rates, and was also regulated by silica; cyanobacteria could fix nitrogen and regulated buoyancy based on light availability, with the same settings as Anabaena (now Dolichospermum) described by Reynolds et al. (2001). For the files used to run the model, see Software Availability.

Calibration
After 5 years of spin-up, 13 years (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) were used for calibration, which was done using the ParSAC software (Table 1, Bruggeman and Bolding 2020), employing a differential evolution method to optimise the maximum likelihood objective function of RMSE between observations and simulations. A single set of parameter values was retrieved from the calibration. The calibration was split into two steps. First, the water temperature data was optimised using 10,000 iterations, by calibrating five parameters based on a previous study ): minimum turbulent kinetic energy, the extinction coefficient of visible light, and scaling factors for heat fluxes, wind speed, and incoming shortwave radiation. In the second step, 400,000 model iterations were done, varying 50 parameters that were determined in a sensitivity analysis (see next section). The objective function compared simulated in-lake concentrations of ammonium, nitrate, phosphate, total phosphorus, chlorophyll-a, water temperature, and oxygen with measured lake data. Physical parameters from the first step were also included in the second step of the calibration, but their ranges were constrained to +/-10% of the value obtained in the first calibration step.
Since the nutrient and chlorophyll-a samples were based on integrated samples for the epilimnion and hypolimnion, for the calibration we assumed these samples to be representative of 3 m and 15 m depth, respectively.
The ParSAC calibration attaches equal weight to each observation, so to avoid attaching too much value to the water temperature and oxygen measurements (which were collected at higher frequency) in the second step, we excluded temperature and oxygen measurements that were not collected on the same day as the nutrients, and we reduced the vertical resolution from 0.5 to 1.0 m. The full time series were used to assess goodness-of-fit.
The results of the calibration can be found in Supplement B.

Sensitivity analysis
To determine what parameters to include in the calibration of the Selmaprotbas model (biogeochemistry), we performed a sensitivity analysis using the ParSAC software We followed a density-based delta-sensitivity method (Borgonovo 2007), which has been described by Andersen et al. (2021)  We then ran the model for all parameter sets in the Latin hypercube and used the delta-sensitivity analysis as described by Borgonovo (2007) on the results to do the sensitivity analysis. In this method, the global importance of a parameter is calculated based on its effect on the entire output distribution, which can be calculated even when parameters are correlated (Borgonovo 2007). A 95-percent confidence interval around the sensitivity values was obtained by 100 resamples using a bootstrapping approach (Plischke et al. 2013). To distinguish sensitive from insensitive parameters, we introduced a dummy parameter in the sensitivity analysis (Andersen et al. 2021). If a parameter's sensitivity value fell within the 95-percent confidence interval of the dummy for all variables (ammonium, nitrate, phosphate, chlorophyll-a, oxygen), that parameter was excluded from the calibration. For parameters that were repeated for the two phytoplankton groups, the parameter was excluded only if it fell in the dummy confidence interval for both groups. If not, the parameter was retained for both phytoplankton groups.
The parameters that were excluded during the sensitivity analysis can be found in Supplement B.

Validation
We termed the "long-term simulation".
All calculations and data handling were done in the R software, version 4.0.1 (R Core Team 2020).

Numerical experiment 1: varying initial conditions before a storm
The aim of our first numerical experiment was to investigate which variables control the response of phytoplankton concentrations to storms. In order to achieve this aim, we induced a one-day storm event for different values of several meteorological and pre-event lake variables. Five variables were chosen that were expected to impact phytoplankton response to storms: -Storm intensity (i.e. wind speed during the event); this represents the magnitude of the disturbance induced by the storm -Mixed layer depth (MLD); MLD prior to the storm is a measure of the vulnerability of the thermal structure to mixing and controls the volume of the epilimnion -Shortwave radiation; incoming solar radiation regulates the availability of light -Surface water temperature; increases the strength of stratification, and therefore the resistance to mixing -Hypolimnetic nutrients; these regulate the potential for nutrient upwelling Ten levels of each variable were taken into consideration in a full factorial design, therefore totalling 10 5 = 100,000 simulations. The storm perturbation had a duration of 24 hours and was initiated 24 hours after initialisation of the model run ( Figure 1a).
We focused on the month of July to generate the weather conditions for this numerical experiment. This is because Lake Erken was always stratified in this month and stratification would have existed for a long enough time to allow for the build-up of hypolimnetic nutrients. In order to have a representative period with natural weather variations, we have selected the year in our dataset with the most generic weather conditions during this month. To that purpose, we calculated the mean and standard deviation for the meteorological driving variables measured in July (wind, pressure, temperature, relative humidity, shortwave radiation, and cloud cover) for each year separately and for the full period . For each year, we then calculated the root mean squared relative error (e.g. Despotovic et al. 2016) between the year and the full period values for both mean and standard deviation. July 2006 had the lowest error value and therefore most closely matched the long-term mean and variance, so these weather conditions were used as baseline in the first experiment ( Figure 1, panels a-d).
This numerical experiment required calculation of generic in-lake conditions to be used as initial conditions, for all model variables. We used the long-term simulation during July to determine these, because we decided that this was more representative of average summer conditions in Lake Erken rather than using the initial profiles for July 2006. Due to the scenarios with different MLD and the importance of MLD for vertical profiles of temperature and solutes, we calculated average profiles for the epilimnion and hypolimnion separately.
Both profiles were defined as ten equidistant values, interpolated from surface to the MLD for the epilimnion, and from the MLD to the maximum depth for the hypolimnion. As an example, when the MLD was 8 m, the epilimnetic profile would consist of ten values, a linear interpolation of the simulated values over 0 -8 m, and when the MLD was 4 m, the epilimnetic profile would still have ten values, but interpolated over 0 -4 m depth. Dates with a density difference between top and bottom of less than 0.15 kg/m 3 were excluded, as well as dates with an MLD shallower than 2 m or deeper than 15 m. The average of all calculated epilimnetic and hypolimnetic profiles during July were taken as initial conditions for our simulations, for each model variable. Because these initial conditions were defined relative to MLD, profiles for any value of MLD could be generated (see dashed lines in Figure   1, panels e-g). The ten different levels of wind speed, MLD, incoming solar radiation, surface temperature, and hypolimnetic nutrient concentration used in the simulations were also determined from the long-term simulation (1999-2020) using data from July. The 5th and 95th percentile of the daily averages were computed, and a linear sequence of ten values between these percentiles was used for the experiment (Figure 1). The only exception was wind speed.
Since the main interest of the present study was high wind speed, the lowest wind speed considered was the median, and the highest wind speed was set to 3 m/s above the maximum recorded daily average wind speed, to anticipate the possibility of increased storm intensity in the future due to climate change (Mölter et al. 2016). Each combination of the independently-changed variables was a separate simulation.
As output of each model run, concentrations of chlorophyll-a, nitrate, and phosphate were volume-averaged over the euphotic zone for the first week after the end of the storm. The depth of the euphotic zone (i.e. the depth where 1% of the light remains) was kept constant and was calculated as -ln(0.01) times the calibrated value of the e-folding depth of visible light ("g2" parameter) in the GOTM model, which gave a euphotic zone depth of 8.0 m. We decided to average over the euphotic zone, because phytoplankton below this zone are unlikely to contribute to primary production and will likely degrade over time, yet they are still present as biomass in the model. Averaging over the full water column would therefore lead to underestimation of the effects of storms on production and chlorophyll-a due to mixed layer deepening. In addition to the average concentrations, average Schmidt stability was also calculated for the first week after the storm.
As the last step of this experiment, we fitted the volume-averaged chlorophyll-a concentrations in the experiment with a random forest model to discern what variables affected the result most, using permutation variable importance. The random forest model contained 1000 trees, and all data from the experiment were used -so no holdout -, as the aim was to find out what variables were most important, not prediction. The fitting of the random forest model and the calculation of the permutation importance were done using the "ranger" R package (Wright and Ziegler 2017).   in other variables that may be specific to Lake Erken were not included. Atmospheric warming not only influences surface temperature and strength of stratification, but can also change the mixed layer depth and lead to an earlier onset of stratification and thus to different vertical profiles of nutrients and oxygen. As such, increased air temperature influences multiple potentially important lake variables that affect the relationship between storms and phytoplankton, which were investigated in isolation in the first numerical experiment.

Numerical experiment 2: inducing storms in long-term scenarios
In each simulation, the volume-averaged chlorophyll-a concentration between 0 and 8.0 m depth (the euphotic zone) was calculated. The maximum difference in concentration between the control and storm scenario, within one week after the storm, was compared between present-day and future climate.

Calibration and validation
During the calibration period, the model simulated water temperature with an RMSE of 0.9 °C. Seasonal cycles in oxygen, nitrate, and phosphate, were also reproduced, as indicated by NSE values above 0 ( Table 3). The main cause for the low fit statistics of chlorophyll-a was the substantial underestimation of the spring chlorophyll-a peak (Supplement E). In the validation period, the model fit worsened slightly for phosphate, nitrate and chlorophyll-a, as indicated by the NSE values. Inspection of the time series (Supplement E) confirmed that seasonal cycles in water temperature were well simulated. Deep-water oxygen concentrations were also simulated accurately, except that under ice, oxygen depletion was underestimated by the model, and in 2014-2016 deep-water anoxia events were missed.
Chlorophyll-a concentrations showed distinct spring and summer peaks, but in almost all years, the spring peak concentrations were underestimated, and sometimes simulated too late. Summer chlorophyll-a levels tended to be close to observed levels, with the exception of some summer blooms. Epilimnetic concentrations of both nitrate and phosphate were simulated to be low in summer, typical of values measured in the lake. However, the increase in nitrate concentrations in autumn was reproduced too early, and winter concentrations of phosphate tended to be underestimated. Chlorophyll-a (µg/l) 6.163 6.820 3. Of the other variables included in the experiment, incoming shortwave radiation had the strongest influence (Table 4). At low incoming radiation (< 150 W/m 2 averaged), the effect of storms on phytoplankton was largely negative, although at moderate wind speeds some increases could still be seen (Figure 2b). When incoming radiation increased above 200 W/m 2 , storms had an overall positive effect on phytoplankton concentration, moderate wind speeds more so than high wind speeds. Hypolimnetic nutrients were less influential than light or surface water temperature (Table 4), but higher concentrations resulted in increased phytoplankton biomass following storms at mixed layers deeper than about 4 m when incoming radiation was high (Figure 2b). At the shallowest mixed layer (2 m), wind speed had a negative effect on chlorophyll-a when hypolimnetic nutrients were high. This was caused by the setup of the experiment (Figure 1), as even under the lowest wind speed a large amount of the hypolimnetic nutrients entered the epilimnion, and high wind speed therefore mostly reduced growth due to mixed layer deepening. Surface water temperature before the storm (and therefore Schmidt stability) also influenced phytoplankton response to storms; a higher initial surface temperature caused a decrease of phytoplankton after storms (Supplement F-2).
Apart from changes in chlorophyll-a, we also looked at changes in Schmidt stability, and nutrient concentrations. The strongest decrease in Schmidt stability was diagnosed at strong wind speeds and shallow MLD (Supplement F), an indication of intense mixing. Although mixing events always entrained nutrients into the epilimnion, noticeable as a peak directly after the storm (not shown), the nutrient concentration averaged over the euphotic zone in the first week after the storm could be lower compared to no-storm conditions, especially for moderate wind speeds (5 -10 m/s, Supplement F-1), due to enhanced phytoplankton uptake. If the initial mixed layer was deeper, the effect of the storms on thermal structure and nutrients was less strong (Supplement F-1).    Overall, impacts of identical storms under the RCP8.5 climate scenario did not have a drastically different effect compared to the present-day situation in our simulations ( Figure  4). Both increasing and decreasing effects of storms on chlorophyll-a were found early and later in the stratified season and also with either moderate or high storm intensity ( Figure   4). Storms with a higher intensity tended to have more effect than moderate-intensity storms, but without shifting towards either more positive or more negative effects on chlorophyll-a. Moreover, the difference in effect between high-and moderate-intensity storms was small (Figure 4).

Model validation
In this study, we applied a coupled physical-biogeochemical model, The goodness-of-fit statistics were worst for the chlorophyll-a dynamics. This was largely due to the significant underestimation of the magnitude of the spring peak, and this peak was also simulated too late in some years. In Lake Erken, under-ice growth and resting stages can play an important role in phytoplankton dynamics (Weyhenmeyer et al. 1999;Yang et al. 2016b), and since these processes were not included in the model, this could be an explanation for why the model did not replicate the magnitude of the spring peak.
However, in summer, during stratification, the long-term simulations matched the observed conditions well in most years, except that short-term summer blooms were not well captured by the model. This may have been due to temporary surface accumulations of buoyant cyanobacteria, which were not reproduced by the model.

Causal factors regulating phytoplankton response to storms
Storm intensity, mixed layer depth, surface water temperature, incoming shortwave radiation, and hypolimnetic nutrient concentrations all regulated phytoplankton response to wind perturbations in the first experiment. The variables also interacted with each other.
The random forest permutation importance suggested that light had the strongest effect and hypolimnetic nutrients the least. However, this also depended on the numerical ranges of the variables used in the experiment, which were based on the summer conditions in Lake Erken. In lakes where the summer conditions are substantially different, the order of variable importance may therefore look differently.
The effect of storm intensity was non-monotonic, in the sense that the wind effect had a maximum around 5-10 m/s. It is possible that such moderate wind speeds caused nutrient upwelling without strongly deepening the mixed layer, and as such promoted growth.
Stronger wind speeds, however, may have caused stronger nutrient upwelling, but also more mixed layer deepening, which led to stringent light limitation. This explanation is supported by the volume-averaged nutrient concentrations over the depth of the euphotic zone after storms (Supplement F); strong wind speeds increased nutrient concentrations in the mixed layer, but this was not accompanied by increased phytoplankton concentrations.
Negative influence of wind on phytoplankton concentrations has been shown in various systems (e.g. Fitch and Moore 2007; Kuha et al. 2016;Jalil et al. 2020). However, a positive effect (relative to low wind speed) of moderate wind speed but a negative effect of high wind speeds has, to our knowledge, not often been shown in the literature. Due to the effect of storms on light and nutrient limitation, this non-monotonic relation of wind and phytoplankton could occur more frequently in stratifying lakes. It should be noted that we do not make a definite distinction between what is a storm and what is not, but rather explore the full range between median wind speeds and wind speeds that are higher than experienced in the past 22 years at Lake Erken. The non-monotonic effect of wind speed suggests that this entire wind speed distribution is relevant for lake ecosystem functioning, not just the extreme storms.
Storms had most effect on the investigated lake variables when the mixed layer was around 8 m deep or shallower prior to the storm. In case of a deeper antecedent mixed layer, the mechanical mixing generated at the water surface is largely dissipated at the depth of the maximum density gradient (Imboden and Wüest 1995), so that storms have less effect on thermal, nutrient, and phytoplankton vertical profiles. This threshold of 8 m did not change substantially when we averaged the response variables over a different depth than the euphotic zone depth (results not shown). A modelling study by Mi et al. (2018) also showed that storm-induced mixing has only a small effect on mixed layer depth and entrainment when stratification is deep and strong. While increases in surface temperature and Schmidt stability by climate warming have been reported in multiple studies (e.g. Fang and Stefan 2009;Kraemer et al. 2015;Pilla et al. 2020), trends in depth of stratification are more uncertain and might change only slightly (Pilla et al. 2020). Our simulations indicated that surface temperature and Schmidt stability in Lake Erken will increase under a warmer climate, but that the mixed layer depth will become shallower, similar to what was found by Ayala et al. (2020). A shoaling of the mixed layer would suggest a larger role of storms in the dynamics of stratified lakes.
Whether the net effect of a storm on phytoplankton concentration was positive or negative also depended on the other variables that were varied as part of the first numerical experiment. The effects of nutrients, light, and temperature could be understood from the perspective of nutrient and light limitation. If light conditions are optimal, nutrients are more likely to be limiting, and therefore nutrient upwelling during a storm is likely to promote an increase phytoplankton biomass. On the contrary, if light is the main limiting factor for growth, further deepening of the mixed layer due to mixing is likely to decrease phytoplankton concentrations. The results of the experiment were largely consistent with this. Higher nutrient concentrations in the hypolimnion slightly promoted higher chlorophyll-a concentrations after a storm, but only under high light conditions. The potential of phytoplankton growth being boosted by storm-induced nutrient upwelling in lakes has been reported by observational studies (e.g. Soranno et al. 1997;MacIntyre and Jellison 2001;Crockford et al. 2015;Giling et al. 2017). Increases in phytoplankton after storms in the first numerical experiment were also seen when incoming shortwave radiation was high. Observed decreases in phytoplankton concentration after a storm have been explained by a combination of dilution due to mixed layer deepening and exposure to more stringent and dynamic light conditions (Ibelings et al. 1994;De Eyto et al. 2016;Kuha et al. 2016). Lastly, a higher pre-storm surface temperature resulted in more negative effects of wind on phytoplankton concentration in the first experiment. As water temperature only has a slight positive effect on cyanobacterial growth rates in the model, this result was mostly caused by the effect of surface water temperature on the thermal structure. Since we kept the hypolimnetic temperature constant, a higher surface temperature caused a higher Schmidt stability and therefore stronger stratification. Stronger stratification resisted mixing, as found by Mi et al. (2018), and as such less nutrients ended up in the epilimnion.

Effects of storms on phytoplankton under a warmer climate
Summer-averaged chlorophyll-a concentrations increased in simulations with warmer air temperatures in the second numerical experiment. This is in line with some previous studies, which report increases in phytoplankton with warming (Markensten et al. 2010;Trolle et al. 2014;Gray et al. 2019). However, trends in chlorophyll-a under atmospheric warming also depend on nutrient availability. The present study and the aforementioned studies were done in meso-/eutrophic systems. Under oligotrophic conditions, however, no change or decreasing trends are more likely because of the more stringent nutrient limitation due to earlier onset of stratification and higher nutrient requirements at increased temperatures (Tadonléké 2010;Kraemer et al. 2017).
The projected increases in deep-water nutrient concentrations are also in line with previous studies. Climate warming mediates an increase in hypolimnetic nutrient levels in deep lakes through an increase in hypolimnetic anoxia (Sahoo et al. 2013;North et al. 2014) and incomplete winter mixing (Salmaso 2005;Yankova et al. 2017). Additionally, an earlier onset of stratification separates the epilimnion and hypolimnion earlier in the year, therefore both increasing nutrient concentrations in the hypolimnion and exacerbating nutrient limitation in the epilimnion. Earlier onset of stratification, lower oxygen concentrations, and increases in hypolimnetic nitrate and phosphate in a warmer climate were indeed reproduced by our simulations.
In the second experiment, model predictions suggested that the response of phytoplankton to summer storms would not deviate strongly from the present-day situation when air temperatures were increased to a level consistent with a RCP8.5 climate scenario. Based on the results of the first numerical experiment and projections for summer-average lake conditions in a warmer climate, either stimulating or reducing effects of storms were expected in the second experiment. The reason for this was that of the variables included in the first experiment, the mixed layer tended to become shallower, and hypolimnetic nutrient concentrations and surface water temperature (and therefore strength of stratification) increased as a consequence of atmospheric warming. The opposing effects of these trends in surface temperature and nutrient conditions may have compensated each other. Another potential reason for the lack of difference in the response to wind events between present climate and RCP8.5 could be that even under this high-emission scenario, the changes in the variables were rather small compared to the ranges that were included in the first numerical experiment (i.e. the variation experienced in July at Lake Erken over the past 22 years). Additionally, no strong response would be expected in years when the mixed layer stayed deeper than about 8 m, the depth below which the first numerical experiment showed a marked reduction of storm effects. It follows from the result of numerical experiment 2 that changes in atmospheric warming alone are not likely to strongly change the response of phytoplankton biomass to storms of similar intensity in Lake Erken.

Implications beyond Lake Erken
The present study was set up for Lake Erken, and thus the range of the scenarios (surface temperature, nutrients, solar radiation, mixed layer depth, and wind speed) and morphometry were specific to this lake. The simulated phytoplankton groups (diatoms and cyanobacteria) were intentionally generic and not unique to Lake Erken, but during calibration the phytoplankton parameters were optimised to match the seasonal patterns of the phytoplankton community in Lake Erken. Different phytoplankton communities respond differently to storms (Stockwell et al. 2020). However, the effects of storms on light and nutrient availability occur widely (Stockwell et al. 2020), and the scenarios tested in the first numerical experiment covered a wide range of lake and atmospheric conditions. Therefore, the processes observed in Lake Erken are expected to be similar in other stratifying, mesotrophic lakes. Although the absolute thresholds found in the study may differ from lake to lake and are prone to model uncertainty, the findings in the present study are applicable to more lakes than Lake Erken alone and may facilitate general understanding of lake responses to storms.
Trends in extreme wind speeds (both frequency and intensity) strongly vary with geographic location (Sainsbury et al. 2018). In the area of Lake Erken, future trends in storm intensity and frequency are uncertain (Mölter et al. 2016). However, in regions where storm intensity and frequency are predicted to increase, such as western Europe, an increased importance of storms for phytoplankton dynamics can be expected. Our first experiment showed that the effect of wind on phytoplankton can be non-monotonic, and that moderate wind speeds have different effects than high wind speeds. Therefore, any shift in the probabilitydistribution of wind speeds is relevant, not just trends in extreme wind speeds.
Air temperatures are increasing globally and this causes a rise in strength of stratification, but trends in mixed layer depth remain uncertain (Pilla et al. 2020). In the present study, mixed layer depth was identified as a key variable in responses to storms, so local trends in this variable may partially control storm impacts. Next to air temperature, mean wind speed (Stetler et al. 2021) and water transparency (Read and Rose 2013) determine mixed layer depth, and may change on local or regional scales. In regions that experience atmospheric stilling  or lake brownification (Jennings et al. 2010), mixed layers might shoal and these lakes may be more strongly impacted by storm events. We only scaled air temperature in our climate warming scenario, by applying a change factor correction to a historical time series of air temperature. This allowed the comparison of identical storms under identical pre-event weather conditions, and the differences could be attributed to air temperature alone. However, a different experimental design using the full output from a climate scenario could give an indication of how the full effect of climate change will influence storm-effects on lakes, including trends in, for example, average wind speed.
Nutrient concentrations and their vertical profiles modulate the effect of storms on phytoplankton. Trends in nutrient loading are mostly controlled by human activities, and developing countries especially may experience increasing trends (Fink et al. 2018). Earlier onset of stratification with climate warming (Woolway et al. 2021) causes more nutrient build-up in the hypolimnion, so if nutrients are limiting in the epilimnion, phytoplankton increase after storms may become more prominent, although this was not observed in our second numerical experiment. In lakes where nutrients are high in the epilimnion throughout the year, a wind episode that deepens the mixed layer is likely to decrease phytoplankton concentration due to dilution and reduced light availability.
This study revealed new insights on the effects of storms on phytoplankton, but only certain aspects of this topic were tested. For example, we focused on summer only. In Europe, winter storms tend to be the most severe, but summer storms may have the most impact on a lake by mixing stratified waters  and it is in summer that phytoplankton blooms occur most frequently. Additionally, we focused solely on wind.
Passing storms tend to affect not only wind speed, but also air temperature and incoming solar radiation. Moreover, precipitation during storms will affect both the quantity and quality of catchment runoff. The retention time of Lake Erken is around 7 years, indicating that catchment runoff during storms is likely to have a minor influence, but in lakes with short residence times (e.g. < 1 year) this change in inflow can be at least as impactful as wind (Klug et al. 2012;Reichwaldt and Ghadouani 2012;De Eyto et al. 2016). In addition, in the first experiment, only a limited number of variables that could affect the lake state prior to and after an event were assessed. Like for the effect of precipitation, the results in the present study do not discard the possibility that such other variables are important too.
The effect of phytoplankton community composition on the response to storms was not systematically explored in the present study. We only calibrated and validated total chlorophyll-a concentration, despite using two separate phytoplankton groups in the model.
The diatoms dominated the spring peak and the summer peak contained more cyanobacteria, which was in line with the seasonal dynamics at Lake Erken. The chlorophylla data was more readily available and at higher frequency than taxonomic data. In order to assess the effect of storms on individual groups, the model would have to be calibrated and validated on group-specific data, and parameter values would have to be informed by a trait-based approach. Traits such as buoyancy regulation, nutrient storage, nutrient acquisition, growth rates, and photoadaptation are likely to be highly relevant for how a phytoplankton community responds to a storm (e.g. Visser et al. 1996;Kasprzak et al. 2017;Stockwell et al. 2020). As such, different communities may show distinct responses to storms under physically and chemically comparable situations. If storms already occur frequently in a system, the phytoplankton community may be adapted to such conditions (Stockwell et al. 2020), and therefore shifts in frequency of storms may be especially relevant when considering community composition. Acquisition of data of sufficient frequency and quality to investigate these topics is a key problem, but a combination of experiments, use of novel monitoring techniques, and modelling may elucidate some of this uncertainty in the near future. Also, our model results pointed at average concentrations and responses, not at phytoplankton bloom or scum formation. Prediction of blooms with data-driven or process-based models still remains a challenge (Rousso et al. 2020), and may require inclusion of processes that are not parameterised in our model, such as life cycles (Hense and Beckmann 2010) or selective grazing by zooplankton (Sommer et al. 2012). This may be part of the reason for why the spring peak and occasional summer spikes in chlorophyll-a were missed by the model used in this study.

Conclusion
High wind speeds (⪆ 10 m/s) always had more negative effects than moderate wind speeds (≈ 5 -10 m/s), but the direction of the effect depended mostly on the level of incoming radiation, surface water temperature, and hypolimnetic nutrients. The effect of storms decreased markedly when the mixed layer depth was about 8 m or deeper. Higher incoming radiation and hypolimnetic nutrient concentrations promoted increases in chlorophyll-a concentrations after storms, whereas increases in surface temperature had a decreasing effect. These outcomes confirmed the conflicting effects of storms on phytoplankton light and nutrient limitation, and provide a mechanistic framework to better understand under what conditions storms tend to either increase or decrease phytoplankton biomass. A simulation forced by a future climate scenario showed earlier onset of stratification and a higher summer chlorophyll-a concentration, averaged over the euphotic zone. However, the response of phytoplankton to storms did not strongly change with warming air temperatures.
Increased understanding of the drivers of storm impacts on lakes can help short-term forecasting, and in some cases may be used to inform lake or reservoir management.
Additionally, it facilitates assessment of how atmospheric trends will affect lakes, specifically those caused by climate change. Different regions are expected to experience different trends in air temperature, (extreme) wind speed, and nutrient loading. Studies evaluating the combined effects of these trends to assess the impacts of storms on lake phytoplankton could further our understanding of the global impact of extreme weather events on lake ecosystems.