Learning to predict spatiotemporal movement dynamics from weather radar networks

Weather radar networks provide wide‐ranging opportunities for ecologists to quantify and predict movements of airborne organisms over unprecedented geographical expanses. The typically sparse spatial distribution of radar measurements poses, however, a major challenge to spatiotemporal predictive modelling. We propose FluxRGNN, a recurrent graph neural network that is based on a generic mechanistic description of population‐level movements across the Voronoi tessellation of radar sites. The resulting hybrid model capitalises on local associations between environmental conditions and animal density as well as on spatiotemporal dependencies inherent to the movement process. We applied FluxRGNN to make 72‐h forecasts of nocturnal bird migration over Western Europe using simulated trajectories and measurements from the European weather radar network. For both datasets, FluxRGNN achieves higher predictive performance than baseline models based on environmental conditions alone. It effectively disentangles local take‐off and landing dynamics from aerial movements and correctly predicts migration directions with an accuracy of 87%. Continental‐scale forecasts of animal density and biomass fluxes have the potential to improve the impact and cost‐effectiveness of wildlife management and conservation efforts. With FluxRGNN this becomes feasible for nocturnal bird migration. In the future, other migration systems could benefit from applying the proposed method to similar static sensor networks.

removal of man-made barriers, obstacles, or sources of pollution (Liechti et al., 2013;Marschall et al., 2011;van Doren et al., 2021), to preventing local transmission of infectious diseases (Acosta et al., 2021;Lisovski et al., 2018), and protecting human infrastructure and agriculture from damage through mass movements (Leskinen et al., 2009;van Gasteren et al., 2019). Focusing such efforts on specific locations and times where interventions are most needed can greatly improve the ecological and societal impact while reducing the associated costs . This is particularly relevant for migratory species whose seasonal movements can span hundreds and thousands of kilometres, making it impracticable to implement management strategies without any form of spatial and temporal prioritisation (Runge et al., 2014). To effectively inform stakeholders about potential hotspots along a migration route, it is essential to have access to both adequate data and methods for spatiotemporal predictive modelling over vast geographical extents.
With the advent of ecological sensor networks and the continuing development of existing earth observation infrastructure, unprecedented opportunities are opening up to monitor and analyse the movements of terrestrial, aquatic, and aerial species at regional, national and continental scales (Farley et al., 2018;Kays et al., 2020).
A prominent example is the use of weather radars to study the aerial movements of birds, bats and insects Shamoun-Baranes et al., 2014). Although designed for meteorological purposes, these radars can provide high-resolution information on the local density and velocity of flying organisms around the radar station. The large spatial coverage of operational weather radar networks has proven particularly valuable for quantifying broad-front nocturnal bird migration at the continental scale (Dokter et al., 2018;Nilsson et al., 2019).
Despite recent technological and methodological advances that have enabled researchers to manage large amounts of data generated by weather radar networks and to translate raw observations into biologically relevant measures Lin et al., 2019;Stepanian et al., 2016), the step towards spatiotemporal predictive models of bird migration remains challenging. On the one hand, this is related to the complexity of the process. In contrast to well-understood physical processes like fluid flow, fully mechanistic models describing the movement of individual birds in response their environment (Aurbach et al., 2020;Liechti et al., 2013) struggle to capture the dynamics of bird migration in sufficient detail. Moreover, these individual-based models remain difficult to calibrate and involve computationally expensive simulations which make them unsuitable for near real-time predictions. To generate reliable and timely forecasts, more flexible and scalable data-driven models are required. On the other hand, however, range-dependent radar biases limit such approaches to vertical profiles (Buler & Diehl, 2009;Dokter et al., 2011), that is point estimates of bird density and velocity per radar and altitude. The resulting spatially sparse and semi-structured datasets make it difficult to model the underlying spatiotemporal movement process in a data-driven way.
Previous studies have explored the use of species distribution models to tackle the aforementioned issues (Costa et al., 2014;Elith & Leathwick, 2009;Soldevilla et al., 2011;van Doren & Horton, 2018). These models exploit the correlation between environmental variables and animal density at a particular site to reduce the complex problem of spatiotemporal prediction to a local regression task. Except for correlations in the predictor variables, they can thus not account for the spatial and temporal dependencies that are inherent to animal movements. This means that, firstly, they are limited in their predictive power because they cannot exploit relevant spatiotemporal patterns in the data, such as the onset and propagation of a migration wave, that may provide valuable information on how animal densities will change over time. And secondly, they cannot be used in contexts where directional information on migratory fluxes is crucial for management strategies, for example with pathogen spreading. Although this shortcoming has been recognised before, attempts to make species distribution models spatially explicit remain tailored towards datasets that can easily be projected on regular grids (De Marco et al., 2008;Domisch et al., 2019;Ovaskainen et al., 2016;Thorson et al., 2015) and are thus not directly applicable to weather radar data. Conversely, geostatistical methods that exploit spatial and temporal correlations in the data have been used for high-resolution interpolation of bird densities based on sparse weather radar measurements (Nussbaumer et al., 2019). While this approach could in principle be extended to generate predictions based on environmental conditions, it does not lend itself to big data applications  or capturing biomass fluxes and associated physical constraints.
The objective of this study is to incorporate both spatial and temporal dependencies into large-scale bird migration forecasts based on environmental conditions and spatially sparse observations obtained from weather radar networks, and to demonstrate the opportunities that state of the art machine learning techniques have to offer in this regard. To this end, we develop a hybrid modelling framework that combines a generic mechanistic description of migratory movements with deep learning techniques to account for (i) the temporal auto-correlation in the measurements of a single radar, (ii) the spatial dependencies between distributed radars that arise from migratory fluxes, and (iii) the influence of past and present environmental conditions on migration timing and direction.
Deep learning, a class of machine learning algorithms that use artificial neural networks to automatically extract complex features and relationships from large datasets, has revolutionised data-driven modelling in many different application domains (Goodfellow et al., 2016). While deep neural networks are extremely flexible, they require large amounts of data and are prone to make inconsistent predictions. Incorporating physical constraints and other prior knowledge about the data-generating process into the model architecture can help overcome these problems, leading to lower model complexity, i.e. lower data requirements, and more realistic and interpretable predictions (De Bézenac et al., 2019). Graph neural networks, for example, can use information about the existence or absence of arbitrary spatial dependencies to model the relations between many similar objects or locations (Battaglia et al., 2018;Zhang et al., 2020). This approach has been successfully combined with recurrent neural networks for temporal sequence modelling to solve real-life spatiotemporal problems such as traffic forecasting on road networks (Yu et al., 2018) or forecasting air and water quality based on environmental sensor networks (Liang et al., 2018).
Building on these advances, we designed a recurrent graph neural network, called FluxRGNN, that encodes spatial and temporal dependencies resulting from a generic mechanistic model of migration. In analogy to fluid dynamics, the mechanistic model describes migratory movements in terms of biomass fluxes between discrete geographical areas, here defined by the Voronoi tessellation of radar locations, leveraging continuity and mass-balance constraints that are inherent to migration processes. Previous studies have successfully applied this idea to estimate seasonal movement rates in bird migration networks from individual tracking data (Kölzsch et al., 2018), to infer residence time from daily counts of marked migrating birds (Drever & Hrachowitz, 2017), and to estimate biomass flows of nocturnal bird migration from weather radar data . Our approach differs from these studies in that we do not infer movement rates retrospectively but make (near real-time) predictions of migratory fluxes based on functional relations learned from data. The resulting hybrid modelling framework benefits from the flexibility and scalability of deep neural networks and at the same time leverages generally applicable system knowledge, providing opportunities for biological or environmental interpretation without relying on hand-crafted features or strong assumptions on species-specific behaviours.
As a proof of concept, we apply FluxRGNN to predict the nocturnal bird migration over Western Europe. We validate the proposed model based on simulated data as well as on observations from 22 weather radars in Germany, the Netherlands and Belgium, and compare its predictive performance with a range of baseline models including local seasonal trends as well as species distribution models.
To show the relevance of both temporal and spatial dependencies, we conduct ablation experiments where we remove different model components and evaluate the respective effect on the predictive performance.

