UvA-DARE (Digital Academic Repository) Built up areas in a wet landscape are stepping stones for soaring flight in a seabird

in order to reduce their transport costs, even in landscapes where thermal uplift isn't prevalent. We explore thermal soaring over land in lesser black-backed gulls by using high-resolution GPS tracking to characterise individual instances of thermal soaring and detailed energy exchange modelling to map the thermal landscape which gulls experience. We determinethatlesserblack-backedgullsareregularlyabletoundertakethermalsoaring,eveninawettemperateland- scape below sea level. By examining the relationship between lesser black-backed gull ﬂ ight, thermal uplift and land use,wedeterminethatbuiltupareas,particularlytownsandcities,providethermaluplifthotspotswhichlesserblack- backed gulls preferentially make use of, resulting in more opportunities for energy saving ﬂ ight through thermal soaring.


H I G H L I G H T S G R A P H I C A L A B S T R A C T
• High-resolution GPS reveals thermal soaring over land in lesser black-backed gulls. • We mapped bird's thermal landscape using a detailed energy exchange model. • Gulls soar to altitudes of hundreds of metres and then glide for several kilometres. • Built up areas create flyways via thermal soaring stepping stones in a wet landscape. • Human-engineering of the landscape impacts the aerial habitats of gulls.

A B S T R A C T A R T I C L E I N F O Editor: Rafael Mateo Soria
Keywords: Energy landscape Gulls Flight Thermal convection Bio-logging Land use The energy exchange between the Earth's surface and atmosphere results in a highly dynamic habitat through which birds move. Thermal uplift is an atmospheric feature which many birds are able to exploit in order to save energy in flight, but which is governed by complex surface-atmosphere interactions. In mosaic landscapes consisting of multiple land uses, the spatial distribution of thermal uplift is expected to be heterogenous and birds may use the landscape selectively to maximise flight over areas where thermal soaring opportunities are best. Flight generalists such as the lesser black-backed gull, Larus fuscus, are expected to be less reliant on thermal uplift than obligate soaring birds. Nevertheless, gulls may select flight behaviours and routes in response to or in anticipation of thermal uplift in order to reduce their transport costs, even in landscapes where thermal uplift isn't prevalent. We explore thermal soaring over land in lesser black-backed gulls by using high-resolution GPS tracking to characterise individual instances of thermal soaring and detailed energy exchange modelling to map the thermal landscape which gulls experience. We determine that lesser black-backed gulls are regularly able to undertake thermal soaring, even in a wet temperate landscape below sea level. By examining the relationship between lesser black-backed gull flight, thermal uplift and land use, we determine that built up areas, particularly towns and cities, provide thermal uplift hotspots which lesser blackbacked gulls preferentially make use of, resulting in more opportunities for energy saving flight through thermal soaring.

