Drivers of phytoplankton responses to summer wind events in a stratified lake: A modeling study

Extreme wind events affect lake phytoplankton by deepening the mixed layer and increasing internal nutrient loading. Both increases and decreases in phytoplankton concentration after strong wind events have been observed, but the precise mechanisms driving these responses remain poorly understood or quantified. We coupled a one‐dimensional physical model to a biogeochemical model to investigate the factors regulating short‐term phytoplankton responses to summer wind events, now and under expected warmer future conditions. We simulated physical, chemical, and biological dynamics in Lake Erken, Sweden, and found that strong wind could increase or decrease the phytoplankton concentration in the euphotic zone 1 week after the event, depending on antecedent lake physical and chemical conditions. Wind had little effect on phytoplankton concentration if the mixed layer was deep prior to wind exposure. Higher incoming shortwave radiation and hypolimnetic nutrient concentration boosted phytoplankton concentration, whereas higher surface water temperatures decreased concentrations after wind events. Medium‐intensity wind events resulted in more phytoplankton than high‐intensity wind. Simulations under a future climate scenario did not show marked differences in the way wind events affect phytoplankton concentration. These findings help to better understand how wind impacts vary as a function of local environmental conditions and how climate warming and changing extreme weather dynamics will affect lake ecosystems.

High wind speeds reshape lake physical and chemical environments in ways that can alter phytoplankton growth rates. At the start of a chain of processes, wind stress at the lake surface induces internal mixing, which deepens the thermocline (Andersen et al. 2020). Mixing affects the vertical distributions of oxygen and nutrients, and in turn phytoplankton growth and vertical distribution. Upwelling of nutrient-rich water during mixing 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 slightly decrease during storms (Mesman et al. 2020;Doubek et al. 2021), potentially reducing growth rates of light-saturated phytoplankton (Trombetta et al. 2019). In addition, 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, wind storms have conflicting effects on nutrient and light *Correspondence: jorrit.mesman@unige.ch This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Additional Supporting Information may be found in the online version of this article.
Author Contribution Statement: J.P.M. drafted the manuscript and did the model simulations. All authors were involved in the conceptualization of the study and reviewed the manuscript.
availability (Stockwell et al. 2020), and the net effect of a wind event 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 wind. Extreme wind events will likely shift in frequency and intensity as a result of climate change, with different regions of the globe experiencing increases or decreases (Mölter et al. 2016;Sainsbury et al. 2018). Concurrent with these changes in meteorological conditions, more energy will be needed to mix the water column as surface water temperatures increase and stratification strengthens (Schmidt 1928;O'Reilly et al. 2015;Pilla et al. 2020). Therefore, in a warmer climate a wind event of a given magnitude and duration may cause less mixing in stratified lakes than in the current climate. The depth of the mixed layer, at the time when a wind perturbation hits, also determines the degree of lake mixing; a deeper pre-event 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 no change (Kraemer et al. 2015;Pilla et al. 2020). Furthermore, local trends in transparency or wind speed are likely to be 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. This will lead to 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 wind event and become available for phytoplankton growth. The combined effects of climate warming on the strength and depth of stratification and on nutrient profiles may therefore cause the response of phytoplankton to high wind speeds to either increase or decrease in a future climate.
Understanding mechanistically how complex lake ecosystems are reshaped by wind events under present and projected future conditions is a challenging task. Process-based models have the advantage of enabling a quantitative comparison between different scenarios and identifying clear causal pathways even in complex systems and under 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. Moreover, extreme events are rare by definition, hard to predict, and thus difficult to study. Extreme winds also act on short timescales (hours to days), while lake monitoring programs often include only weekly or monthly sampling. Therefore, biogeochemical data describing responses to wind events are scarce. In the present study, we used the General Ocean Turbulence Model (GOTM) model coupled to a biogeochemical model to study the impact of wind events on phytoplankton dynamics. Onedimensional process-based hydrodynamic models (including GOTM) can simulate physical effects of extreme wind with reasonable accuracy (Mesman et al. 2020), providing some confidence that processes related to the transport and mixing of biogeochemical particles and solutes during wind events can be accurately simulated. Our approach can help identify potential regulating factors for the responses of phytoplankton to wind events, and how climate warming may affect these responses.
We used Lake Erken, a mesotrophic dimictic lake in Sweden, as a case study, because of the available long-term time series of physical and biological variables that we used to calibrate our model. However, the findings of this study are generalizable to the processes regulating phytoplankton responses to wind across stratifying lakes. Phytoplankton communities in such lakes are shaped by their need for both nutrients and light, which exhibit opposing gradients with lake depth. We assessed scenarios covering a broad range of atmospheric and lake conditions, which reflect conditions present in many temperate, stratifying lakes.
We investigated (1) how wind speed, thermal structure, light availability, and nutrient availability control phytoplankton response to wind events during summer stratification, and (2) how climate warming may influence the response of phytoplankton to wind events. We performed two numerical experiments to address these topics. In the first experiment, we repeatedly simulated a wind event while changing wind speed, incoming shortwave radiation, and preevent MLD, surface water temperature, and hypolimnetic nutrient concentration in a full factorial design. In the second numerical experiment, we compared the response of phytoplankton to wind perturbations between present-day and future-climate air temperatures, at different times of the year and at different intensities of wind events. These simulations may help to disentangle and better understand the dynamic response of phytoplankton to wind events in a changing world.

Site description
Lake Erken (Fig. 1) is a mesotrophic lake in Sweden (59 50 0 37 00 N, 18 35 0 38 00 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 yr (Blenckner et al. 2002). Lake Erken is dimictic, 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 (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 are circulated through the water column after the autumn turnover .
In most years, Lake Erken experiences a distinct spring phytoplankton bloom followed by a clear water phase, and then a second phytoplankton peak in summer or autumn (Yang et al. 2016b), similar to many monomictic and dimictic lakes across the globe (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, the cyanobacterium Gloeotrichia echinulata is a major bloom-forming species (Karlsson-Elfgren et al. 2003;Yang et al. 2016b).

Data collection
We used meteorological, water temperature, and biogeochemical data from 1999 to 2020. Meteorological data (wind speed, air temperature, air pressure, relative humidity, shortwave radiation, cloud cover, and precipitation) were collected every hour at a weather station on a small island in the lake (Fig. 1). Moras et al. (2019) replaced missing meteorological data with 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 estimated by the HYPE model (Hydrological Predictions for the Environment, Lindström et al. 2010) from 2004 onwards, and validated using measured data from the main tributary of the lake. Discharge and temperature data in the main tributary were automatically monitored and summarized to daily values, while phosphate, total phosphorus, nitrate, and particulate organic matter concentrations in the inflow were measured once or twice per month. Because the discharge data were available only since 2004, the first 5 yr of the discharge data 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 thermistor chain roughly 500 m northeast of the weather station, every 0.5 m down to 15 m depth (Fig. 1). Starting in 2017, these data were collected year-round. Schmidt stability (Schmidt 1928;Idso 1973) and MLD were calculated from the water temperature profiles. Schmidt stability was calculated using the "rLakeAnalyzer" R package (Winslow et al. 2019) and the MLD was defined as the depth where water density had increased by 0.15 kg m À3 relative to the surface (similar method as in Wilson et al. 2020).
Water samples to determine nutrient concentrations (phosphate, total phosphorus, nitrate, and ammonium) and chlorophyll a (Chl a) were collected above the deepest point of the lake (21 m), close to the site of the weather station, every 2 weeks during the ice-free season. During stratification, separate integrated samples of the epilimnion and hypolimnion were taken. The integrated samples were created by (1) determining the depth of the thermocline from temperature and oxygen profiles, (2) taking water samples at 2-m intervals using a tube sampler, and then (3) creating a mixed sample for epilimnion and hypolimnion separately, weighting each sample according to the lake bathymetry. In winter, if the ice was sufficiently safe for sampling, a single integrated sample was taken through a hole in the ice approximately every month. Nutrients were analyzed using standard laboratory techniques (Ahlgren and Ahlgren 1976;Goedkoop and Sonesten 1995). Chl a was concentrated by filtration on glass fiber filters and then analyzed using spectrophotometry (Ahlgren and Ahlgren 1976).

Model description and setup
GOTM is a one-dimensional 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, the biogeochemical equations are applied to each layer in GOTM, including water-atmosphere and water-sediment exchange, and GOTM regulates the transport of biogeochemical substances (e.g., oxygen and nutrients) between the layers. Using FABM, GOTM was coupled to a modified version of the SELMA model, which 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 (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, mineralization of detritus, phytoplankton growth regulated by nutrients and light, and grazing by zooplankton (Neumann et al. 2002). Chl a content and phytoplankton biomass are linked through a fixed chlorophyll-to-carbon ratio, which was required to match our observations of Chl a to the biomass estimates of the model. We modified the SELMA model code by (1) expressing phytoplankton and zooplankton concentration 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 Chl a concentration directly to the carbon content 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 setup to other biogeochemical models and a more complete description of potentially relevant processes. A detailed common-language description of the model has been supplied in Supplement S1 of Supporting Information and the model code is publicly available (see Software and Data 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 h and 0.5 m thick layers. A fourth-order Runge-Kutta time integration scheme was used in 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 to run the model, see Software and Data Availability.

Calibration
After 5 yr of spin-up, 13 yr (2004-2016) were used for calibration. We used the parsac software (Table 1, Bruggeman and Bolding 2020), employing a differential evolution method to optimize the maximum likelihood objective function of the root mean square error (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 were optimized using 10,000 iterations, by calibrating five parameters based on a previous study Table 1. Software used or referred to in this study. Supplement S1 of Supporting Information contains more information on how GOTM and Selmaprotbas are coupled to each other, and how the Selmaprotbas model was derived from the SELMA and PROTBAS models.

Abbreviation
Full name Description References  (Ayala et al. 2020): 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 run, 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, Chl 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 AE 10% of the value obtained in the first calibration step. Because the nutrient and Chl a samples were based on integrated samples for the epilimnion and hypolimnion, for the calibration we assumed these samples to be representative of 3 and 15 m depth, respectively. The parsac calibration attaches equal weight to each observation. 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 S2 of Supporting Information.

Sensitivity analysis
We performed a sensitivity analysis using the parsac software to determine what parameters to include in the calibration of the Selmaprotbas model (Bruggeman and Bolding 2020;Andersen et al. 2021). Parsac applies the Sensitivity Analysis Library (Herman and Usher 2017) in Python. All parameters in the Selmaprotbas model, scaling factors for inflow discharge and concentrations, and the five calibrated parameters in GOTM were included in the sensitivity analysis (totaling 55 parameters), and responses in simulated mean values of ammonium, nitrate, phosphate, total phosphorus, Chl a, and oxygen were assessed.
We followed a density-based delta-sensitivity method (Borgonovo 2007), which has been described earlier for the GOTM-FABM-PCLake model (Andersen et al. 2021). First, a Latin hypercube sampling (McKay et al. 2000) was applied to generate a number of parameter sets equal to 200 times the number of parameters. This number of parameter sets was based on a convergence test, where we found that values for sensitivity started to converge at this sampling size (results of convergence test not shown). We then ran the model for all parameter sets in the Latin hypercube and applied the deltasensitivity analysis to the results. 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% confidence interval around the sensitivity values was estimated 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% confidence interval of the dummy for all variables (ammonium, nitrate, phosphate, Chl 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 S2 of Supporting Information.

Validation
We calculated RMSE, mean absolute error (MAE), mean error (ME), and Nash-Sutcliffe efficiency (NSE, Nash and Sutcliffe 1970) from the simulated and measured values of water temperature, nutrients, oxygen, and Chl a as measures for the goodness-of-fit. All measured values at each depth were compared to the corresponding simulated values, by linear interpolation of simulated values if necessary. The last 4 yr of the time series (2017-2020) were used for validation, and thus not used to train the model. The measures for the goodness-of-fit were compared between the calibration and the validation period to assess the quality of the simulation. The simulation from 1999 to 2020 using the observed weather conditions and the calibrated parameter values is 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: Variables controlling phytoplankton response to wind events
Our first numerical experiment investigated which variables control the response of phytoplankton concentrations to wind events. Consequently, we induced a 1-d wind event for different values of several meteorological and pre-event lake variables. Five variables were chosen that were expected to impact the phytoplankton response to wind: • Wind speed during the event: Represents the magnitude of the disturbance induced by the event. • MLD prior to the wind event: A measure of the sensitivity of the thermal structure to further deep-water entrainment and controls the volume of the epilimnion. • Shortwave radiation during the event: Regulates the availability of light. • Surface water temperature: Increases the strength of stratification prior to the event, and therefore the resistance to mixing. • Hypolimnetic nutrients: Regulate the potential for nutrient upwelling.
Ten levels of each variable were taken into consideration in a full factorial design, therefore totaling 10 5 simulations. The wind perturbation had a duration of 24 h and was initiated 24 h after initialization of the model run (Fig. 2a). We focused on the month of July to generate the weather conditions for this numerical experiment. Lake Erken was always stratified in July and stratification would have existed for sufficient time to allow for the build-up of hypolimnetic nutrients. To have a representative period with natural weather variations, we selected a year in our dataset with the most generic weather conditions during July. 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 (c) incoming shortwave radiation, and (d) cloud cover. The solid lines were used in all simulations and the dashed lines indicate the different scenarios used in the numerical experiment, which were varied independently of one another. In panels (e-g), the initial vertical profiles of (e) water temperature, (f) nitrate, and (g) phosphate are shown. The initial profiles varied in both MLD (dashed lines) as well as the epilimnetic value for water temperature and the hypolimnetic value for nitrate and phosphorus (solid lines). Because the experiment had a full factorial design, each of the panels (e-g) could contain 100 lines (varying both MLD and absolute value), but instead variation in either MLD (keeping the absolute value constant) or absolute value (keeping MLD constant) is shown. Nitrate and phosphate were part of the same scenario (i.e., "nutrients") and therefore did not vary independently. These initial profiles were based on the average conditions during July in the long-term simulations. The values of the levels used in the scenarios can be found in Supplement S3 in the Supporting Information. cloud cover) for each year separately and for the full period . For each year, we then calculated the root mean squared relative error (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 we used these weather conditions as baseline in the first experiment ( Fig. 2a-d).
For all model variables, we based initial conditions on the average long-term July profiles for the epilimnion and hypolimnion separately. Both profiles were defined as 10 equidistant values, interpolated from surface to the MLD for the epilimnion, and from the MLD to the maximum lake depth for the hypolimnion. As an example, when the MLD was 8 m, the epilimnetic profile consisted of 10 values, a linear interpolation of the simulated values over 0-8 m, and when the MLD was 4 m, the epilimnetic profile still had 10 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 , and dates with an MLD < 2 or > 15 m were excluded. 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 Fig. 2e-g).
The 10 different levels of wind speed, MLD, incoming solar radiation, surface temperature, and hypolimnetic nutrient concentration used in the simulations were also taken from the long-term simulation (1999-2020) using data from July. The 5 th and 95 th percentiles of the daily averages were computed, and a linear sequence of 10 values between these percentiles was used for the experiment (Fig. 2). The only exception was wind speed. Because 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 5 m s À1 above the maximum recorded daily average wind speed, to anticipate the possibility of increased extreme wind speeds in the future due to climate change (Mölter et al. 2016). Moreover, energy transfer from wind to the lake increases exponentially with wind speed (Wüest et al. 2000). Because we used a constant wind speed for 24 h (Fig. 2), our simulated wind events would cause less mixing than a brief event with the same daily mean but varying wind speeds. Using a higher wind speed than the recorded maximum partially compensated for this. Each combination of the independently changed variables was a separate simulation.
The model output of each run-concentrations of Chl a, nitrate, and phosphate-was volume-averaged over the euphotic zone for the first week after the end of the wind event. 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 a euphotic zone of fixed depth, 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 mixed layer deepening on production and Chl a. In addition to the average concentrations, average Schmidt stability was also calculated for the first week after the event.
As the last step of this experiment, we fitted the volumeaveraged Chl 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. All data from the experiment were used (no holdout), because we wanted to identify the most important variables rather than make predictions. 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).
Numerical experiment 2: Effect of a warmer climate on phytoplankton response to wind events Whereas the first numerical experiment investigated the effect of individual variables, the second experiment focused on discerning the net effect of a warming climate on phytoplankton response to wind events. The long-term simulation, from 1999 to 2020, under observed weather conditions, was taken as the baseline. Every year (N = 22) a 24-h wind perturbation was added to this baseline. To avoid cumulative effects, only one perturbation was applied to each 22-yr simulation. For the rest of the period, observed weather data were used. Based on average seasonal patterns of stratification, this perturbation occurred early summer (9 th of June, corresponding to the first day-of-the-year when the averaged Schmidt stability over all the years exceeded 50% of the maximum) or midsummer (4 th of July, 80% of maximum Schmidt stability), and with moderate (7.2 m s À1 , i.e., 95 th percentile of daily average wind speed) or high (9.0 m s À1 , 99 th percentile) intensity (Table 2; for the determination of these thresholds, see Supplement S4 of Supporting Information).
We repeated our design for a climate scenario in which the air temperature was scaled to the period 2041-2070 of an Representative Concentration Pathway (RCP) 8.5 emission scenario (IPCC 2014) of the regionally downscaled output of the HadGEM2-ES global climate model (Collins et al. 2008). This climate model output was created as part of the EURO-CORDEX experiment ( Table 2, Jacob et al. 2013). The measured Lake Erken air temperature was scaled according to the delta-decile method (Perroud and Goyette 2010); for every month and every decile of daily averaged air temperature (0 th to 10 th percentile, 10 th to 20 th percentile, etc.) the increase in temperature was calculated and applied to the historical time series. The other meteorological conditions, including relative humidity and wind speed, were kept the same. We chose this approach to isolate the effects of warming alone and to draw conclusions for a wider range of lakes, as local trends in other variables (such as cloud cover) 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 MLD and lead to an earlier onset of stratification and thus to different vertical profiles of nutrients and oxygen (Kraemer et al. 2015;Woolway et al. 2021). As such, increased air temperature influences multiple potentially important lake variables that affect the relationship between wind events and phytoplankton, which were investigated in isolation in the first numerical experiment.
In each simulation, the volume-averaged Chl a concentration between 0 and 8.0 m depth (the euphotic zone) was calculated. The maximum difference in concentration between the control and the perturbation scenarios (Table 2), within 1 wk after the event, was determined. These values were then compared between the present-day and the future climate scenarios.

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 Chl a was the substantial underestimation of the spring Chl a peak (Supplement S5 of Supporting Information). In the validation period, the model fit worsened slightly for phosphate, nitrate, and Chl a, as indicated by the NSE values. Inspection of the time series (Supplement S5 of Supporting Information) 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. Chl a concentrations showed distinct spring and summer peaks, but the spring peak concentrations were underestimated in almost all years, and sometimes simulated too late. Summer Chl a levels tended to be close to observed levels, with the exception of some summer blooms, when the model underestimated observed values. Simulated summer epilimnetic concentrations of both nitrate and phosphate were low, 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.
Numerical experiment 1: Variables controlling phytoplankton response to wind events MLD, surface water temperature, hypolimnetic nutrient concentration, and incoming shortwave radiation all interacted with wind speed to affect phytoplankton response to wind events. Wind speed had a nonlinear effect on the average phytoplankton concentration in the first week following the event: moderate wind speeds increased Chl a concentrations, but strong winds (on the order of 10 m s À1 or higher) had a Table 2. Design for the second numerical experiment. Starting from the long-term simulation (1999-2020), 24-h wind perturbations were added either early or later in summer, and with moderate or high intensity. These wind perturbations were added to every year in separate simulations, to avoid accretion of effects. This design was repeated for a future-climate scenario, in which air temperatures were scaled to the level of 2040-2070 according to an RCP 8.5 climate scenario.  Table 3. RMSE, MAE, ME, and NSE for the calibration ("Cal.," 2004-2016) and validation ("Val.," 2017-2020) periods. For RMSE, MAE, and ME, values close to 0 indicate an optimal fit, whereas for NSE a value close to 1 indicates an optimal fit. These metrics are calculated for the epilimnion in case of nitrate, phosphate, and Chl a, and for the full water column for oxygen and temperature. less positive, or even a reducing effect on Chl a concentration (Fig. 3b, inset). When the initial MLD was deeper than about 8 m, even strong winds had substantially less effect on Chl a concentrations than under shallower mixed layers (Fig. 3b).

Variable
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 daily averaged), the effect of wind on phytoplankton was largely negative, although at moderate wind speeds some increases were still evident (Fig. 3b). When incoming radiation increased above 200 W m À2 , wind had an overall positive effect on phytoplankton concentration in the euphotic zone, moderate wind speeds more so than high wind speeds. Hypolimnetic nutrients were less influential than light or surface water The subpanels each represent a scenario with a different hypolimnetic nutrient concentration (low on the left, high on the right; N and P concentrations are indicated on top) and a different average incoming shortwave radiation (low on top, high on bottom, values indicated on the left). Each subpanel has the wind speed during the event on the x-axis and the MLD on the y-axis (see labels in the inserts), and each small rectangle inside each subpanel corresponds to one simulation. The change in main panel (b) was calculated relative to the simulation with a wind speed of 3.7 m s À1 (which is why the left column of each subpanel always has a value of 0). The numerical experiment was performed using 10 levels of shortwave radiation and nutrient levels. Only four by four subpanels are shown for visualization purposes, because the results for intermediate input values can be interpolated from the ones presented. temperature (Table 4), but higher concentrations resulted in increased phytoplankton concentration following wind events at mixed layers deeper than about 4 m when incoming radiation was high (Fig. 3b). At the shallowest mixed layer (2 m), wind speed had a negative effect on Chl a when hypolimnetic nutrients were high. Surface water temperature before the event (and therefore Schmidt stability) also influenced phytoplankton response to wind; a higher initial surface temperature caused a decrease of phytoplankton after a wind perturbation (Supplement S6 of Supporting Information).
Apart from changes in Chl a, we also examined changes in Schmidt stability, and nutrient concentrations. The strongest decrease in Schmidt stability was diagnosed at strong wind speeds and shallow MLD (Supplement S6 of Supporting Information), an indication of intense mixing. Although mixing always entrained nutrients into the epilimnion, noticeable as a peak directly after the wind event (not shown), the nutrient concentration averaged over the euphotic zone in the first Table 4. Permutation-based variable importance of wind speed, MLD, hypolimnetic nutrients, incoming shortwave radiation, and surface water temperature for predicting change in volume-averaged Chl a over the euphotic zone (i.e., the values shown in Fig. 3b) in the first numerical experiment, based on a fitted random forest model. The out-of-the-bag R 2 was 0.992. The importance values represent the degree to which the mean squared error of the fitted random forest model increases if the input column for that variable is randomly permuted. The greater the error after permutation, the more the model relies on the variable for making predictions, and therefore the more important the variable is to the model. week after the event could decrease due to enhanced phytoplankton uptake, especially for moderate wind speeds (5-10 m s À1 , Supplement S6 of Supporting Information). If the initial mixed layer was deeper, the effect of the wind on thermal structure and nutrients was less strong (Supplement S6 of Supporting Information).

Permutation-based importance
Numerical experiment 2: Effect of a warmer climate on phytoplankton response to wind events Under the warming scenario, Lake Erken surface temperatures in summer increased by roughly 1.6 C while deep-water temperatures remained similar to present-climate conditions, and therefore the Schmidt stability increased (increase in median by 43%, Fig. 4). The mixed layer shoaled by 1.2 m (14%), and stratification tended to both form earlier and vanish later in the year (medians are 3.5 d earlier and 6.5 d later in the RCP8.5 scenario, respectively). Median deep-water oxygen concentrations decreased by 20% (1.17 mg L À1 ), while deep-water nitrate and phosphate concentrations increased (nitrate increase in median of 50%, 14.4 μg L À1 ; phosphate increase in median of 28%, 8.8 μg L À1 ). Average summer Chl a concentration tended to be higher in the warming scenario compared to present-day conditions (0.56 μg L À1 higher, 10% increase, Fig. 4).
Overall, impacts of identical wind events under the RCP8.5 climate scenario did not have a drastically different effect compared to the present-day situation in our simulations (Fig. 5). Both increasing and decreasing effects of wind events on Chl a were found early and later in the stratified season and also with either moderate or high intensity (Fig. 5). Events with a higher wind speed tended to have more effect than moderate-intensity events, but without shifting toward either more positive or more negative effects on Chl a. Moreover, the difference in effect between high-and moderate-intensity wind events was small (Fig. 5).

Discussion
We applied a coupled physical-biogeochemical model, GOTM-Selmaprotbas, to investigate the effect of wind events on lake phytoplankton under present and projected conditions in Lake Erken, Sweden. Moderate wind speeds (≈ 5-10 m s À1 ) promoted phytoplankton concentration more strongly than high wind speeds and also had increasing effects compared to low wind speeds when there was sufficient incoming light. This shows that the effect of wind events on phytoplankton is nonlinear. The effect of wind on phytoplankton decreased when MLD exceeded 8 m, highlighting the effect of thermal structure on the susceptibility of phytoplankton to wind perturbations. Higher availability of light and nutrients caused a larger increase of phytoplankton after wind events, whereas a stronger stratification decreased phytoplankton concentrations after wind perturbation. In the second numerical experiment, we found that under warmer atmospheric conditions, the response of phytoplankton to wind events is similar to that under present-day climatic conditions. These findings help to improve our understanding of how wind events influence phytoplankton and how climate histogram represents the maximum difference in volume-averaged Chl a concentration in the euphotic zone between simulations with and without wind events on the y-axis. Positive values indicate that wind events had an increasing effect on Chl a concentration, and negative values that the event had a decreasing effect. The back-to-back histograms compare the present-climate (bars to the left, in light gray) against the RCP8.5-scaled climate scenario (bars to the right, in black). The four panels compare these outcomes in simulations when wind events are moderate or severe (left-and right-hand panels, respectively), and when they occur early in summer (9 th of June) vs. later in the summer (4 th of July; top vs. bottom panels). Each panel represents the outcome of multiple simulated events, one for each year in the long-term simulation (N = 22). On the x-axis are the number of simulated events that showed the corresponding change in Chl a compared to baseline conditions. change (including changing extreme weather patterns) will affect lakes.

Model validation
The model successfully reproduced the seasonal cycles of all variables. Water temperature was simulated well, and with an RMSE of 0.9 C the fit was similar to those in previous studies in Lake Erken (Moras et al. 2019;Ayala et al. 2020). Periods of summer hypolimnetic hypoxia were captured well by the model in most years, although the clear underestimation of under-ice oxygen consumption did point toward an incomplete description of winter oxygen dynamics by the model. Also, the dynamics of dissolved nutrients (nitrate and phosphate) were reproduced rather well, with concentrations close to observed values in summer, the period under study. However, winter concentrations of phosphate and nitrate tended to be underestimated and the increase in nitrate in autumn was simulated too early, potentially due to a missing recalcitrant organic matter component, such as macrophytes near the shore. In addition, expanding the model description of sediment dynamics may lead to improvements in oxygen and nutrient simulations. The statistics describing the nutrient simulations for our simulations of Lake Erken are similar to recent applications of the GLM-AED2 model in Lake Mendota, USA, (Ladwig et al. 2021) and of the CE-QUAL-W2 model in Rappbode Reservoir, Germany, (Mi et al. 2020).
The goodness-of-fit statistics were worst for the Chl 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). These processes were not included in the model, which could explain why the model did not replicate the magnitude of the spring peak. Moreover, as chlorophyll-to-carbon ratios tend to be higher during the high-nutrient conditions in winter/spring (Riemann et al. 1989;Jakobsen and Markager 2016), the use of a fixed ratio in our model may have partially caused the underestimation of the spring peak. However, in summer, during stratification, the long-term simulations matched the observed conditions well in most years, which supports the use of the model in the numerical experiments, that were conducted during the summer months. Short-term changes in the monitored Chl a concentrations, however, were not well captured by the model, as the model points at average concentrations and responses. Temporary surface accumulations of buoyant cyanobacteria or sudden increases in pelagic biomass due to emerging resting stages would not be simulated, which may have been the cause of the missed summer phytoplankton peaks.
The choice for a one-dimensional model enabled the high number of scenarios used in this study. However, this also meant that wind direction and associated three-dimensional processes that cannot be resolved in a one-dimensional model (e.g., near-shore upwelling) were not included in the model, despite their potential importance in how extreme wind affects lakes (MacIntyre and Jellison 2001;Roberts et al. 2021). During calibration, a one-dimensional model partially compensates for the lack of description of such processes, but a complete simulation would require a three-dimensional model. Especially lakes with a complex bathymetry may therefore show different mixing dynamics in response to wind events than assessed here.

Causal factors regulating phytoplankton response to wind events
As mentioned before, the effect of wind speed was nonmonotonic, with a maximum increase in phytoplankton concentration after wind events of around 5-10 m s À1 . Such moderate wind speeds may cause nutrient upwelling without strongly deepening the mixed layer, and as such promoted growth. Stronger wind speeds, however, caused stronger nutrient upwelling, but also more mixed layer deepening and consequently light limitation (Kuha et al. 2016;Jalil et al. 2020). This explanation is supported by the volume-averaged nutrient concentrations over the depth of the euphotic zone after the wind perturbations (Supplement S6 Supporting Information); strong wind speeds increased nutrient concentrations in the mixed layer, but this was not accompanied by increased phytoplankton concentrations. A positive effect (relative to low wind speed) of moderate wind speed combined with a negative effect of high wind speeds has, to our knowledge, not yet been shown in the lake literature. In marine studies, a maximum phytoplankton concentration was found at moderate wind speeds (Millet and Cecchi 1992;Fitch and Moore 2007), similar to our study. This nonmonotonic effect of wind speed suggests that the entire wind speed distribution is relevant for aquatic ecosystem functioning, not just the extremes.
Another result was that wind had the largest effect on the investigated lake variables when the mixed layer was around 8 m or shallower prior to the event. 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 wind has less effect on thermal, nutrient, and phytoplankton vertical profiles. The 8-m threshold did not change substantially when we averaged the response variables over a shallower or deeper depth than the euphotic zone depth (results not shown). Wind-induced mixing also had only a small effect on MLD and entrainment when stratification was deep and strong in a modeling study by Mi et al. (2018).
Whether the net effect of a wind event on phytoplankton concentration was positive or negative also depended on the other variables that were varied as part of the first numerical experiment. Higher nutrient concentrations in the hypolimnion promoted higher Chl a concentrations in the euphotic zone after the event, but only under high incoming light conditions. Wind-induced nutrient upwelling is indeed known to potentially boost phytoplankton growth in lakes (Soranno et al. 1997;MacIntyre and Jellison 2001;Crockford et al. 2015;Giling et al. 2017). The negative effect of wind on Chl a at the shallowest MLD at high nutrients was caused by the setup of the experiment (Fig. 2), as even under the lowest wind speed a large amount of nutrients entered the epilimnion, and high wind speeds therefore mostly reduced growth due to mixed layer deepening. Phytoplankton concentration after wind events also increased in the first numerical experiment 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-event 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: because we kept the hypolimnetic temperature constant, a higher surface temperature caused a higher Schmidt stability and therefore stronger stratification. Stronger stratification resists mixing (Mi et al. 2018), and as such less nutrients ended up in the epilimnion.

Effects of wind events on phytoplankton under a warmer climate
Summer-averaged Chl a concentrations increased in simulations with warmer air temperatures in the second numerical experiment, consistent with previous studies in mesotrophic/ eutrophic systems (Markensten et al. 2010;Trolle et al. 2014;Gray et al. 2019). The projected increases in deep-water nutrient concentrations are also in line with previous studies. Climate warming promotes 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). In addition, 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 wind events would not deviate strongly from the present-day situation when air temperatures were increased to a level consistent with a RCP8.5 climate scenario. The first experiment suggested that surface temperature increases would cause lower phytoplankton concentrations after wind events in the RCP8.5 scenario and that higher hypolimnetic nutrient concentrations would cause higher phytoplankton concentrations. The opposing effects on phytoplankton of these trends may have largely canceled each other out in the second experiment, resulting in no marked difference in the response of phytoplankton to wind events between the present-climate and RCP8.5 scenarios. Another potential reason for these relatively small differences could be that even under this high-emission scenario, the changes in the lake 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 yr). In addition, 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 wind effects. Thus, changes in atmospheric warming alone are not likely to strongly change the response of phytoplankton concentration to wind events of similar intensity in Lake Erken.

Implications beyond Lake Erken
The range of the scenarios (surface temperature, nutrients, solar radiation, MLD, and wind speed) and morphometry were specific to Lake Erken. The simulated phytoplankton groups (diatoms, dominant in spring, and cyanobacteria, dominant in summer) were intentionally generic and not unique to Lake Erken, but during calibration the phytoplankton parameters were optimized to match the seasonal patterns of the phytoplankton community in this lake. However, the effects of extreme wind events 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 we found are likely to differ from lake to lake and are prone to model uncertainty, our findings are applicable to other lakes as well and may facilitate general understanding of lake responses to wind events.
Trends in extreme wind speeds (frequency, duration, and intensity) vary with geographic location (Webster et al. 2005;Sainsbury et al. 2018). Future trends in storm intensity and frequency in the region surrounding Lake Erken are uncertain (Mölter et al. 2016). However, in regions where storm intensity and frequency are predicted to increase, such as western Europe, the importance of wind for phytoplankton dynamics may increase.
Air temperatures are increasing globally and are expected to strengthen stratification, but trends in MLD remain uncertain (Pilla et al. 2020). In the present study, MLD was identified as a key regulator in responses to wind, so local trends in this variable may affect impacts of wind events. Next to air temperature, mean wind speed (Stetler et al. 2021) and light penetration (Read and Rose 2013) determine MLD, and these 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 become more strongly impacted by wind events. Nutrient concentrations and their vertical profiles modulate the effect of wind events 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 and a later overturn with climate warming (Woolway et al. 2021) causes more nutrient build-up in the hypolimnion (Isles et al. 2017), so if nutrients are limiting in the epilimnion, phytoplankton increase after wind events may become more prominent. However, 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.
We focused solely on wind, but meteorological extremes such as passing storms tend to affect not only wind speed, but also air temperature, humidity, 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 yr, indicating that catchment runoff during storms is likely to have a minor influence. However, in lakes with short residence times (e.g., < 1 yr), storm-induced inflow can be at least as impactful as wind (Barbiero et al. 1999;Reichwaldt and Ghadouani 2012;De Eyto et al. 2016).

Conclusion
The primary result of our research was that high wind speeds (⪆ 10 m s À1 ) always had more decreasing effects on phytoplankton concentration than moderate wind speeds (≈ 5-10 m s À1 ), although the magnitude of the effect also depended on the level of incoming radiation and antecedent surface water temperature and hypolimnetic nutrients. The effect of wind decreased markedly when the MLD was about 8 m or deeper. Higher incoming radiation and hypolimnetic nutrient concentrations promoted increases in Chl a concentrations after wind events, whereas increases in surface temperature had a decreasing effect. These outcomes confirmed the conflicting effects of wind events on light and nutrient limitation of growth (Diehl et al. 2002;Kasprzak et al. 2017;Stockwell et al. 2020), and provide a mechanistic framework to better understand under what conditions wind events tend to either increase or decrease phytoplankton concentration. A simulation forced by a future climate scenario suggested that the response of phytoplankton to wind events did not strongly change with warming air temperatures despite earlier onset of stratification and a higher summer Chl a concentration.
Increased understanding of the drivers of wind impacts on lakes can help short-term forecasting and in some cases may be used to inform lake or reservoir management. In addition, this facilitates assessment of how atmospheric trends will affect lakes, specifically those caused by climate change. Future trends in extreme wind events include changes in intensity, duration, and frequency, and 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 wind events on lake phytoplankton could further our understanding of the global impact of extreme weather events on lake ecosystems.

Data availability statement
The GOTM code for the lake-branch is publicly available at https://github.com/gotm-model/code/tree/lake (last accessed: 29 July 2021). The code for the Selmaprotbas model can be found at https://github.com/jorritmesman/selmaprotbas (doi: 10.5281/zenodo.5094298), including instructions on how to couple it to GOTM. The model setup, including the configuration and input files, can be found at https://github.com/ jorritmesman/Erken_GOTM_SP_setup. Gap-filled hourly meteorological data for Lake Erken until 2018 have been published by Moras et al. (2019). Lake data until 2016 have been made available on https://www.ieg.uu.se/erken-laboratory/lakemonitoring-programme/ (last accessed: 29 July 2021). For the most recent data, contact Don Pierson (don.pierson@ebc.uu. se). The HadGEM2-ES model output for the EURO-CORDEX experiment was downloaded from the ESGF node of the German Climate Computing Centre (DKRZ).