| MATERIAL S AND ME THODS
Given a sequence of bird density and velocity estimates from a weather radar network together with past and future environmental conditions, we predict a sequence of future bird densities resulting from directed migratory movements across the observed spatial domain. We use a generic mechanistic description of population-level migratory movements to define the dependency structure of bird densities at the radar locations. The specific functional relations are then estimated based on data using a set of deep neural networks which are coupled according to the mechanistic model structure.
This results in a recurrent graph neural network architecture that is completely differentiable and can thus be fitted to weather radar network observations using standard optimisation methods for machine learning. Figure 1 shows an overview of the overall modelling framework.

| Migratory movements as fluid flow
The macroscopic patterns of broad-front bird migration arise from a large number of individuals pursuing similar strategies to move from one geographical area to another. In analogy to the continuum assumption in fluid dynamics, we abstract away from these individual trajectories and consider bird migration as a continuous process in space and time. This process can be described in terms of bird density [birds km −3 ] and a dynamically changing velocity field v [km h −1 ] along which migrants move . The time evolution of the total number of individuals in any given volume V i is then given by the continuity equation in integral form In words, any change in abundance in V i can be explained by (i) migratory movements into or out of the volume, or (ii) individuals entering and leaving the sky during departure and stop-over. The former is expressed by the surface integral ∫ V i ( v) ⋅ ndA i which defines the net migratory flux across the volume surface V i as the rate at which migrants leave minus the rate at which migrants enter the volume through infinitesimal surface elements dA i . These rates depend on the bird density as well as the velocity field v in relation to the corresponding surface normal vector n. The take-off and landing dynamics are encapsulated in the volume integral ∫ V i s dV i , where s is a local source/sink term that captures the number of migrants that are appearing or disappearing at a given location within the volume boundaries.