Introduction
The distinctive circling flights of birds soaring on thermals, whereby birds gain height without flapping their wings before undertaking sustained glides, is a long observed movement phenomenon that has contributed to Science of the Total Environment 852 (2022) 157879 our understanding of avian movement at a range of scales, from migration (Leshem and Yom-Tov, 1996) to fine-scale flight control (Williams et al., 2018). Over land, the thermal convection which supports this type of flight is driven by surface-atmosphere interactions in the atmospheric boundary layer. Birds have to navigate this complex aerial environment so the ability to detect and utilise thermals offers a highly useful opportunity to save energy for many soaring birds (Hedenstrom, 1993). Therefore, surfaceatmosphere interactions, especially those linked to thermal lift, may influence bird flight decisions and flight costs (Bohrer et al., 2012;Mandel et al., 2011;Santos et al., 2021).
When the Earth's surface absorbs solar energy (radiation), the energy is partitioned towards heating the soil (ground flux), evaporation (latent heat flux) and heating the air (sensible heat flux). As air is heated, air density decreases and causes the air to rise, creating thermal lift. If less radiation is taken up by ground flux and latent heat flux, more energy is available for sensible heat flux, resulting in greater thermal production. Thermal lift is thus strongly dependent on atmospheric conditions (mainly solar radiation and air temperature) and surface properties (mainly land use, soil type and soil water content) that determine the energy partitioning over the three components ground flux, latent heat flux, and sensible heat flux (Wilson et al., 2000). Where and when absorbed radiation is higher, soil water content is lower, and land cover generates less evaporation, we expect thermal production to be higher.
The degree to which birds utilise thermal lift and the characteristics of thermal soaring flight is dependent on multiple interacting internal and external factors, relating to bird morphology, meteorological conditions, and the behavioural response of birds to their environment (Pennycuick, 2008;Williams et al., 2018). Bird morphology, particularly relating to wing dimension and body mass, has been used to understand differences in thermal soaring characteristics between species (Pennycuick, 1983;Spaar, 1997). By then integrating meteorological information, either using proxies of thermal strength such as boundary layer depth, or by modelling vertical updraft, it becomes possible to study the effect of atmospheric conditions on thermal soaring across morphologically different species (Shamoun-Baranes et al., 2003), or to gain insight into behavioural variation in thermal soaring among morphologically similar species (Bohrer et al., 2012). Behavioural variation during thermal soaring can be probed further to investigate how birds make decisions in relation to time and energy; from adjusting air speeds during gliding (Katzner et al., 2015;Vansteelant et al., 2017;Nourani et al., 2021), to timing departures from thermals (Harel and Nathan, 2018;Shepard et al., 2011), to route formation and use of social information (van Loon et al., 2011). Overall, soaring birds have been shown to adjust many aspects of their flight behaviour in relation to their morphology, internal motivation and various proxies of thermal uplift.
Not only are the interactions which govern thermal soaring in birds' complex, but so too are the interactions which govern the distribution of thermals. The distribution of thermal uplift across a landscape forms part of the energy landscape , as soaring birds can obtain a lower cost of transport by optimising soaring flight in relation to the thermal landscape, which is a particularly important strategy in areas where thermal availability is patchy or sparse. Many animals have been shown to establish movement corridors through landscapes which reduce their transport costs (Panzacchi et al., 2016), including in relation to atmospheric characteristics (Brandes and Ombalski, 2004;Sage et al., 2019). This principle could also be applicable to the thermal landscape, whereby certain areas will provide hotspots of thermal uplift which birds can utilise for circling and then glide in-between. Selection of routes based on thermal uplift has been observed on a large scale in obligate soaring birds during migration (Bohrer et al., 2013;Scacco et al., 2019) but not at finer scales or during daily foraging movements, due in part to the complexities of modelling thermal uplift on small spatial scales and across highly spatio-temporally variable landscapes. However, fine-scale changes in the spatio-temporal thermal landscape may cumulatively impact flight costs and energetic trade-off decisions. The ways in which birds make these kinds of energy saving cost-benefit decisions will also depend on species morphology, behavioural specialisations and internal motivation; an obligate soaring bird may depend more on thermal uplift hotspots in relation to thermals than a flight generalist that can more readily switch to flapping flight when alternatives are not available.
Lesser black-backed gulls Larus fuscus are flight generalists with the capacity to utilise flapping as well as soaring flight (Shamoun- . They are not generally associated with thermal soaring over land, as they do not have the typical morphology or flight characteristics of obligate soaring species. However, lesser black-backed gulls are known to undertake soaring flight in response to conditions associated with thermal convection and orographic lift Shepard et al., 2016), amounting to around 30 % of their flight time (Sage et al., 2019;Shamoun-Baranes et al., 2016). Previous attempts to directly measure the thermal uplift experienced by lesser black-backed gulls have been limited by medium-range weather models, which cannot resolve complex micro-scale spatio-temporal distributions of thermal uplift . The degree to which lesser black-backed gulls depend upon thermal convection is therefore not known; their movement is not likely to be halted by an absence of thermal uplift, but atmospheric conditions may impact their flight mode and thereby their flight costs (Sage et al., 2019;Shamoun-Baranes et al., 2016;. Throughout the breeding season, lesser black-backed gulls are central place foragers making daily trips to foraging areas up to tens of kilometres from their breeding colony Garthe, 1997;Thaxter et al., 2013), with individuals exhibiting a range of foraging strategies in diverse marine, intertidal and (increasingly over recent years), terrestrial habitats Corman et al., 2016;Isaksson et al., 2016;Spelt et al., 2019;Tyson et al., 2015). During repeated terrestrial foraging trips lesser black-backed gulls have been shown to anticipate and adjust their flight behaviour and routes in response to orographic uplift (Sage et al., 2019). Lesser black-backed gulls may therefore also be able to perform similar behavioural adjustments in response to thermal uplift, although the thermal uplift likely occurs with different spatio-temporal dynamics than orographic uplift which may necessitate different behavioural strategies. As lesser black-backed gulls repeatedly transverse and utilise diverse habitats during their daily flights, they provide a novel system in which to study the phenomenon of thermal soaring in a mosaic landscape. Thereby, physical, environmental, and behavioural limitations of avian flight behaviour in relation to thermal soaring can be examined at a fine spatio-temporal scale, based on the scale at which they likely respond to variation in the thermal landscape.
The aim of this paper is to provide a comprehensive insight into thermal soaring in lesser black-backed gulls, which both describes the characteristics of thermal soaring in lesser black-backed gulls and explores the role of a mosaic landscape in shaping soaring behaviour. We break our study down into three sub-questions: First, what are the characteristics of thermal soaring in lesser black-backed gulls and under what uplift conditions do we observe thermal soaring? For this question we identify explicit bouts of thermal soaring from high-resolution GPS and accelerometer data and model the thermal uplift encountered during thermal soaring. Second, how does thermal uplift across a mosaic landscape influence flight mode and flight routes? Here we model thermal uplift across centrally-placed trips and a random flight path model to determine whether gulls select routes and adapt their flight mode in response to available thermal uplift. Third, to what extent does the landscape support thermal soaring? Here we model landscape wide thermal soaring potential throughout the study period, deemed thermal soaring suitability. We examine whether gulls select areas with high soaring suitability and examine how the relationship between land use and thermal soaring suitability influences route choices.

Ethics statement
Field work to fit GPS-loggers to lesser black-backed gulls was conducted with approval from and in accordance with the Dutch ethics committee on animal experiments (DEC) of the Royal Netherlands Academy of Arts and Sciences (KNAW). Permission to work in the Kelderhuispolder gull colony was granted by the Staatsbosbeheer.

Study area and study period
The lesser black-backed gulls monitored for this study breed on the island of Texel, Netherlands, on an interface between the North Sea, Wadden Sea, and mainland North-Holland. The time-frame of this study, 1 May-31 July 2018, coincides with the breeding season, during which time gulls make trips lasting several hours from their colony to terrestrial and marine locations to feed or rest Tyson et al., 2015). This study focuses on the land area of North-Holland used by gulls for terrestrial foraging, a flat landscape predominantly at or below sea level characterised by a mosaic landscape of arable and urban land, varying soil types, and temperate climate. It is a typical Dutch engineered landscape with ground water tables that are artificially controlled at 20-40 cm below the surface for pastures and 40-100 cm for arable land.

