Thermal soaring over the North Sea and implications for wind farm interactions

Seabirds use several flight modes at sea, including thermal soaring, in which thermal uplift is used to gain altitude and save energy. An increase in flight altitude may have consequences for wind farm interactions if it results in birds spending more time within the rotor-swept zone (RSZ). To understand conditions under which thermal soaring occurs and potential implications for wind farm interactions, we investigated thermal soaring in relation to atmospheric conditions in June and July at 2 study areas in the North Sea, west and north of the Dutch coast. We developed algorithms that identified thermal soaring in GPS tracks of lesser black-backed gulls Larus fuscus and radar tracks of seabirds. By combining species-specific 3-dimensional information on flight behaviour from bio-logging with the continuous spatiotemporal coverage of radar positioned at wind parks, we obtained a more comprehensive overview of thermal soaring at sea than either method would obtain alone. Our results showed that birds flew at higher altitudes during thermal soaring than non-soaring flight, increasing the proportion of flight time within the RSZ. Thermal soaring occurred inside offshore wind farms to a similar degree as outside. Thermal soaring was positively correlated with positive temperature differences (ΔT) between sea surface and air (at 2 m above sea level), and north and north-westerly winds. We show that the probability of thermal soaring over the North Sea, inside and outside wind farms, increases with larger temperature differences, resulting in increased time spent within the RSZ and an increased collision risk for seabirds.


Thermal soaring over the North Sea and implications for wind farm interactions
Section S1.Species distribution in study area The ESAS 5.0 (European Seabirds at Sea) database is a ship survey database (Reid & Camphuysen 1998) that incorporates all observational data made through a standardized protocol (Camphuysen et al. 2004) and is maintained by the NWO-NIOZ Royal Netherlands Institute for Sea Research.The database was queried on 4 January 2022 for all bird species observed in flight in the months of June and July in the area near Luchterduinen radar (Table S1) and Gemini radar (Table S2).For Luchterduinen, data was sampled between 52.15-52.65°N,3.5-4.5°E,which describes a rectangular area around the radar (52.427827°N, 4.185345°E).For Gemini, data was sampled between 53.75-54.25°N,5-7°E, which describes a rectangular area around the radar (54.036983°N, 6.041655°E).For readability, all species with an observation rate below 0.1% have been omitted from the table.
Gemini has 657 km² surveyed (2.237 km steamed on effort), with data collected between 1987 and 2012.Luchterduinen has 1607 km² surveyed (5354 km steamed on effort), with data originating from the same period.Figure S1 shows the range and altitude from the radar at which the probability of detection for an object of 1 standard avian target (SAT, Figure S1a) and 0.125 SAT (Figure S1b) by the horizontal S-band antenna of the RobinRadar 3D Fixed system is > 80%.A SAT is a theoretical object used as a standard for evaluating the performance of avian radar systems and approximates the physical features of a carrion crow (Corvus corone) with a radar cross section (RCS) of -16 dB m 2 and a mass of 500 g.A 0.125 SAT object correlates to a RCS of -25 dB m 2 and a mass of 62.5 g, the size of a song thrush (Turdus philomelos).The dynamic filter activity in each radar image directly affected detection probability of birds by increasing the threshold at which objects are detected.The filter is always active, ranging from 0 (no filtering) to 1 (complete filtering).The number of bird observations per hour is negatively related to the filter activity (black dots, Figure S2) and could affect the outcome of modelling efforts through inclusion of unreliable observation hours.Firstly, the estimated effect of the hourly averaged dynamic filtering on bird observations was modelled through Generalized Additive Modelling (GAM).Hourly bird count was used as the dependent variable, with hourly average filter activity as predictor (thin plate regression spline smoother, k = 5 to prevent overfitting), and assuming a Gaussian distribution of the model error.The estimated effect showed a decrease in bird count around filter activity = 0.1 (blue line, Figure S2), which levelled out around filter activity = 0.35.The filter activity at which this levelling out process starts was taken as the threshold above which hourly bird counts were considered highly affected.This value was found by calculating the second derivative over the estimated effect (red line, Figure S2) and finding the filter activity at the maximum of the curve (red dashed line, Figure S2).This threshold was 0.327 for Luchterduinen radar and 0.311 for Gemini radar; observation hours in which the average filter activity was higher than the threshold were excluded from further analysis.