| Spatial discretisation
To model the dynamics of the entire system, we partition the spatial domain of interest into a number of nonoverlapping volumes and apply Equation 1 to each volume respectively, which results in a coupled system of ordinary differential equations. Based on the finite volume method from computational fluid dynamics (Versteeg & Malalasekera, 2007), we then numerically integrate the net migratory flux through the boundary surface of V i : where F j→i is a numerical flux term that approximates the true mi- The predicted terms are combined into animal density forecasts according to a mechanistic description of population-level movements. and the neighbouring volume V j . Importantly, F is conservative, i.e.
F j→i = − F i→j . This ensures that our model is consistent in the sense that any migrants entering volume V i through f ij must leave V j through the same face.
A major advantage of the finite volume method is its flexibility regarding the size and shape of the so-called control volumes or cells, on which the fluxes of a conserved quantity are modelled.
This makes it particularly suitable for the spatially sparse and semistructured observations from weather radar networks. Here, we define the set of control volumes based on the two-dimensional Voronoi tessellation (Ata et al., 2009;Beni et al., 2010) of radar locations (see Figure 2b), where each Voronoi cell represents the column of air above its polygon base. This ensures that each control volume corresponds to exactly one point measurement and no unobserved volumes are introduced as it would be the case with a regular grid.
The total number of migrants in each volume V i can then be approximated by using the vertically integrated bird density (VID) ̂ i [birds km −2 ] measured by the corresponding radar as an estimate of the average VID in the Voronoi cell with area A i : that approximates the volume integral ∫ V i s dV i for each volume respectively as the difference between the number of migrants taking off, N source(i) , and the number of migrants landing, N sink(i) .

| Temporal discretisation
Since weather radars take measurements at discrete time points, we further discretise the finite volume formulation by applying an explicit time integration scheme. That means, given measurements of the initial abundances N Discretising the continuous dynamics can result in numerical instability or untenable computational costs when the time resolution Δt is too low or too high with respect to the spatial tessellation and movement speed (Moin & Mahesh, 1998). The irregular geometry of the Voronoi tessellation makes it particularly difficult to find the optimal Δt that reduces the risk of migrants traversing multiple cells within one time step and thereby violating the continuity assumption, while ensuring sufficiently large fluxes between neighbouring cells and preventing high computational costs due to many small system updates. To achieve a good trade-off, we advise choosing Δt such that the speed v max and the minimum distance between any two sensors d min is close to 1.

