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

2.1. 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⁻³], flight speed [m s⁻¹] 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).

2.2. 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⁻²] 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° and 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⁻¹ (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.

2.3. 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⁻²] is

where v = [v_lon, v_lat] is the bird velocity field [km h⁻¹] along latitude and longitude and W is the source/sink term [birds h⁻¹ km⁻²] and $\mathrm{\nabla }$ 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) withwhere $\mathit{\Phi }=\rho \mathbf{v}=[{\mathit{\Phi }}_{\mathrm{lon}},{\mathit{\Phi }}_{\mathrm{lat}}]$ is the flux term expressed in [birds km⁻¹ h⁻¹] 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 spatio-temporal 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⁻²] and (2) fluxes in latitude–longitude [birds h⁻¹ km⁻¹] at the boundaries of the study area.

2.4. 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).

3.1. 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.

• Download figure
• Open in new tab
• Download PowerPoint

3.2. Year-round accumulation of migratory birds on the ground

Summing the flow at the daily (or nightly) scale allowed us to characterize the year-round changes in numbers of migratory birds on the ground (figure 3). A seasonal sum of the bird movements (take-off, landing, entering, and leaving) is provided in supplementary material, figure 3.1. The number of birds on the ground rose steeply in March with almost 200 million more birds entering in the study area than leaving it, i.e. to migrate further North or East. Numbers then declined from August onwards and plummeted in mid-October. The number of birds on the ground became negative in autumn because our method did not explicitly account for reproduction and mortality, thus the new generation is included in the birds leaving the area in autumn. The spring migration period was shorter and more condensed (March–May) than the autumn migration (August to mid-November) (figure 3), with 50% of all take-offs taking place during 19 nights in spring and 29 nights in autumn. At peak migration, we estimated 118 (Q5–Q95: 99–137) million birds taking off in a single night in spring (30 March–1 April) and 148 (133–164) million in autumn (17–18 October). Prior to these two peak migration events, we observed that the accumulation curve of birds on the ground flattened, indicating temporarily reduced migratory traffic, possibly due to unfavourable weather conditions (Zugstau) [34,35].

• Download figure
• Open in new tab
• Download PowerPoint

3.3. 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.

• Download figure
• Open in new tab
• Download PowerPoint

In spring, 494 (Q5–Q95: 453–540) million birds entered the study area through the southern transects (Alps, Spain and Atlantic) and, at the same time, 251 (Q5–Q95: 228–273) million left it in the northern transects (UK, North and East), creating a surplus of 243 (Q5–Q95: 200–290) million birds that remained within the study area (figure 4). Similarly, in autumn, 544 (Q5–Q95: 480–610) million more birds departed than arrived: 314 (Q5–Q95: 284–344) million birds entered through the northern transects while 858 (Q5–Q95: 797–917) left through the southern transects. The ratio of the autumn deficit to the spring surplus is 2.2 (Q5–Q95: 1.8–2.8), meaning that for one bird staying in spring, two birds left the ground in autumn.

Compared to spring, birds took a more easterly route in autumn, with proportionally more birds flying through the Alps transect (autumn/spring = 114/53 = 2.15) than through the Atlantic transect (233/168 = 1.4). Moreover, nearly the same number of birds crossed the East transect in autumn and spring (150/138 = 1.1). Overall, this pattern could be indicative of a clockwise loop migration where birds migrate to their breeding areas via the Iberian Peninsula in spring and fly to their non-breeding areas further east in autumn.

The seasonal fluxes per transect summarize the number of birds entering and leaving and can therefore cover some fine-scale features of migration. Firstly, some transects showed a more unidirectional flow of migrants whereas the entering and leaving fluxes were more balanced for other transects. 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 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.

4.1. 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).

4.2. 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].

4.3. 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⁻² (for an area of $19.8\hspace{0.17em}\mathrm{million}\hspace{0.17em}{\mathrm{km}}^2$ for USA and Canada). Despite the differences of scale and ecological context, we found a comparable average density of 286 birds km⁻², 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 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.