Figure S2
The relation between hourly bird counts and radar filter activity in Luchterduinen radar.Hourly bird counts (left y-axis) at different filter activities (x-axis).Black dots depict all observation hours available.The blue line shows the estimated effect of average hourly filter activity on hourly bird counts (GAM).The red line shows the second derivative of this estimated effect (right y-axis), with the dashed red line depicting the value at which the second derivative was at its maximum, which was used to find the filter activity threshold for data exclusion.

Section S5. Overview of synoptic conditions in North Sea areas surrounding Luchterduinen and Gemini wind farm in support of thermal soaring
Soaring has been observed in both radar and GPS data on some specific days during the summer of 2019 and 2020.In the period 20-07-2020 -22-07-2020 high a peak in thermal soaring activity occurred in the west and north study areas around the Dutch coast (Figure 3, main text).To explore the synoptic conditions preceding, during and after this event (18-07-2020 -23-07-2020) we examined synoptic weather charts showing surface pressure and synoptic weather patterns for Europe at 6-hour temporal resolution, downloaded from the KNMI data center (https://www.knmi.nl/nederland-nu/klimatologie/daggegevens/weerkaarten),and the time series of sea surface temperature, air temperature (2 m) and wind speed (10 m) extracted from ERA5 for the two locations.Synoptic weather charts are presented in Figure S3 (taken at 12:00 UTC), and all weather charts for the period are presented in supplementary video S1.
On 18-07-2020 a stationary front is present over the North Sea.By the end of 19-07-2020 a cold front has formed along the Dutch coastline.The cold front leaves behind a trough of low atmospheric pressure at the end of the day and into the following days (solid blue lines), which brings cold maritime wet air from the northwest.The presence of such cold air creates suitable conditions for the development of thermals near the Dutch coast; the air temperature drops and remains below the sea surface temperature for several days (Figure S4.a-b), which drives a flux of sensible heat from the sea to the atmosphere (Markowski & Richardson 2010).Note that during these days the wind speed is not particularly high (Figure S4.c), which provides beneficial conditions for thermal soaring as it reduces the likelihood of thermal disturbance and break up.

Section S6. Remaining temporal autocorrelation in the logistic regression model residuals
Temporal autocorrelation in the model residuals of the proportion of thermal soaring per hour as function of temperature difference, and wind u-and v-components was addressed through a first order autoregressive covariance structure.However, this covariance structure could not completely resolve temporal autocorrelation in the models for the west area radar data.More complex covariance structures (second-/third order autoregressive and moving average autoregressive) did not improve the residuals further.Below is the remaining temporal autocorrelation plot for the west area radar logistic regression models.

Figure S1
Figure S1 Probability of detection by the horizontal S-band antenna of the RobinRadar 3D-fix at different ground range and altitude from the radar position (0,0 on the axes).(a) The area shows the range/altitude combinations at which the probability of detection > 80%.Probability of detecting a target of size 1 SAT.(b) Probability of detecting a target of size 0.125 SAT.Colour scale shows theoretical detection probability ranging from 100% (purple) to 80% (blue).Figures provided by RobinRadar

Section S4 .
Accounting for detection bias caused by dynamic radar filtering

Figure S5
Figure S5Autocorrelation function for the logistic regression model of proportion of circling tracks as a function of temperature difference between sea surface and air in the west area radar data.The x-axis shows the lag in hours from the observation, the y-axis shows the autocorrelation function.Autocorrelation is present for lag 1-6 and lag 20.

Figure S6
Figure S6Autocorrelation function for the logistic regression model of proportion of circling tracks as a function of temperature difference between sea surface and air in the west area radar data.The x-axis shows the lag in hours from the observation, the y-axis shows the autocorrelation function.Autocorrelation is present for lag 1-27, with severe auto-correlation for lag 1-9.

Table S1
Number of birds observed in flight and percentage of total per species in the months of June and July from 1987 -2012 in the area near wind farm Luchterduinen between 52.15-52.65°N,3.5-4.5°E.Data was queried on 4 January 2022.For readability, all species with an observation rate below 0.1% have been omitted

Overview of UvA-BiTs logger data per individual per year Table S3
Summary table of data collected for each tagged lesser black-backed gull per year.The same individual number over a different year indicates a returned individual.Total time indicates the sum of all time intervals for the recorded GPS measurements.