F I G U R E 2
Weather radar locations in Germany, Belgium, the Netherlands, France and Switzerland with (a) 500 simulated trajectories of birds migrating from a line transect between 5848′N, 1436′E and 5012′N, 1954′E towards the destination area (surrounded by dashed line) in southern France, and (b) the corresponding 2D Voronoi tessellation using 25 unobserved boundary cells around the region of interest.

| Learning temporal dependencies
We assume both S (t→t+1) i and F (t→t+1) j→i to depend on previous environmental conditions and migratory movements in the respective volume V i . To learn these temporal dependencies, we use an encoder-decoder setup consisting of two recurrent neural networks, encoder and decoder , with a long short-term memory (LSTM) architecture (Hochreiter & Schmidhuber, 1997). LSTMs hold an internal state, or "memory", that is updated repeatedly as new information is fed to the network. In this way, temporal patterns, such as the accumulation of migrants due to poor migration conditions, can be extracted from sequential data. To map a sequence of past observations at time points t 0 − C, … , t 0 to a sequence of predictions for future time points t 0 + 1, … , t 0 + H, it is common to use two LSTMs, an encoder and a decoder, stringed together. The encoder iteratively processes the observations from past time points, here environmental conditions x i , bird velocities v i and abundances N i , and at every time step t → t + 1 incorporates the relevant information into its internal state Importantly, the same neural network is applied simultaneously to all cells. This greatly reduces the overall model complexity and the amount of data required to learn the network parameters. To account for spatial variability, the sensor position p i is included as an additional input to encoder .
After processing all past observations, the final encoder state z (t0) i is fed into the decoder, which iteratively combines it with additional information about the current time t + 1, here environmental forecasts x (t+1) i , as well as the previous (predicted) abundance N (t) i , to update its own internal state which is a high-dimensional latent representation of the information available about the past and present of V i , that is relevant for predicting animal densities at time t + 1.
For each time step of the decoder, this latent representation is translated into an estimate of the source/sink term S (t→t+1) i , by applying a simple feed-forward neural network, or multi-layer perceptron (MLP), S . The net source/sink term is then given by Here, the source term N (t→t+1) source(i) is a direct output of S , whereas the sink term N (t→t+1) sink(i) is computed as the product of the total amount of migrants in the cell, N (t) i , and the fraction of migrants leaving the cell, sink(i) , with the latter being estimated by S . This decomposition ensures that not more than the available number of migrants is leaving the cell through a local process.

| Learning spatial dependencies
In contrast to S (t→t+1) i , the flux terms F (t→t+1) j→i describe spatial movements between neighbouring cells and thus depend on the state of both V i and V j , as well as on their relative position to each other. To learn these spatial dependencies, we introduce another feedforward neural network F , which models the fraction of migrants in V j that move to the neighbouring cell V i by accounting for the com- The spatial dependency structure resulting from the mechanistic description in Section 2.1 can be represented as a graph , namely the Delaunay triangulation corresponding to the Voronoi tessellation of radar locations. In this representation, the fluxes F (t→t+1) j→i can be interpreted as messages being sent along the edges of , that is, across the Voronoi faces. All incoming messages are then aggregated at each node, that is, Voronoi cell, and combined with local information to form a prediction of animal densities at the next time point. The specific functions associated with messages and node updates are captured by the neural networks encoder , decoder and F . This procedure follows the message passing framework underlying many graph neural network architectures (Gilmer et al., 2017). The overall model, defined by Equations 2-4, can thus be formalised as a recurrent graph neural network that explicitly incorporates the equations of the underlying mechanistic model into its architecture.

| Boundary conditions
Due to the finite size of weather radar networks, the number of observed cells and thus the modelled domain is limited. To ensure that our model can account for all possible migratory movements, including movements across the domain boundary, we introduce a set of unobserved boundary cells B which enclose the observed control volumes (see Figure 2b). At every prediction step t → t + 1, we extrapolate the decoder states h (t) k and predicted abundances N (t) k to these boundary cells by averaging over all observed neighbouring cells. Assuming that all relevant environmental variables are available for the otherwise unobserved boundary cells, the boundary flux terms F (t→t+1) j→k and F (t→t+1) k→j with j ∈ B and k ∈ B can then be computed as before. Since bird densities are spatially correlated (Nussbaumer et al., 2019), the introduced errors are expected to be negligible compared to the gains in terms of physical consistency and additional information on environmental conditions.

