Quantifying year-round nocturnal bird migration with a fluid dynamics model

To understand the influence of biomass flows on ecosystems, we need to characterize and quantify migrations at various spatial and temporal scales. Representing the movements of migrating birds as a fluid, we applied a flow model to bird density and velocity maps retrieved from the European weather radar network, covering almost a year. We quantified how many birds take-off, fly, and land across Western Europe to (1) track bird migration waves between nights, (2) cumulate the number of birds on the ground and (3) quantify the seasonal flow into and out of the study area through several regional transects. Our results identified several migration waves that crossed the study area in 4 days only and included up to 188 million (M) birds that took-off in a single night. In spring, we estimated that 494 M birds entered the study area, 251 M left it, and 243 M birds remained within the study area. In autumn, 314 M birds entered the study area while 858 M left it. In addition to identifying fundamental quantities, our study highlights the potential of combining interdisciplinary data and methods to elucidate the dynamics of avian migration from nightly to yearly time scales and from regional to continental spatial scales.

RN, 0000-0002-8185-1020; SB, 0000-0002-0844-164X; LB, 0000-0002-8182-0152; GM, 0000-0002-8820-2808; FL, 0000-0001-9473-0837; BS, 0000-0002-7736-7527 To understand the influence of biomass flows on ecosystems, we need to characterize and quantify migrations at various spatial and temporal scales. Representing the movements of migrating birds as a fluid, we applied a flow model to bird density and velocity maps retrieved from the European weather radar network, covering almost a year. We quantified how many birds take-off, fly, and land across Western Europe to (1) track bird migration waves between nights, (2) cumulate the number of birds on the ground and (3) quantify the seasonal flow into and out of the study area through several regional transects. Our results identified several migration waves that crossed the study area in 4 days only and included up to 188 million (M) birds that took-off in a single night. In spring, we estimated that 494 M birds entered the study area, 251 M left it, and 243 M birds remained within the study area. In autumn, 314 M birds entered the study area while 858 M left it. In addition to identifying fundamental quantities, our study highlights the potential of combining interdisciplinary data and methods to elucidate the dynamics of avian migration from nightly to yearly time scales and from regional to continental spatial scales.