GPS data
The tracking data used in this study was retrieved from a long-term tracking dataset of lesser black-backed gulls breeding on Texel, tracked using the UvA Bird Tracking System since 2008 . Data consisted initially of 22 individuals, before being filtered according to data interval resolution and trip definitions (see Sections 2.5 and 2.6). GPS measurement intervals ranged between 3 and 600 s, with 300 s being the default interval setting, and a high resolution (3 s) interval setting taking over when batteries are fully recharged. Following every GPS measurement, tri-axial acceleration was measured at a frequency of 20 Hz for 1 s and used for behavioural classification. Behaviour was classified using a random forest classifier described in Shamoun- . For the purpose of this study the class 'soaring' was used to classify thermal soaring flight (see Section 2.5) and all flight behaviours were combined to identify flight (Section 2.6).

Thermal uplift model
Thermal uplift was modelled using a simplified estimate of the convective velocity scale w* (m s −1 ) (Stull, 1988;Bohrer et al., 2012): where g is the gravitational acceleration, z is the boundary layer height (m), T is the potential temperature (see Supplementary material S.1), c p is specific heat capacity at constant pressure, and H is the surface sensible heat flux (J s −1 m 2 ) (Deardorff and Willis, 1974;Holtslag and Boville, 1993). A dry air assumption was made, neglecting the contribution of water vapor flux to w*, which is presumed to be small in the convective boundary layer over land. Convective velocity scale w* is a coefficient that scales with the magnitude of turbulent vertical transport in the planetary boundary layer (Deardorff, 1970) and indicates where in horizontal space vertical transport is likely to be larger. The mosaic landscape of North-Holland is expected to contain high w* gradients due to discrete variations in land use, generating large differences in H, which is a critical driver of the spatio-temporal variation in w*. This difference in H arises from variation in attributes of changing land uses, such as in albedo and surface roughness, which in combination with the availability of water dictate the amount of energy that is converted to H. Therefore, in order to capture the spatio-temporal variation in w*, H was explicitly modelled using a set of models that solved the energy balance at the Earth's surface at a high spatio-temporal resolution (at a 5 m spatial and 1 h temporal scale).
The energy balance equation partitions radiation into latent heat, warming of the ground and sensible heat (Budyko, 1961) where H is the sensible heat flux, R n is the net radiation, λE is the latent heat produced by evaporation and G is the ground flux. In this research R n was retrieved from the ERA5 reanalysis data set (Hersbach et al., 2018(Hersbach et al., , 2020, leaving G and E to be modelled in order to supply H as a residual. In the following section the fundamental principles behind the calculation of G and λE in the model are described. A more elaborated description can be found in Supplementary material S.1.

Latent heat flux
Hydro-meteorological research has shown that evapotranspiration (the sum of direct water evaporation and the transpiration that evaporates through stomata of plants) can be predicted with a crop factor approach (de Bruin and Lablans, 1998). Such methods use a reference evapotranspiration over a well-watered grassland, multiplied by dimensionless crop factors to obtain the evapotranspiration for various types of land use. Crop factors can change throughout the growing season to account for the development of plants. The LGN2018 land use map (Rip and Hazeu, 2019) with a spatial resolution of 5 m was used and linked to dynamic dual crop factors (Allen et al., 1998) that allowed for the separation of plant transpiration and soil water evaporation, alongside the inclusion of effects of soil water deficit on transpiration by plants and evaporation on bare soil (see Supplementary material S.1. for details). Soil water deficits were taken from the NHI model (de Lange et al., 2014). The evapotranspiration (E in kg m −2 h −1 ) was multiplied by the specific heat of vaporization (λ in Jkg −1 ) to calculate the latent heat flux, the energy needed for evapotranspiration. The leaf area index (LAI) of land use classes was quantified with the LAI300 data set (Baret et al., 2016).

Ground flux
To estimate the ground flux G, the energy absorption by the soil and the heat transfer to greater soil depths was calculated with a dynamic model based on the diffusion and continuity equations for heat exchange (Koorevaar et al., 1983), with a vegetation layer, a soil surface layer and eight soil layers, down to a depth of 3 m (details are given in Supplementary material S.1). The soil stores energy when the radiation is high during the day, while reducing the sensible heat flux, and releases energy during the night (Panagos, 2006).

Sensible heat flux model data retrieval and model implementation
When the net radiation, the ground flux and the latent heat flux are known, the sensible heat flux can be calculated with Eq. (2). Confirmation of the calculated sensible heat flux was carried out by comparing the dynamic modelling results at a cloudless moment in June with the Surface Energy Balance Algorithm for Land (SEBAL) (Bastiaanssen et al., 1998), for which MODIS infrared satellite images are used (details can be found in Supplementary material S.2).

Characterising thermal soaring
The flight behaviour of lesser black-backed gulls in relation to thermal lift was examined on two temporal GPS measurement scales in conjunction with the research aims. High-resolution data was used to identify and characterise thermal soaring in lesser black-backed gulls in relation to thermal lift. Following this, a larger dataset of GPS measurements at a wider range of temporal resolutions comprising all daily terrestrial movements was considered in order to examine how thermal lift influences flight mode and route decisions.
In order to identify specific instances of thermal soaring, highresolution GPS and accelerometer data were selected (location interval < 10 s), as they provided high enough temporal frequency to identify thermal circling. All points were grouped into periods of uninterrupted high-resolution measurements, deemed 'bouts' (where uninterrupted was defined as temporal gaps between positions not exceeding 600 s, allowing for small data gaps). Within a bout, points were assigned 1 of 3 categories: circling, gliding, and other, i.e. neither circling nor gliding were identified. To allocate these classifications, an algorithm was applied to the tracks which calculated moving averages over 5 locations at a time, defining climb rate, turning angle, ground speed, and proportion of time spent soaring (where soaring flight points had been defined by the accelerometer behavioural classification, see Section 2.3). Where climb rate was >0, turning angle between subsequent measurements was >30°, ground speed >5 m s −1 , and proportion of time spent soaring was >0.6, a location was defined as circling. A soaring proportion > 0.6 ensures that the majority flight behaviour was soaring while allowing for incidental flapping, which was commonly observed during exploration of the gull GPS data. Where climb rate was ≤0, ground speed was >5 m s −1 and proportion of time spent soaring >0.6, a location was defined as gliding. Gliding points which were not preceded by circling points were removed, as they could not be conclusively connected to thermal soaring behaviour. Consecutive circling points were considered a circling bout, consecutive gliding points were considered a gliding bout, and a circling bout with its directly subsequent gliding bout was considered a thermal soaring bout.
Summary statistics were determined for every thermal soaring bout. The following metrics were calculated for circling: mean climb rate, total time spent circling, thermal entry altitude, thermal exit altitude, circling radius, proportion of boundary layer height reached. For gliding, total time spent gliding and glide ratio were calculated. A complete summary of thermal soaring metrics and details of how each of these variables were defined can be found in Supplementary material S.3.
To determine the thermal uplift associated with circular soaring, an estimate of w* was made for each circling bout. First the origin location of each thermal was estimated. Birds regularly enter thermals at altitudes above ground level; when thermals are tilted by horizontal winds (Treep et al., 2016), the area at the Earth's surface which triggers a thermal may not be the area over which circular soaring is observed. To backtrack the origin location of a thermal, a vector was drawn between the first and last circular soaring point within a circling bout and geometrically retraced to its origin. Convective velocity scale w* was then calculated for each circling bout, using Eq. (1) spatio-temporally annotated onto the thermal origin location (boundary layer height (m) and pressure measurements used to calculate potential temperature in Eq. (1), were retrieved from ERA5 reanalysis dataset, further details in Supplementary material S.1).
Initially, w* was calculated on the same spatial scale as provided by H, a 5 m grid size. However, thermal plumes are generated over a larger spatial footprint (Stull, 1988) of tens to hundreds of metres, therefore the expected area for which H will influence w* is larger. We additionally estimated the median radius of circling bouts by calculating a linear regression through circling GPS coordinates to determine the centre of a thermal, then calculating the median horizontal distance between circling points and the estimated centre. This gave a mean radius of 27.4 m (IQR: 22.0-36.1), indicating a minimum thermal diameter of at least twice this value (while in practice circling radius is a property of bird morphology and birds are unlikely to utilise the entire width of a thermal (Pennycuick, 2008). Therefore, H was aggregated to a lower spatial resolution by taking the mean of H over a larger grid size of 105 × 105 m and applied to Eq. 1 to estimate w* at a 105 × 105 m spatial and hourly temporal resolution. 105 m was used rather than 100 m to allow for an equal number of spatial cells on each side of the central cell containing the GPS location.
2.6. Flight behaviour in response to thermal uplift 2.6.1. Defining flight routes and destinations To examine how thermal uplift influences flight mode and route decisions, GPS and accelerometer data of all location intervals were included in determining daily movement trips. A comparative analysis was made with concurrent randomly generated trips, following Sage et al. (2019). At this scale it was not possible to identify individual instances of thermal circling or gliding. Instead, points were presumed to be associated with thermal soaring (either by circling or gliding) if their behaviour was classified as soaring and if the altitude exceeded 54 m above ground level. This altitude was chosen based on the result of calculating the median and interquartile range in minimum circling altitudes from Section 2.5. As a general measure of the efficacy of this classification for thermal soaring, a comparison was made to the high-resolution GPS data used in Section 2.5 by looking at the number of high-resolution data points at altitudes above 54 m with a soaring accelerometer classification. Above 54 m, 98 % of high-resolution GPS locations with a soaring accelerometer classification had also received a thermal soaring classification (either circling or gliding) and this altitude cut off was therefore deemed to be an effective measure for estimating thermal soaring in all the GPS locations.
Movement data was classified into trips, defined by an individual travelling further than 500 m from the colony (53.00 N, 4.72 E), staying away from the colony for at least 1 h and travelling a cumulative distance of at least 4 km before returning. Destinations of each trip were also identified and used as a constraint in the generation of random trips. Destinations within a trip were defined where the distance moved over a 15-minute interval either side of a GPS location was <100 m. A total of 176 trips were identified across 17 individuals throughout the study period, ranging from 1 to 24 trips per individual. For each trip, summary metrics were calculated: the total trip duration calculated as the sum of all GPS intervals, the total trip distance calculated as the sum of all distances between GPS intervals.

Generating random routes
A null model of lesser black-backed gull movement was generated to provide a comparison to the gull's true movements. A random route was generated for every real route, connecting the destinations along the route using a random trajectory model devised by Technitis et al. (2015). The resulting random route had the same trip duration and destinations as the real route, but trajectories could deviate from the true track. Speed was kept constant within each random route for ease of calculation. A more in-depth description of the model can be found in Sage et al. (2019).

Modelling thermal uplift for real and random routes
For the purpose of measuring flight response to thermal uplift, only locations designated as flight in the accelerometer classification were considered part of a real route, and every equivalent location from the random trajectory model was considered the random route. For all real and random routes, a value of w* was modelled, using the same annotation method as in Section 2.5. Since the temporal resolution was not always high enough to identify circling bouts, the spatial origin of a thermal could not be determined, so H and z i were calculated for the area directly underneath each GPS location. As in Section 2.5, H was spatially aggregated such that w* was calculated over a 105 by 105 m spatial and hourly temporal resolution.
In order to determine whether gulls adjusted their flight behaviour in relation to thermal lift in real time, the w* experienced during real routes was compared to that of random routes and to that of thermal soaring flight (as defined at the start of Section 2.6). Distributions of w* across these three categories were presented as relative probability density distributions, weighted by GPS measurement time interval, and statistically compared via a two sample Kolmogorov-Smirnov (KS) test. Additionally, to examine the w* experienced during thermal soaring compared to the w* experienced specifically during circling, the distribution of w* in the high-resolution circling bouts was presented as a relative probability distribution and compared to thermal soaring flight via a KS test.

Environmental influence on opportunities for thermal soaring
In order to examine how opportunities for thermal soaring are distributed in space and in relation to land use, a spatial heatmap of thermal soaring opportunities was generated based on availability of suitable thermal uplift throughout the study period, termed the thermal soaring suitability landscape. Thermal soaring suitability was calculated as the proportion of time each spatial cell generated sufficient uplift for soaring. To define sufficient thermal uplift for soaring, a thermal soaring threshold was calculated using the relative probability density distributions of w* in real routes and thermal soaring flight detailed in Section 2.6. Above a given value of w* the relative probability of thermal soaring increased beyond that of the real routes. This meant the probability of drawing a thermal soaring location from ranges above the given value of w* was higher than the probability of drawing any type of flight behaviour from within the same range. The given value of w* at which this occurred was deemed the thermal soaring threshold and indicated the point beyond which lesser black-backed gulls are more likely to use thermal soaring modes of flight than any nonspecific flight mode.
To generate the thermal soaring suitability landscape, w* was calculated for all cells in the study area based upon Eqs. (1)-(2) on an hourly basis across the entire study period. Every value in space and time was assigned a 1 or 0 according to whether it exceeded the thermal soaring threshold, then for each cell a proportion value was calculated, i.e. the proportion of all time across the study period for which that cell exceeded the thermal soaring threshold. The thermal soaring suitability landscape indicated the expected opportunities for thermal soaring in horizontal space. Therefore, in order to determine whether gulls displayed a preference for areas with higher thermal uplift suitability, in their flight routes and during thermal soaring, all real and random routes were spatially annotated with thermal soaring suitability, aggregated by taking the mean of a 105 × 105 m cell centred on the GPS location (as in Section 2.5). The distributions of thermal soaring suitability for real and random routes, as well as thermal soaring locations were presented as relative probability distributions weighted by GPS measurement time interval. These distributions were statistically compared via a two-sample KS test, as in Section 2.6.
In order to relate opportunities for thermal soaring to land use, all real and random routes were also annotated with the LGN2018 land use map to assess the types of land use utilised by gulls and their reliability for producing thermal uplift. In order to include all land uses which may influence bird behaviour, the proportion of all land types in the 105 × 105 m grid around each GPS point was calculated. The thermal soaring potential associated with each land use was then weighted by the proportion of the grid for each GPS point which it covered, as well as the time interval of the GPS point to which it was associated.
All data preparation, modelling, and analysis were carried out in ArcGIS Pro, Rstudio version 4.0.2, Matlab 2017a and Matlab R2019b.

Characterising thermal soaring
Across 144 bouts of high-resolution GPS data, 873 instances of thermal circling were identified from 6 individuals. An example of a thermal soaring bout is shown in Fig. 1. Key metrics of the thermal soaring observed from the high-resolution bouts, which can be then compared to equivalent metrics in other species, are provided in Table 1

Flight behaviour in response to thermal uplift
Across 176 trips identified in 17 individuals, median trip duration was 3.7 (IQR: 2.3-6.0) hours, corresponding to a median total point-to-point trip distance of 78 (IQR: 49-143) km. Based on the accelerometer data, 38 % of time in flight was spent soaring, and 21 % of flight time was spent soaring at altitudes above 54 m, deemed to be thermal soaring flight.
A comparison of the w* experienced across all flight routes, the associated random routes, all thermal soaring points, and high-resolution circling bouts are shown in Fig. 2. Overall, the distribution of w* in the real and random routes displayed a similar bimodal distribution, with the right-hand distribution being associated with positive w*, or rising air. The real and random distributions were statistically similar (KS test D = 0.063, p = 0.27). Thermal soaring had a much lower probability of occurrence at negative w* and a much greater probability of occurrence at higher positive values of w*. The distribution of w* during thermal soaring was statistically distinct from both real and random routes (KS test between thermal soaring and random routes: D = 0.18, p < 0.0001, KS test between thermal soaring and real flight: D = 0.17, p < 0.0001). The distribution of w* during thermal soaring compared to high-resolution circling was also statistically distinct (D = 0.15, p < 0.0001), with circling predominantly occurring at higher values of w*. Above w* values of 0.97 m s −1 the probability of soaring was greater than the probability of all flight, so this value was deemed the thermal soaring threshold. The probability drawn from the distribution of thermal soaring when w* ≥ 0.97 m s −1 was 0.81, whereas the probability drawn from the distribution of flight routes in this same range was 0.51, and 0.45 for random routes.

Environmental influence on opportunities for thermal soaring
Based on Section 3.2, the thermal soaring threshold above which lesser black-backed gulls were estimated to show greater likelihood of thermal soaring was taken to be 0.97 m s −1 and was therefore used to model the thermal soaring suitability landscape. The landscape generated is shown in Fig. 3, and illustrates the regularity with which thermal soaring opportunities occurred across the study area. The thermal soaring suitability (i.e. the proportion of time per cell for which the thermal soaring threshold was exceeded) ranged from 0 to 0.465. The spatial distribution of thermal soaring suitability formed a mosaic landscape, where a network of urban areas and roads generated a higher thermal soaring suitability (Fig. 3). Fig. 3a shows the spatial distribution of high-resolution thermal circling locations (backtracked to an estimated thermal origin point) as described in Section 2.5 superimposed on the thermal soaring suitability landscape. Circling locations were spread across the network of flight routes (Fig. 3b), but showed clustering in certain areas, particularly urban areas. All the measured flight routes of lesser black-backed gulls over the thermal soaring suitability landscape (Fig. 3b) reveal several areas of high thermal suitability where paths converge, with some convergence areas overlapping with the clustering of circling locations seen in Fig. 3a.
Comparisons of thermal soaring suitability in real, random, and soaring routes revealed common peaks in thermal soaring suitability as well as some notable differences (Fig. 4a). All distributions were statistically distinct from each other (KS test between real routes and random points: D = 0.17, p < 0.0001, KS test between random routes and thermal soaring flight D = 0.25, p < 000.1, KS test between real flight and thermal soaring flight D = 0.18, p < 0.0001). All distributions peaked at values of thermal soaring suitability around 0.27, with the greatest peak being observed in the random routes. This indicates that the landscape contains a significant amount of area which is producing thermal uplift above the thermal soaring threshold 27 % of the time. A second notable peak was observed towards the highest values of thermal soaring suitability, around 0.45. This peak was larger in real routes than random, and even larger in thermal soaring points. Above a thermal soaring suitability value of 0.3, the likelihood of real routes being encountered over random routes is always greater. Above this value, the probability of a sample being drawn from the distribution of thermal soaring was 0.58, for real routes was 0.56, and for random routes was 0.44. This also means the majority of flight movements and of thermal soaring occurs in areas exceeding the thermal soaring threshold at least 30 % of the time.
The relationship between land use, thermal soaring suitability and gull movements is shown in Fig. 4b-d. The distribution of thermal soaring suitability for each land use class is shown for random routes (Fig. 4b), real flight (Fig. 4c) and thermal soaring flight (Fig. 4d). The peak at thermal soaring suitability around 0.27 is linked to agricultural land types, particularly agricultural grass, bulbs, potatoes and other agricultural uses. The second peak at thermal soaring suitability around 0.45 is predominantly associated with urban areas, particularly urban grassland, buildings, roads and railways. The difference between the thermal soaring suitability experienced during real routes vs random routes can, according to Fig. 4b and c, be attributed to a shift in land use. Gulls are less likely to fly over agricultural areas associated with the peak at 0.27, and more likely to fly over   urban areas associated with the peak at 0.45. This difference is also observed between random routes and thermal soaring flight, with the difference associated with urban areas being even higher. In Fig. 5 we visualise the changes in altitude and w* experienced over the course of one long trip in which thermal soaring was abundant. Four specific thermal soaring bouts are illustrated in relation to w* in the surrounding landscape modelled at same hour in which the soaring bouts occurred.

Discussion
Lesser black-backed gulls are frequently able to undertake thermal soaring, with clear links between their flight behaviour and the atmospheric conditions they experience. We estimate that up to 21 % of their flight time is spent in thermal soaring, which presents a considerable opportunity for energy savings in flight. The North-Holland landscape is a patchy mosaic landscape predominantly reclaimed from former sea; it is very flat with high ground water levels. Nevertheless, it regularly generates atmospheric conditions suitable for thermal soaring in lesser black-backed gulls. The thermal soaring suitability of this landscape can be ascribed in part to qualities arising from human engineering of the land, particularly the presence of urban areas that serve as thermal uplift hotspots.

Characteristics of thermal soaring in lesser black-backed gulls
By characterising individual instances of thermal soaring in lesser blackbacked gulls we gain an understanding of their flight capabilities and strategies, particularly in comparison to other species more associated with thermal soaring in a range of environmental conditions. Despite not being associated with terrestrial obligate soaring birds in terms of flight morphology or strategy, in many respects lesser black-backed gulls proved quite comparable to them in climb rate and glide ratio. Climb rate during circling was 1.16 m s −1 , which compares well with studies of obligate soaring birds (Harel et al., 2016;Shepard et al., 2011;Spaar and Bruderer, 1996). Glide ratio in lesser black-backed gulls was 15.2, again comparable to large obligate soaring birds  and also to aerodynamic theory (Pennycuick, 2008), which predicts a minimum glide ratio of 15.6 and best glide ratio of 19.4 (for a reference lesser black-backed gull with mass = 0.806 kg, wing area = 0.191 m 2 , wing-span = 1.34 m, aspect ratio = 9.42).
Where characteristics of thermal soaring flight in lesser black-backed gulls deviate most from other species relates mostly to environmental differences, e.g. the scale of thermals encountered. Soaring behaviour in lesser black-backed gulls has previously been linked to boundary layer height , which generally dictates thermal height. Typical boundary layer heights reach around 400-700 m in our wet landscape, while in other studies of thermal soaring over land boundary layer heights or thermal depths of 1000-2000 m are normal (Leshem and Yom-Tov, 1996;Shamoun-Baranes et al., 2003;Shannon et al., 2002;Williams et al., 2020). Lesser black-backed gulls also use a slightly lower proportion of the boundary layer, 0.48, than has been reported for some larger soaring species (Shamoun-Baranes et al., 2003). This could indicate that either morphology is limiting the utilisation of the top of boundary layer where thermal convection typically becomes weaker, or that a behavioural trade-off exists which gives preference to leaving thermals earlier where lift rates are higher (Harel and Nathan, 2018;Shepard et al., 2011). These behavioural trade-offs are dependent on an individual's internal drivers, which vary throughout the year and over different movement scales. For example, birds may alter behavioural trade-offs when leaving thermals depending on whether they are undertaking daily foraging movements (such as in this research), or long migratory flights, where time and energy are under Graph shows the change in altitude with increasing distance from colony, coloured according to the real time w* experienced (annotated to each GPS point for 105 m by 105 m area directly below each point). Grey rectangles indicate where circling was identified, also indicated by an increase in altitude, while decreases in altitude generally indicate gliding. Letter annotations link to 4 circling events shown on the below maps, where routes are plotted over a map of the local landscape w* at the time of the bout, given in local time (UTC + 2). different constraints. During shorter trips the relative benefit of saving energy during trips may be lower.
Overall, thermal soaring performance in lesser black-backed gulls is in many ways comparable to obligate soaring species, but dependence upon thermals still varies in relation to morphology, environment and internal drivers. The need for lesser black-backed gulls to moderate thermal soaring behaviour in anticipation of predictable thermals may be lower than in obligate soaring birds, since in the absence of thermals it is relatively less costly for gulls to switch to flapping flight.

Influence of thermal uplift on flight mode and route choice
Gulls clearly tend to use soaring flight when thermal uplift experienced is strong (Fig. 2); as thermal uplift increases, the probability of thermal soaring also increases relative to the overall distribution of thermal uplift available. The probability of circling is nearly absent at w* values below the theoretical minimum sink rate for a lesser black-backed gull, 0.469 m s −1 (Pennycuick, 2008), supporting the expectation that maintained soaring would not be possible in updrafts beneath the minimum sink rate. When w* ≥ 0.97 m s −1 the probability of thermal soaring is higher than for all flight, which provides an indication of the range of w* where soaring is most likely to take place. When the circling part of thermal soaring is isolated and examined in relation to w*, a stronger selection for thermal uplift is seen compared to thermal soaring as a whole (Fig. 2). This is expected, as circling flight is the component of thermal soaring which directly depends on thermal uplift while gliding can occur in the absence of thermal uplift, (although thermal uplift is still profitable as it reduces sink rate during gliding, see Spaar and Bruderer (1997)).
While we clearly see that thermal uplift influences flight mode, it is less clear whether thermal uplift influences route choice, as thermal uplift experienced during real routes is comparable to random routes (Fig. 2). This could be due to the way hotspots of thermal soaring suitability in a mosaic landscape drive the available uplift. Lesser black-backed gulls showed a greater likelihood of flying, and specifically soaring, over areas mostly associated with urban land uses, predominantly urban grass, buildings, roads and railways (Fig. 4). These areas also provide the suitable thermal uplift for soaring most regularly (Fig. 4), supporting the idea that gulls perceive which areas are more likely to provide suitable thermal uplift, which then influences their route choice and flight mode. These urban areas are distributed in a patchy way across an otherwise predominantly agricultural landscape (Fig. 3), creating hotspots of reliable thermal convection that birds can use as stepping stones to soar upwards before gliding onwards to another hotspot. This mechanism is reflected in the high densities of circling observed over thermal hotspots in the landscape (Fig. 3a). An example is provided in Fig. 5 to demonstrate the use of thermal uplift hotspots as stepping stones: an individual undertakes several thermal soaring bouts within one trip and generally encounters higher w* when climbing than when gliding (also see red lines in Fig. 2). If hotspots are regularly spaced throughout the landscape, a soaring bird can travel efficiently between them with minimal deviation to its route. North-Holland has regularly interspersed hotspots of thermal uplift suitability (Fig. 3), especially inland to the west, where many flight paths converge. Lesser black-backed gulls appear to avoid the more wet reclaimed land (polders) to the east where thermal uplift suitability is generally lower (Fig. 3b), compared to random routes (Supplementary material S.4).

Support for thermal soaring in a patchy wet landscape
Land use plays a key role in generating the mosaic thermal uplift landscape. The majority of North-Holland is dedicated to a range of agricultural land uses, where most absorbed energy is used for evapotranspiration and suitable thermal uplift is limited to around 27 % of the time (Fig. 4). Urban hotspots, with thermal suitability above 40-45 % of the time, are predictable resources for thermal soaring throughout the breeding season as there is limited water for evaporation and most energy is used to heat the air (Fig. 3a). Conversely, arable areas are influenced by different types of vegetation, which affect evapotranspiration through land cover characteristics and ground flux through changes to the root zone and soil moisture, causing complex, and thus difficult to predict, spatio-temporal dynamics in w* throughout the season (Supplementary material S.5). As a result, while arable areas are abundant in the landscape and sometimes generate strong thermal uplift, they do not provide reliable hotspots through which routes can be formed (also see Supplementary video S.6).
Potential for high thermal uplift is expected to be lower in flat or wet regions (Kerlinger, 1989;Stull, 1988), however, we now see that the landscape of North-Holland generates a surprising amount of thermal uplift that supports thermal soaring in lesser black-backed gulls, particularly when considering climate, land use, and water table in comparison to other studied environments (Kerlinger, 1989;Leshem and Yom-Tov, 1996;Mandel et al., 2011;Santos et al., 2017). The range of thermal uplift experienced by lesser black-backed gulls (Fig. 2) is comparable to that experienced by golden eagles, Aquila chrysaetos, and turkey vultures, Cathartes aura in the Eastern United States (Bohrer et al., 2012), despite significant differences in topography, land use, and weather conditions. Urban thermal hotspots are the main driver for these high levels of thermal uplift in North-Holland (Fig. 3) and were identified as a result of modelling thermal uplift at a very high spatial resolution, meaning at coarser scales the range and distribution of thermal uplift found in this study may not have been identified (Supplementary material S.3). The spatial scale used to examine thermal uplift in this study was tailored towards understanding how lesser black-backed gulls use thermals during circling behaviours in relation to fine scale landscape variation. It is only by increasing the spatial and temporal resolution over which we model thermal uplift that we are able to resolve these differences between landscape features and examine how birds utilise them.
Lesser black-backed gulls have previously been shown to select flight routes to encounter beneficial orographic uplift, suggesting lesser blackbacked gulls have knowledge of their aerial surroundings (see Sage et al., 2019). Like thermal uplift, orographic uplift depends on landscape features, but in relation to wind strength and direction rather than earth surface energy balance. Thus, while gulls may learn to use the landscape to make use of these two sources of uplift and similar landscape features may provide thermal or orographic uplift (for example urban areas (Shepard et al., 2016)) the dynamics of the two uplift landscapes differ. Thermal uplift increases with solar radiation throughout the day, generally peaking at midday before declining, while changing weather systems, overcast days or precipitation result in less thermal uplift (see video Supplementary material S.6 for dynamic visualisation of the entire study period). Future studies which integrate both thermal and orographic soaring strategies and examine their spatiotemporal rhythms (Bohrer et al., 2012;Duerr et al., 2012;Santos et al., 2017) may be able to identify when and where good conditions for each strategy overlap or diverge, and determine how gulls make decisions between different sources of energy saving flight.

Conclusions
Lesser black-backed gulls can exploit the knock-on atmospheric effects of a human engineered landscape by favouring areas of high thermal soaring suitability and using thermal soaring flight when conditions are suitable to lower their cost of transport. The North-Holland landscape provides regular hotspots of urban thermal uplift orientated in such a way that birds can easily make connections across them, avoiding, for example, low wetland areas with less thermal soaring suitability. A mosaic landscape can still create movement corridors if the movement behaviour (in this case thermal soaring) is not dependent on continuous landscape features. The thermal soaring strategies we identify in this study are likely to occur in other landscapes or among other soaring species, but the prominence of urban hotspots in relation to the rest of the landscape is likely to influence the degree to which birds converge upon them. Furthermore, the connectivity of the landscape for birds depends not only on the structure of the landscape but also on time and energy trade-offs birds make in relation to their internal drivers and navigational capacity. Urban areas provide attractions to opportunistic species like lesser black-backed gulls, particularly for foraging (Pais de Faria et al., 2021; Spelt et al., 2019) and soaring raptors depending less upon urban food resources, or migratory soaring birds that possess less local knowledge, likely have different energetic trade-offs to make and may use different strategies to navigate the same thermal landscape as local breeding gulls. Therefore, by studying a wide array of energy landscapes and movement systems it will become more feasible to identify the point of trade-off for certain behaviours, such as the limits to which gulls will go in order to encounter thermal soaring hotspots in a mosaic landscape, the conditions under which gulls switch between different energy saving flight modes, or even whether thermal soaring opportunities may lead to range expansions by facilitating longer trips with lower movement costs (Bocedi et al., 2014). These trade-offs can also be examined under different ecological circumstances, e.g. during reproduction, migration, or exploratory movements to understand how internal and external drivers interact to influence movement decisions. Integration of all this information will then give insight into the true impact of human developed landscapes on the flight environment and its consequences for many flying species.
CRediT authorship contribution statement ES, WB and JSB conceptualised and designed the study. WB and WvD conceptualised the sensible heat flux model and developed the methodology with WvD carrying out the modelling, and validation of the sensible heat flux model. ES implemented the model on GPS data and carried out all other formal analysis. Field work was carried out by KCJC and long term tracking research by KCJC and JSB. Funding was acquired by JSB. Writing was led by ES, supported by WB and JSB, with all authors reviewing and editing the manuscript.

Data availability
Data will be made available on request.

Declaration of competing interest
The authors declare no competing interests.