| Data
In autumn more than 800 million nocturnally migrating birds fly across western Europe, heading towards their wintering grounds in southern Europe and Africa . Along this flyway a multitude of weather radars continuously monitor the sky and provide information on bird densities and velocities. We use this infrastructure to train and evaluate our model using both simulated and real-world data of nocturnal bird migration across Western Europe between 1 August and 15 November 2015, 2016 and 2017.

| Time resolution
For weather radars in Western Europe, the average distance to the closest radar is 118km, with a minimum of 61km and a maximum of 179km. The maximum ground speed of migratory birds measured by these weather radars is approximately 20ms −1 , or 72kmh −1 . We thus chose a time resolution of 1 h for all experiments to address the trade-off discussed in Section 2.1.3.

| Simulated data
For validation purposes, we simulated the movements of 100,000 migrating birds using an individual-based model. For simplicity, we limited movements to the two-dimensional plane. A detailed description of the behavioural rules and model parameters is provided in Supporting Information C.
The simulated trajectories (see Figure 2) were converted into artificial radar measurements using the locations of 48 operational weather radars distributed over Germany, the Netherlands, Belgium, France and Switzerland. For every cell in the corresponding Voronoi tessellation, we generated hourly measurements of the number of migrating birds and the average velocity vector, assuming that each radar perfectly observes its Voronoi cell. In addition, we computed source(i) and N (t→t+1) sink(i) based on the number of birds moving between Voronoi cells and the number of birds departing and landing in each cell respectively.

| Weather radar data
To evaluate our modelling framework in a real-world setting, we used data from 22 C-band Doppler weather radars in Germany, the Netherlands, and Belgium. The vol2bird algorithm (Dokter et al., 2011) with a radar cross section of 11cm 2 per bird was applied to the raw polar volumes to extract vertical profiles of bird densities and velocities based on data points within a range of 5-25 km from the radar location. Finally, the profiles were vertically integrated from 200m above the radar elevation up to 5000m above sea level and aggregated in time to obtain point measurements of the migratory activity above each radar at hourly time intervals. Assuming a uniform distribution of animal densities within each Voronoi cell, vertically integrated densities were extrapolated to the entire volume by multiplying with the corresponding cell area (see Section 2.1.2).
Since we were interested in nocturnal migrations, all measurements between civil dawn and dusk were set to zero.

| Environmental data
To account for the effect of dynamically changing atmospheric conditions on migration timing and direction, we use a range of weather variables retrieved from the ERA5 reanalysis dataset (Hersbach et al., 2020) as inputs to our model. For the simulated data, this includes the same u and v wind components at 850hPa (ca. 1500m a.s.l.) that were used as input to the individual-based model. For the real-world weather radar data, we additionally include 10m surface winds, temperature, cloud cover, specific humidity, total (cumulative) precipitation, surface pressure, and sensible heat flux. All variables were aggregated into a single value per Voronoi cell and time point, using the average of 100 randomly sampled locations.
Furthermore, we include a set of time-related variables to account for more general daily and seasonal trends. This includes the day of the year, solar position and its derivative, a binary variable indicating if it is day or night, and two binary variables marking the time points of civil dusk and dawn. Together, these variables form the vector of environmental conditions x (t) i for cell V i and time point t , which is used as input to both the temporal and spatial components of our model. beyond the training scope. Sequences with more than 10 % missing data were excluded, resulting in a training set of 4897 sequences for the simulated data and 3772 sequences for the radar data.

| Training procedure
During training, the mean squared error (MSE) over all cells and time points of the training set is computed to determine the mismatch between the predicted and observed bird densities. The model parameters are then adjusted using the Adam optimiser (Kingma & Ba, 2015) to reduce the mismatch. This procedure is repeated until the MSE for the validation set stops improving by more than 10 −6 every 50 iterations, or until a maximum of 300 iterations has been reached.
To facilitate a more stable training procedure, we used scheduled sampling (Bengio et al., 2015) and rescaled all input variables linearly to a range between 0 and 1. For all spatial features, such as sensor coordinates, wind, and animal velocity vectors, each dimension was rescaled by the same factor to retain directional information.