Background
The sheer numbers of migratory birds create huge biomass flows [1][2][3] that impact ecosystem functions and human economy, agriculture and health through the transport of energy, nutrients, seeds, and parasites [4]. To understand these influences on ecosystems and make use of, or avoid, the resulting services and disservices, we need year-round and continental-wide monitoring of migratory fluxes and their quantification at various spatial and temporal scales. Continental networks of weather radars are increasingly becoming essential tools to monitor large-scale migratory movements [5]. However, most studies so far have focused on specific stages of the migration journey: migratory flights (e.g. [2,[6][7][8][9]) or stopovers (e.g. [10][11][12]). Yet, none have explicitly considered and differentiated between the three successive stages of take-off, flight and landing: we therefore lack a comprehensive model of the entire migratory journey. Broad-scale bird migration has already been modelled with agent-based models [13,14] and collective graphical models [15]; however, both these approaches rely on strong assumptions on the mechanisms of migration, can be hard to calibrate and often result in large uncertainties.
To integrate migratory take-off, flight and landing into a single framework, we adopted a methodology from fluid mechanics. While novel in aeroecology, fluid mechanics methods have long been used in movement ecology (e.g. [16][17][18][19]), with a few examples in bird movement ecology [20,21]. For instance, the spread of invasive bird species has been modelled with reaction-diffusion equations [22][23][24][25]. In general, these approaches aim to fit parameters describing fluid mechanisms (e.g. diffusion coefficient) to a sparse dataset of observations and learn the underlying physics of the phenomena. Here, we also modelled the nocturnal broad-fronted bird migration as a conserved quantity based on the continuity equation (e.g. [26]). However, we employed a complete dataset (in space and time) of bird density and velocity interpolated from weather radar measurements (figure 1) [8], allowing us to apply directly the advection equation without assuming an underlying process of movement (e.g. movement proportional to the gradient of a quantity). In addition, we added a sink/ source term (equivalent to the reaction term) to quantify how many birds take-off and land. Indeed, since we assume that the biomass of birds moving from one grid cell to another is conserved, any change of bird density in the air must be explained by movements to and from the ground.
With this framework, we can quantify how many birds take-off, fly, and land at any given time and location. In subsequent steps, we used the resulting maps of take-off and landing to (1) track waves of bird migration between nights across Europe, (2) estimate the accumulation (i.e. changes in numbers) of birds on the ground throughout the year and (3) quantify the seasonal flow in and out of the study domain through several transects.

Data
We used the data from 37 weather radars in France, Germany, The Netherlands and Belgium operating between 13 February 2018 and 1 January 2019. This dataset is currently the longest available time series over a large part of Western Europe. It consists of vertical profiles of bird density [birds km −3 ], flight speed  . First, we interpolate vertical profile time series of bird density and velocity field measured by weather radar data into continuous spatio-temporal maps following [8]. 2. Flow model ( §2.3). Then, using the interpolated data in a flow model allows us to estimate the number of birds entering, leaving, taking off from and landing in each grid cell at each time step. 3. Migration processes ( §2.4). The resulting maps of take-off and landing birds allow us to investigate the spatio-temporal variation of stopover, the accumulation of birds on the ground, and the geographical variation in the seasonal fluxes of migrating birds. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210194 [m s −1 ] and flight direction [°] which were generated with the vol2bird software [27] and are available on the ENRAM repository [28] at a 5 × 200 m (0-5000 m a.s.l.) resolution. Similar to previous studies [7,8], the vertical profiles were cleaned as follows (electronic supplementary material, 1.2). First, we eliminated high-reflectivity contamination (e.g. rain and ground scatter) using a dedicated graphical user interface. Then, we removed contamination from slow-moving targets with low reflectivity such as insects or snow based on standard deviation of radial velocity and air speed [29]. Finally, we vertically integrated bird density and flight speed (i.e. volumetric to areal) while (1) accounting for the impact of local topography on the surveyed volume and (2) simulating bird density in the volume of air below the altitude surveyed (electronic supplementary material, 1.3).

Interpolation and simulation
Since the radars provide point observations (averaged over a 5-25 km radius around the radar location), we interpolated bird density [birds km −2 ] into a spatio-temporal grid using the methodology developed in [8]. The bird velocity field (i.e. the vector field of birds' flight speed and direction) was interpolated for the two north-south and east-west components separately using a similar methodology. Adjustments of the interpolation method to a year-round dataset and to a velocity field are detailed in electronic supplementary material, 2.
The interpolation grid was defined between latitudes 43°and 55°and longitudes −5°and 16°, with a resolution of 0.25°a nd between 13 February 2018 and 1 January 2019 with a resolution of 15 min in time. Similarly to [8], grid cells were excluded if (1) they were located over a water body or above 2000 m a.s.l., (2) they were more than 150 km away from a weather radar, (3) they spanned over day time (i.e. from sunrise to sunset), or (4) rain intensity exceeded 1 mm h −1 (interpolated from ERA5 dataset from [30]). Nights without data were excluded from the interpolation (5 nights in early April, and 34 nights in July-August). The resulting interpolation maps can be visually explored at https://birdmigrationmap.vogelwarte.ch/2018/.
To correctly calculate aggregated measures (e.g. average bird density or sum of birds take-off) and their uncertainties, we generated 500 geostatistical simulations of bird density [8,31]. Indeed, since the estimation (i.e. Kriging) provides 'the most likely value' at any point in space and time, estimation maps are smoother than the real process, and aggregating such smooth maps (i.e. sum or mean in time and/or space) would produce unreliable values [31,32]. On the other hand, simulations represent 'one-of-many possible value', so that together they adequately capture the range of possible spatio-temporal patterns of the process (e.g. how bird density changes in time and space) and we can use them to compute the mean and uncertainty of aggregated operations.

Flow model
Based on the principle of mass conservation, the continuity equation (equation (2.1)) describes the transport of a conserved quantity (e.g. bird density): the rate of change of this quantity is equal to its flux into and out of a given volume (e.g. sky). The equation can also include a source/sink term, which accounts for the appearance (and disappearance) of the quantity (e.g. take-off and landing). The differential form of the continuity equation for bird density ρ [birds km −2 ] is where v = [v lon , v lat ] is the bird velocity field [km h −1 ] along latitude and longitude and W is the source/sink term [birds h −1 km −2 ] and r denotes the vector differential operator. The continuity equation is discretized with a forward time centred space (FTCS) scheme [33]. The source/sink term can be computed for each cell (i, j, t) with and discretized in longitude, latitude and time with the indexes i, j, t, respectively. Δlon, Δlat and Δt are the grid resolution in longitude, latitude and time, respectively. We applied this model to bird migration by using the spatiotemporal maps of bird density (ρ) and flight speed vector (v) derived from geostatistical simulations ( §2.2). The local fluxes were computed for each grid-cell by multiplying the density with the flight vector and then linearly interpolated to the grid cells' boundaries for both the longitudinal and latitudinal components. As the grid was defined in equal latitude and longitude intervals, we account for the varying resolution of Δlon in km along the latitude axis in the discretization. Finally, using equation (2.2), the source/sink term was computed for each grid cell at each time step as the change of bird density over time minus the spatial difference of fluxes. The source/ sink term W was composed of birds taking-off and landing (within the study area) which can be separated according to the sign of W. Indeed, as the reference of the mass balance was the sky, positive values of W correspond to birds taking-off while negative values correspond to landing. Additionally, the values of the fluxes at the study area's boundaries were extracted as the number of birds entering ( positive) and leaving (negative) the study area. We used the 500 simulations to produce space-time maps of (1) take-off and landing [birds km −2 ] and (2) fluxes in latitude-longitude [birds h −1 km −1 ] at the boundaries of the study area.

Migratory processes
The resulting maps were processed to address specific ecological questions. We were particularly interested in characterizing and quantifying nightly migration pulses and stopovers, the accumulation of birds on the ground, and the seasonal migration flows. To achieve this, we processed each of the 500 simulations as follows: -The nightly migratory pulses and stopovers were calculated by summing the take-off and landing movements separately over each night, and by visually comparing the maps of landing in the morning with those of take-off the following evening.
-The year-round accumulation of migratory birds on the ground was quantified by first aggregating the four fluxes (take-off, landing, entering, leaving) over the whole study area and for each night. Then, the nightly change in the number of birds on the ground was computed as the difference between landing and take-off, or equivalently, between entering and leaving. The cumulative sum of these daily changes corresponds to the number of birds that remained on the ground. We arbitrarily set the starting value of the accumulation to zero because the initial number of (resident and/or wintering) birds on the ground is unknown. -The seasonal flow of bird migration is quantified by summing the fluxes of birds entering and leaving the study area over spring (February-June) and autumn (August-December). To capture the variability of movements across Europe, we defined six transects along the boundary of the study area according to the major flyways: UK, the North, the East, the Alps, Spain and the Atlantic ( figure 4).

Nightly migratory pulses and stopovers
For illustration purposes, we selected a well-defined migration wave spanning from 6 to 10 April, during which birds moved from south-western France to north-eastern Germany (figure 2). The nightly averaged bird density and flight speed was highest between the main take-off and landing areas. More importantly, one night's landing and the following night's take-off were in good agreement, demonstrating that the data and proposed methodology can accurately track a wave of migration over several days. This agreement was particularly striking in this example because most birds continued their migration during subsequent nights. The crossing of the study area in approximately 4 nights corresponds to nightly migratory bouts of around 300 km. The same nightly take-off and landing maps are available for the entire year in the electronic supplementary material.

Seasonal flow
The bird migration in the study area (both in spring and in autumn) was mainly directed between Spain and eastern Germany (figure 4). Indeed, even the migration through the Atlantic transect mostly comprises birds crossing the Bay of Biscay from/to Spain. For instance, 91% of all movements across the Spanish transect were in-movements in spring, i.e. most birds entered the study area rather than left it, and similarly, 90% of all movements in autumn were out-movements of birds leaving the study area. By contrast, movements across, for example, the  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210194 Alps transect were less unidirectional in both seasons with only 63% of all movements in spring being movements into the study area (entering) and similarly, 72% of all movements in autumn were movements out of the study area (leaving) (electronic supplementary material, figure 3.2). Secondly, the timing of migration differed between transects (electronic supplementary material, figure 3.2), with, for example, Spain and Atlantic transects seeing more than half of their migration before mid-March, while only 20-30% of birds had crossed the East and North transects at that time.

Discussion
In this study, we presented a novel methodology inspired from fluid dynamics to model the flow of nocturnal migrants, from take-off, during nocturnal flight, to landing. The model produces high-resolution maps that enable us to investigate the dynamics of migratory movements at various spatial and temporal scales. We used the largest dataset available on the ENRAM data repository to characterize and quantify nightly, seasonal and year-round migration patterns over most of Western Europe.

Model
The model presented in this study builds on the methodology developed in [8], which interpolates point observations of bird densities measured by weather radars into continuous maps. We used these maps as the input for a flow model by considering bird migration as a fluid. This allows us to extract more dynamic information about bird movements and, in particular, their take-off and landing.
The approach used in this study models bird flow (i.e. average bird movement) rather than individual birds. Indeed, the weather radar data consist of bird density and flight speed averaged over a 25 km radius. Therefore, the estimated flows cannot capture the properties or behaviour of individual birds. For instance, the speed of individual birds is typically higher than the speed of the overall flow, and the flight directions of individual birds are more variable than the overall flow direction. Similarly, the modelled flow is unable to track separately multiple bird populations simultaneously migrating in different directions.
Throughout each step of the methodology, we identified the sources of errors and tracked the corresponding uncertainties to reliably estimate the ranges of the model output, as detailed below. Despite our best efforts to clean the data (electronic supplementary material, 1.2), there is an inherent uncertainty in the weather radar data (e.g. ground scattering, measurement errors, radar biases). We partially accounted for the data uncertainty at low altitude by generating uncertainty ranges in the vertical integration (see electronic supplementary material, 1.3). For more details on the data quality of weather radars, we refer the readers to the assessment and comparison found in [8,36,37]. We handle these unknown errors in the geostatistical framework (i.e. interpolation) by fitting a nugget effect in the spatio-temporal model (more details in electronic supplementary material, 2). The nugget effect essentially accounts for the random noise in the data (e.g. due to the data error or uncertainties), and then allows the interpolated values to diverge from datapoints. In addition to the data error (i.e. difference between 'true' passage and measured passage), the nugget effect also models small-scale variability and/or discontinuity in bird densities caused by geographical features (mountains, deserts, seas) or weather conditions (e.g. rain). We propagated these uncertainties throughout our methodology by generating 500 simulations of bird density representing the range of possible values according to the uncertainty ( §2.2). Each simulation is run in the flow model, and provides a distribution of the possible values for each output (e.g. number of birds on the ground).

Stopover
In this study, we demonstrated how waves of bird migration at the regional scale can be tracked over multiple nights ( figure  2). Our flow model links birds on the ground with birds in the air and can thus quantify the fluxes of take-off, flight and landing. Looking ahead, this example suggests that a forecast system based on a flow model of bird migration could accurately predict birds' landing during the night, and perhaps, on a longer term, take-off and landing over a few days.
Our method can compute the rates of both take-off and landing at higher spatial and temporal resolution than earlier approaches (e.g. 3 h after sunset in [38], interpolation at civil twilight in [39,40] or at maximum density within 2 h after sunset in [13])-a feature that will be particularly useful in follow up studies that link movements and stopovers to geographical features or short, intense weather events.
Although our model can identify the places and times where birds stopover, other aspects of stopover dynamics such as stopover duration or survival remain to be tackled in future multidisciplinary studies. The main obstacle to addressing stopover dynamics is the inability to track birds during the day, i.e. which of the birds landing one day are the ones taking-off the following day(s). A similar problem appears at the seasonal scale, where we cannot separate birds that are wintering, breeding or passing from the birds landing or taking-off. A potential solution would be to explicitly model stopover duration with a residence time model [20].

Accumulation and seasonal flow
Using our novel methodology and an almost continuous 1 year dataset, we assessed the relative changes in the number of birds on the ground. We estimated that in autumn 2018, 858 million birds (Q5-Q95: 797-917) migrated southward through Spain and over the Alps (including the Atlantic transect) ( figure 4). This number includes only nocturnal migrants, but both long-and short-distant migrants. The only previous quantification of migrant bird population estimated that between 1.52 and 2.91 billion long-distance migrants leave the entire European continent in autumn [3]. Our estimation is in agreement with these numbers if we consider that the area of origin of birds migrating through our southern transects corresponds to roughly one-third of the European continent (i.e. from the British islands to Scandinavia, Finland to Poland, and our study area). In North America, the number of birds migrating south of the USA in autumn was estimated to around 4.72 billion birds [2], which corresponds to an average density of 236 birds km −2 (for an area of 19:8 million km 2 for USA and Canada). Despite the differences of scale and ecological context, we found a comparable average density of 286 birds km −2 , again assuming that a third of the European bird population migrates through the southern transect. Note that these densities include only the nocturnal migrant birds royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210194 crossing the predefined transects (south of USA and southern transects in Europe).
The ratio between autumn and spring fluxes can be used to estimate an index of net recruitment, accounting for both reproduction and mortality [2]. For the USA, [2] estimated a ratio of 1.36 for a transect along the southern border and 1.60 for a transect along the northern border. In our study area in Europe, the resulting indices are 1.74 (Q5-Q95: 1.55-1.94) in the southern transects (Alps, Spain and Atlantic; figure 4) and 1.26 (1.09-1.42) in the northern transects (UK, North and East; figure 4). However, such derived values such as recruitment critically depend on birds taking similar migration routes in both spring and autumn. If migration routes vary between seasons, e.g. when birds take a more easterly route in autumn, recruitment numbers become distorted. Instead of computing the ratio of migratory birds flying across non-representative transects, we can take advantage of the flow model to estimate a ratio of migratory birds entering and leaving an area of interest, and thereby relate the recruitment index computed over this area to environmental characteristics. For the entire study area, a recruitment index of 2.26 (1.80-2.81) resulted from the ratio between the relative number of birds that have left in autumn (i.e. leaving minus entering) (544 million; figure 3) and the relative number that have arrived in spring (i.e. entering minus leaving) (243 million; figure 3). However, as the fluxes of wintering and breeding bird populations cannot be distinguished (see discussion on stopover), this recruitment index also depends on the number of wintering birds that leave the study area in spring and return in autumn with offspring. Therefore, while this recruitment index can characterize the migratory bird population growth, it cannot separate the influence of breeding and wintering populations. A possible avenue to address this challenge would be to combine breeding and/ or wintering bird atlas data with our accumulation of birds on the ground. This approach could provide absolute numbers of breeding, passing and wintering birds along with their corresponding recruitment indices.