| Predictive performance
For both simulated and radar data, we evaluated FluxRGNN's predictive performance in terms of commonly used metrics for forecasting problems, namely root mean square error (RMSE) and Pearson r correlation coefficient, evaluated based on normalised hourly bird densities. Since FluxRGNN is trained using a stochastic learning algorithm, we repeated the training process five times using varying random seeds and report the mean and standard deviation over these trials for each evaluation metric, respectively.

| Comparison to baseline models
We compared the predictive performance of our modelling framework to three baseline models, ranging from a simple historical average (HA), to a generalised additive model (GAM) similar to Kranstauber et al. (2022), capturing daily and seasonal trends based on temporal features alone, to a more complex species distribution model based on gradient-boosted regression trees (GBT) similar to van Doren and Horton (2018). For GBT, we performed a cross-validated grid search to determine the best hyperparameter settings (see Supporting Information B).

| Ablation experiments
To investigate the relevance of the different model components, we performed three ablation experiments, a common approach in deep learning research where individual components are removed and the performance of the diminished model is compared to the full model (Deneu et al., 2021). In particular, we considered model variants (i) without spatial dependencies, i.e. without migratory fluxes, (ii) without encoding past conditions, and (iii) without introducing unobserved boundary cells.

| Simulation study
FluxRGNN predicts bird densities by combining the source and sink of a cell with all fluxes across its boundary. Since weather radar measurements do not provide full information on these biomass fluxes, we used the simulated data to evaluate if FluxRGNN can accurately disentangle take-off and landing dynamics from aerial movements without being explicitly trained to do so. We aggregated all terms over the first predicted night, using only sequences from the test dataset that start at the first hour of a night, and quantified the agreement between simulation and prediction in terms of the Pearson r correlation coefficient. For aerial fluxes between neighbouring cells, we determined the distance of all faces to the domain boundary by computing the shortest path from either of the two Voronoi cells to any boundary cell in the associated adjacency graph and analysed how the correlation coefficient varies with this distance. In addition, we computed the accuracy with which the correct sign, that is, movement direction, was predicted.

| Radar case study
We demonstrate the use of FluxRGNN for ecological purposes based on predictions of bird densities and fluxes using the best performing model instance trained on the real-world weather radar data. In particular, a map of averaged predicted nightly fluxes for the autumn migration in 2017 is presented. We compare the apparent patterns to the general distribution of flight directions measured by all weather radars during the same period. Finally, we showcase an example prediction for a single radar and its Voronoi cell, illustrating the various biological and environmental insights FluxRGNN can provide.

| Model implementation
FluxRGNN is implemented in Python 3.8 using the packages pytorch and pytorch-geometric. The code, including all model implementations and scripts for training, testing and data analysis, are publicly available at https://github.com/Fiona Lippe rt/FluxRGNN.

| Comparison to baseline models
Due to the recurrent nature of FluxRGNN, the predictive performance generally decreases as a function of the forecasting horizon (see Figure 3). For both simulated and radar data, the two evaluation metrics show a strong deterioration within the first 5 h and level off afterwards. With the simulated data, another drop in performance occurs after 24 h, indicating that multi-night predictions (i.e. on average >24 h into the future) are more difficult than single-night predictions. With the radar data, this effect is less prominent. More details on nightly patterns of evaluation metrics are available in Supporting Information D.
In contrast, the considered baseline models do not include any temporal dependencies. Their predictive performance is thus the same regardless of the forecasting horizon. Up to a forecasting horizon of 24 h, FluxRGNN outperforms all baseline models in terms of RMSE and Pearson r, for both datasets respectively. With the radar data, FluxRGNN's Pearson r eventually converges towards that of the species distribution model GBT, but does not drop below. In all other cases, FluxRGNN retains its leading position up to the maximum considered forecasting horizon. Given the relatively stable performance for forecasting horizons between 24 and 72 h, it can be expected that even for longer horizons FluxRGNN will produce useful predictions.  Interestingly, this correlation increases with the distance to the domain boundary, reaching values up to 0.87 ± 0.02 (see Figure 5a).

F I G U R E 3
Root mean square error (RMSE) and Pearson r as a function of the forecasting horizon, evaluated for baseline models (HA, GAM, GBT), FluxRGNN and three variants with different model components removed, for simulated (left) and radar data (right) respectively. Shaded regions correspond to the standard deviation over 5 model instances using different random seeds.
The per-cell quantities, that is, source, sink, influx and outflux, show a reasonable agreement between simulation and prediction.
The average Pearson correlation coefficient exceeds r = 0.5 for all four components (see Figure 5b), suggesting that our model learns both the spatial and temporal aspects of bird migration and does not focus on one over the other. Figure 5b shows that outbound components, that is, sink and outflux, reach lower correlation coefficients than inbound components. This can be attributed to a few cells close to the domain boundary (mostly over Switzerland and France) where predicted bird densities are generally less accurate, which in turn constrains the quality of predicted outfluxes and sink terms (see Sections 2. 2.1 and 2.2.2). Excluding all cells that border at least one boundary cell from the analysis increases the agreement between simulated and predicted outfluxes (from r = 0.69 ± 0.05 to r = 0.79 ± 0.04) and sink terms (from r = 0.54 ± 0.05 to r = 0.64 ± 0.02 ). We provide a more detailed breakdown of correlation coefficients by cell in Supporting Information D.  in Hannover, Germany. The shape of nightly migration pulses F I G U R E 4 Visual comparison of (a) simulated and (b) predicted nightly migratory fluxes, averaged over all sequences for autumn 2017. In only 0.5 % of the cases, simulated birds move to a nonadjacent cell within one time step, indicating that the chosen time resolution of 1 h is appropriate for the given movement process and radar network structure.

F I G U R E 5 Pearson correlation (a)
between simulated and predicted nightly fluxes between neighbouring cells, and (b) between simulated and predicted nightly source and sink, as well as nightly in-and outfluxes per cell. Error bars show standard deviations over 5 model instances.
for radar measurements and predictions is largely congruent. In the first predicted night, the radar detects migration occurring mainly within the first few hours after dusk. FluxRGNN predicts a very similar pattern. Due to the underlying mechanistic description, FluxRGNN can, however, provide additional insights into this pattern. After an initial phase of birds taking off from the ground (source > sink), our model suggests that the subsequent reduction in bird density is caused by birds leaving the cell and flying towards the south. In contrast, the relatively strong migration during the third forecasted night is predicted due to a long period of influx from the north-east (deboo and deumd) into night). As soon as the conditions improve these birds depart all together and fly towards south-west, leading to a strong influx into the region around dehnr.

| DISCUSS ION
FluxRGNN leverages information on environmental conditions together with spatial and temporal dependencies between spatially sparse measurements of weather radars to predict populationlevel movements over large geographical expanses. In contrast to other spatiotemporal deep learning approaches for sparse and irregular data (Liang et al., 2018;Seo et al., 2020), our framework is a hybrid between mechanistic and data-driven modelling. In this way we can effectively disentangle local take-off and landing dynamics from aerial movements across the Voronoi tessellation of radar locations. Our case study illustrates how this leads to comprehensive bird migration forecasts, which can be used to anticipate the spread and accumulation of migrants over time. Such insights are crucial to, for example, develop targeted containment strategies for the spread of pathogens like avian influenza (Acosta et al., 2021).

| Bird density predictions
As expected, FluxRGNN predicts bird densities more accurately than methods based on local correlations between bird density and environmental conditions alone. The results of our simulation-based ablation experiments indicate that FluxRGNN's spatial and temporal components are key to this predictive power. Nonetheless, when applied to real radar data, spatial fluxes do not lead to an improved predictive performance. This inevitably raises the question if a spatially explicit approach is beneficial for bird density predictions in practice, or if existing local regression models (e.g. van Doren & Horton, 2018) are sufficient as long as no information on biomass fluxes is needed. Yet, we find that radar measurements exhibit clear spatiotemporal correlations along the major migration direction (see Supporting Information D). This suggests that the propagation of migration waves is indeed captured in the data and could in principle be exploited for better prediction accuracy. Therefore, considering also the increasing agreement between predicted and simulated fluxes with distance to the domain boundary, we are confident that the spatial model component will become more effective as the radar network size increases.

| Migratory fluxes
The general migratory patterns predicted by FluxRGNN in our radar Poland, and leaving towards France. The comparatively strong predicted outflux over the Alps remains to be investigated further. We anticipate more realistic estimates of these fluxes when including additional environmental variables, such as landscape characteristics, and expanding the spatial domain southwards.

| Limitations
While modelling migratory movements on the Voronoi tessellation of weather radars allows us to explicitly infer biomass fluxes, it limits the spatial resolution of the model. This is in contrast to existing local regression models that can produce high-resolution bird density maps, such as BirdCast (van Doren & Horton, 2022). Moreover, the assumption that bird densities are uniformly distributed across each Voronoi cell makes the resulting predictions sensitive to potential radar biases. Applying a simple linear interpolation to the model inputs and outputs could alleviate this issue and translate the discontinuous Voronoi cell estimates into spatially smooth bird density maps.
Generally, current constraints to better exploit the potential of FluxRGNN are not computational but due to limited data availability. In particular, computational limits are not reached because scaling-up the spatial domain of FluxRGNN is straightforward by design: using GPU acceleration the model can make predictions at all sensor locations in parallel, and at the same time model complexity is contained by sharing parameter values among all sites. On the other hand, data is not easily accessible. Especially in Europe, obtaining unfiltered radar data from the many national meteorological services remains cumbersome. A major effort should be made to make this data openly accessible and interoperable (Shamoun-Baranes et al., 2021).

| Broader impact and future work
This study focused on modelling nocturnal bird migration based on vertically integrated observations from the European weather radar network. The proposed method can, however, easily be extended to work with full vertical profiles by defining multiple Voronoi tessellation layers between which birds can move vertically. Moreover, since the underlying mechanistic description is extremely general, FluxRGNN can in principle be applied to any kind of static sensor network, such as acoustic monitoring networks (Blumstein et al., 2011;Marques et al., 2009;Spillmann et al., 2015) and camera traps (Kays et al., 2009;Rowcliffe et al., 2008), providing data on densities or abundances of a moving population. Importantly, no a priori estimates of biomass flows are required since FluxRGNN infers these quantities in an unsupervised way. Moreover, while we use velocity measurements as additional inputs to the encoder neural network, this is not strictly necessary. The same model architecture can be used to estimate spatial fluxes from density measurements and environmental conditions alone. In general, there are two major requirements for FluxRGNN to be transferable to other domains. First, the space-time resolution of the data should match the scale of the underlying process. In particular, the continuity assumption has to be met in terms of cell sizes, time resolution and movement speed (see Section 2.1.3). Second, the spatial and temporal coverage of the data needs to be sufficient. For seasonal movements, at least 3 years of data should be available to avoid overfitting to a specific season.
Spatially, we recommend using sensor networks with no less than 20 sites to limit the effects of unobserved boundary cells. For systems that meet these requirements, FluxRGNN has the potential to generate comprehensive movement forecasts that inform stakeholders about relevant hot spots and biomass fluxes and thereby help improve the ecological and societal impact of wildlife management efforts.
Finally, we see various opportunities for leveraging newly emerging methods from machine learning to further advance the application of model-data fusion techniques to animal movement systems. Building directly on FluxRGNN, we see especially opportunities to increase spatial resolution and identify processes by post hoc analyses. For example, recent work by Iakovlev et al. (2020) shows promising results on deep learning unknown partial differential equations from sparse data points. Refining FluxRGNN to such a continuous setting will allow for more detailed predictions capturing also fine scale patterns related to topography, habitats and weather conditions, similar to van Doren and Horton (2018).
Furthermore, post hoc analysis methods from the field of explainable artificial intelligence (Adadi & Berrada, 2018) could be applied to better understand the relations between movement patterns and environmental conditions learned by the neural network components of FluxRGNN. In combination with the underlying mechanistic model, this could help fill knowledge gaps existing for many migration systems. Science Foundation (NSF 1927743). Computational support was provided by SurfSara.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/2041-210X.14007.

DATA AVA I L A B I L I T Y S TAT E M E N T
Preprocessed weather radar data, simulated data and environmental data, together with trained models and results are available at https://doi.org/10.5281/zenodo.6874789 .