A geostatistical approach to estimate high resolution nocturnal bird migration densities from a weather radar network

Quantifying nocturnal bird migration at high resolution is essential for (1) understanding the phenology of migration and its drivers, (2) identifying critical spatio-temporal protection zones for migratory birds, and (3) assessing the risk of collision with man-made structures. We propose a tailored geostatistical model to interpolate migration intensity monitored by a network of weather radars. The model is applied to data collected in autumn 2016 from 69 European weather radars. To cross-validate the model, we compared our results with independent measurements of two bird radars. Our model estimated bird densities at high resolution (0.2°latitude-longitude, 15min) and assessed the associated uncertainty. Within the area covered by the radar network, we estimated that around 120 million birds were simultaneously in flight [10-90 quantiles: 107-134]. Local estimations can be easily visualized and retrieved from a dedicated interactive website: birdmigrationmap.vogelwarte.ch. This proof-of-concept study demonstrates that a network of weather radar is able to quantify bird migration at high resolution and accuracy. The model presented has the ability to monitor population of migratory birds at scales ranging from regional to continental in space and daily to yearly in time. Near-real-time estimation should soon be possible with an update of the infrastructure and processing software.

collision risks with human-made structures and aviation to inform stakeholders.
However, the great majority of migratory landbirds fly at night (Winkler, 1999), 24 challenging the quantification of the sheer scale of bird migration. 25 Radar monitoring has the potential to quantify migratory movements of birds at the 26 continental scale (Drake & Bruderer, 2017). Initially limited to single dedicated 27 short-range measurements, the use of existing weather radar networks provide 28 continuous monitoring over large geographical areas such as Europe or North America 29 (Gauthreaux, Belser, & van Blaricom, 2003;Shamoun-Baranes et al., 2014), and led to 30 an upswing of radar aeroecology (Bauer et al., 2017;Dokter et al., 2018;Van Doren & 31 Horton, 2018). One important challenge in using networks of weather radars is the 32 interpolation of their signals in space and time. Recent studies (Dokter et al., 2018;Van 33 Doren & Horton, 2018) have used relatively simple interpolation methods as they 34 targeted patterns at coarse spatial and/or temporal scales. However, these methods are 35 insufficient if higher spatial or temporal resolution is wanted such as for the 36 fundamental and applied challenges outlined above. 37 To achieve high resolution interpolation of migration intensity derived from weather 38 radars ( 20km-15min), we propose a tailored geostatistical framework able to model the 39 spatio-temporal pattern of bird migration. Starting from time series of bird densities Locations of weather radars of the ENRAM network, whose fall 2016-data were used in this study (yellow dots), and their key characteristics (right panel). We used data from two dedicated bird radars -in Switzerland and France -for validation (red dots).
Based on the reflectivity measurements of these weather radars, we used the bird 59 densities as calculated and stored on the repository of the European Network for the 60 Radar surveillance of Animal Movement (ENRAM) 61 (github.com/enram/data-repository)(see  for details on the 62 conversion procedure). We inspected the vertical profiles and manually cleaned the bird 63 densities data (see detailed procedure in Supporting Information S1).

4/20
Since we targeted a 2D model, we vertically integrated the cleaned bird densities from the radar elevation and up to 5000 m above sea level. Because we aimed at 66 quantifying nocturnal migration, we restricted our data to night time, between local 67 dusk and dawn (civil twilight, sun 6 • below horizon). Furthermore, as rain might 68 contaminate and distort the bird densities calculated from radar data, a mask for rain   Figure 3. Illustration of the proposed mathematical model decomposition into trend, amplitude, curve and residual. Note that the power transformation was not applied on this illustration.
where s lat and s lon are latitude and longitude of location s, w lat is the slope coefficient 106 in latitude and w 0 is the value of the trend at the origin. Because no longitudinal trend 107 is observed in the data (Figure 1a), only latitude is used to parametrize the trend 108 function (see Figure S3-2 in Supporting Information S3). It is worth noting that if 109 longer periods are considered, Eqn. 2 should be replaced by a more complex parametric 110 function in order to handle the emerging patterns of long term non-stationarity.

111
Curve 112 The second non-stationarity visible in the dataset is the nightly pattern (Figure 2e-g) 113 that results from the onset and sharp increase of migration activity after sunset, and its 114 slow decrease towards sunrise (e.g. Bruderer & Liechti, 1995). Similar to the trend, this 115 7/20 pattern needs to be extracted from the original signal to avoid non-stationarity. This is 116 done using a curve template c for all nights and locations, defined as the polynomial of 117 where a i are the coefficients of the polynomial and N N T (Normalized Night Time) is 119 the standardized proxy of the progression of night defined as where t ↓ (s, t) and t ↑ (s, t) are the times of civil dusk and dawn respectively. N N T is 121 defined such that the local sunrise or sunset occur at N N T = −1 and N N T = 1, 122 respectively.

123
Amplitude 124 After removing the non-stationarities with the trend and the curve, the amplitude A 125 models the nocturnal bird densities at the daily scale as a stationary space-time random 126 process. Its value is therefore constant within a night at a given location but varies 127 between locations and between nights. It accounts for the correlation at the scale of

Bird migration mapping
After parametrisation, the geostatistical model can be used to interpolate bird density 138 observations derived from weather radars to produce high resolution maps. The full 139 mathematical description of the procedure is detailed in Supporting Information S2 and 140 a brief outline is given below.

141
The estimation of the bird density at any unsampled location Z(s 0 , t 0 ) * is performed 142 are modelled as random processes, the estimations of A(s 0 , t 0 ) * and R(s 0 , t 0 ) * are 147 performed by kriging (e.g. Chilès & Delfiner, 1999;Goovaerts, 1997). An important 148 advantage of using kriging is that it expresses the estimation as a Gaussian distribution, 149 thus providing not only the "most likely value" (i.e. mean or expected value) but also a 150 measure of uncertainty with the variance of estimation. Tracking back the uncertainty 151 provided by kriging to the final estimation Z(s 0 , t 0 ) * is non-trivial but possible through 152 the use of a quantile function (see Supporting Information S2 for details). Consequently, 153 the estimation Z(s 0 , t 0 ) * is expressed as the median and its uncertainty range is defined 154 as the quantiles 10 and 90. A continuous space-time estimation (with uncertainty) map 155 is then computed by repeating the procedure for estimating a spatio-temporal point 156 s 0 , t 0 on a discrete set of points (i.e. grid).
In the case study presented in this paper, both estimation and simulations maps are 168 calculated on a spatio-temporal grid with a resolution of 0.2 • in latitude (43 • to 68 • ) and 169 longitude (-5 • to 30 • ) and 15 minutes in time, resulting in 127x176x2017 nodes. Over this large data cube, the estimation and simulation are only computed at the nodes  traffic rates (MTR) and average speed of birds aloft (Nilsson et al., 2018;Schmid et al., 193 2019). 196 Cross validation

197
In our study, the normalized error of estimation over all radars has a near-Gaussian which shows that the model provides a reliable uncertainty range.  Figure 4. Comparison of the estimated bird densities (black line, 10-90 quantiles uncertainty range in grey) and the bird densities (red dots) observed using dedicated bird radars at two locations in (a) Herzeele, France (50 • 53'05.6"N 2 • 32'40.9"E) and (b) Sempach, Switzerland (47 • 07'41.0"N 8 • 11'32.5"E). Note that because of the power transformation, model uncertainties are larger when the migration intensity is high. It is therefore critical to account for the uncertainty ranges (light gray) when comparing the interpolation results with the bird radars observations (red dots).

Application to bird migration mapping 228
The main outcome of our model is to estimate bird densities at any time and location 229 within the domain of interest. This is illustrated by the estimation of bird densities time 230 series at specific locations (e.g. Figure 4), and by the generation of bird densities maps 231 at different time steps ( Figure 5).  Figure 5. Maps of bird densities estimation every hour of a single night (4/5 October). Civil sunset and sunrise limits are visible on the first and last snapshots. The highest bird densities are in the corridor from Northern Germany to southwestern France. Rain limits migration in Southern Poland, Czech Republic and Southern Germany. The sunrise and sunset fronts are visible at 18:00 and 05:00 with lower densities close to the fronts. A rain cell above Poland blocked migration on the Eastern part of the domain. In contrast, a clear pathway is visible from Northern Germany through to Southwestern France.
While the estimation represents the most likely value of bird density at each node of 233 the grid (e.g. Figure 5), a simulation generate several realizations (i.e. equiprobable 234 outcomes of the migration process) that reproduce the space-time patterns of migration 235 (e.g. Figure 6). As a consequence, only realizations are able to reproduce number during 236 peak migration as noted when comparing the colorscale of Figure 5 and Figure 6. This 237 becomes important when assessing for instance the total number of birds aloft.  advantageous if the correlation between bird densities and these covariates is stronger 281 than the spatio-temporal correlation of the nearby radars.

282
In addition to quantifying bird migration at high resolution, we can also deduce the 283 spatio-temporal scales at which migration is happening from the covariance function of 284 the model ( Figure  acquisition and for spread of warning systems on peak migration events.

291
As a proof of concept, we used three weeks of bird densities data available on the 292 ENRAM data repository (see Data Accessibility). As more data from weather radars 293 become available, our analyses can be extended to year-round estimations of migration 294 intensity at the continental scale, in Europe and in North America. We also importantly 295 pre-processed the bird density data, i.e. restricted our model to nocturnal movements 296 and applied a strict manual data cleaning. This is because the bird density data  Supporting Information S1 Data pre-processing 347 Supporting Information S1 describes the full procedure applied to manually clean the raw time series of bird densities.

351
1. The data of 11 radars are discarded because their quality was deemed insufficient by visual inspection. The reasons 352 for this poor quality are various: S-band radar type, altitude cut, poor processing or large gaps. The same radars 353 were removed in . In addition, we also excluded radars of Bulgaria and Portugal (4 radars) from 354 this study because of their geographic isolation and the necessity of spatial coherence in the methodology presented. 355 2. If rain is present at any altitude bin, the full vertical profile was discarded (blue rectangle in Figure S1-1). A 356 dedicated MATLAB GUI was used to visualise the data and manually set bird densities to "not-a-number" in such 357 cases.

358
3. Zones of high bird densities are sometimes incorrectly eliminated in the raw data (red rectangle in Figure S1-1). To 359 address this,  excluded problematic time or height ranges from the data. Here, in order to keep 360 as much data as possible, the data have been manually edited to replace erroneous data either with "not-a-number", 361 or by cubic interpolation using the dedicated MATLAB GUI. 4. Due to ground scattering, the lowest altitude layers are sometimes contaminated by errors, or excluded by the initial 363 automatic cleaning procedure. This is solved by copying the first layer without error into to the lowest ones. 6. Finally, the data recorded during daytime are excluded. Daytime is defined at each radar by the civil dawn and dusk 367 (sun 6 • below horizon).

368
The resulting cleaned vertical-integrated time series of nocturnal bird densities are displayed in Figure S1-1d.
where b i are the coefficients of the polynomial.

382
Both the amplitude and the residual are modelled as stationary processes which can be described with a covariance function (also called auto-covariance). Let us denote the generic standardized random variableX for either the standardized amplitudeÃ or the standardized residualR. The covariance CX is a positive definite function depending only on the lag-distance (∆s, ∆t). In our model, we use covariance functions of Gneiting type (Gneiting, 2002) CX (∆s, ∆t) = cov X (s, t),X (s + ∆s, t + ∆t) The transformed variable Z (s 0 , t 0 ) p is reconstructed by combining the deterministic parts (trend and curve) with the kriging estimation of the amplitude and residual as in Eqn. 1. Because A and R are normally distributed, Z p is also 5/20 normally distributed and its mean µ p Z and variance σ 2 Z p are sufficient to describe its distribution with -9) and as A and R are independent, (S2-10)

S2.1 Probability distribution function of bird densities 402
Because of the power-transform, the probability distribution function (pdf) of Z, denoted by f Z (z), is non-trivial.

403
However, the quantiles of this pdf can be derived analytically from the quantiles of the pdf of Z p as detailed hereafter. 404 We introduce the normally distributed variable X = Z p with a pdf where µ Z p and σ 2 Z p are the mean and variance of Z p derived from Eqn. S2-9 and S2-10.

406
One possible way to compute the pdf of Z p consists in computing the mean of several functions of this random 407 variable. For any given measurable function g, Using the change of variable z = x 1/p , which leads to dx = pZ p−1 dz, Eqn. S2-12 becomes This equation allows identifying the pdf of Z as This last equation provides the analytical pdf of bird densities Z, as the pdf f X (x) = f Z p (z p ) is fully known The probability distribution function of Z is non-symmetric and skewed, and therefore cannot be conveniently described 414 with this expected value and variance. Instead, we use the quantile function Q Z (ρ; s 0 , t 0 ), which returns the bird density 415 value z corresponding to a given quantile ρ The quantile function allows to describe Z because the quantile value ρ is preserved through power transform.

417
Therefore, the quantile function of Z is computed with where F Z p (Z p ) is the cumulative distribution function of Z p (s 0 , t 0 ).

S3 Model parametrisation 420
Supporting Information S3 presents the method of parametrisation and discusses the meaning of model parameters in 421 terms of bird migration.

422
Power transform

423
The value of power transformation p is inferred by maximizing the Kolmogorov-Smirnov criterion of the p-transformed 424 observation data Z(s, t) p . The Kolmogorov-Smirnov test (Massey, 1951) is testing the hypothesis that data Z(s, t) p are 425 normally distributed. The optimal power transformation parameter was found for p = 1/7.4 and the resulting Z p 426 histogram is illustrated in Figure S3-1 together with the initial data Z. The fitted distribution shows that bird densities is highly skewed: the lowest 10% are below 1 bird/km 2 while the 428 upper 10% are above 50 bird/km 2 with density up to 500 bird/km 2 . A power transformation on such skewed data creates 429 significant non-linear effects in the back-transformation. For instance, the symmetric uncertainty of an estimated value in 430 the transformed space (quantified by the variance of estimation) will become highly skewed in the original space. the initial trend is fitted to the average of the transformed bird densities for each radar, (2) the amplitudes are computed 447 from the de-trended data, and (3) the parameters of the curve are derived by fitting a polynomial on the data corrected 448 from the trend and amplitude effects. After convergence of the algorithm, the resulting misfit value becomes the value of 449 the residual. The resulting planar trend is shown in Figure S3-2a together with the average transformed bird densities of 450 each radar. The trend displays a strong North-South gradient, which can be explained by the larger migration activity in 451 southern Europe during the study period. A 2-dimensional planar trend was initially tested in order to accommodate the 452 northeast-southwest flyway. However, this more complex model did not significantly improved the fit to data, and has 453 therefore been discarded. The de-trended values illustrated in Figure S3-2b are more stationary with the exception of 454 Finland and Sweden.  also noted this difference between both countries, but excluded the fact that 455 this difference is due to errors in the data since the southern Swedish radar shows consistent amounts of migratory 456 movements with a neighbouring German site. Figure  Transformed bird densities  The calibrated covariance functions provide information about the degree of spatial and temporal correlation of the bird 478 migration process. The spatial covariance of the amplitude ( Figure S3-5a) shows that the nocturnal bird densities are 479 well-correlated for locations separated by less than 500 km, and completely uncorrelated for more than 1500 km. The 480 temporal covariance has an asymptotic behavior and never decreases under 0.2 ( Figure S3-5b). This non-zero sill can be 481 due to either a remaining temporal non-stationarity in the dataset, or systematic errors in radar observations (caused by 482 e.g. different types of technology or local topography affecting migration). Note that since the covariance is evaluated 483 only on a discrete 1-day lag-distance, the shape of the covariance between 0 and 1 is artificially created to fit the Gneiting 484 function. Overall, the temporal correlation of the amplitude is weak with only 40% for the covariance for lag-1. It is important to recall that since the weather radars are relatively well-spread (Figure 1b), the spatial covariance of both the 486 amplitude and residual is poorly constrained for lag-distances below 100 km, and consequently the importance of the distribution is close to zero which indicates that the estimation is unbiased (i.e. in average, the estimation is neither 495 underestimating (mean below 1) nor overestimating (mean above 1)). However, its variance is below 1, which indicates a 496 slight overestimation of the uncertainty range (i.e. in average, the uncertainty ranges are too wide). Next, the normalized kriging error is assessed for each radar ( Figure S4-2). The resulting distributions indicate that 498 the goodness of the estimation is different for each radar. In Figure S4-  The cross-validation is further illustrated in Figure S4-4 for a specific radar located in Boostedt, Germany 501 (54 • 00'16"N, 10 • 02'49"E) indicated with gold circle in Figure S4-3. For this radar, the general pattern of the signal is well 502 estimated for both the nightly amplitude and the intra-night variation. Bird densities are underestimated during the peak 503 migration occurring at 4th and 5th of October, is but the estimated value remains within the uncertainty bounds. The main block of the website is a map with a standard interactive visualization allowing for zoom and pan. On top of 511 this map, three layers can be displayed:

512
• Layer 1 corresponds to bird densities displayed in a log-color scale. This layer can display either the estimation map, 513 or a single simulation map by using the drop-down menu (1a).

514
• Layer 2 corresponds to the rain (rainy areas are in light blue), which can be hidden/displayed with a checkbox (1b) 515 • Layer 3 corresponds to the bird flight speed and direction, displayed by black arrows. The checkbox allows to 516 17/20 display/hide this layer (1c). The last item on the top-right menu is the link menu (1d).

517
Block 2: time series 518 The second block (hidden by default on the website) shows three time series, each in a different tab (2a):

519
• Densities profile shows the bird densities [bird/km 2 ] at a specific location.

520
• Sum profile shows the total number of bird [bird] over an area.

521
• MRT profile shows the mean traffic rate (MTR) [bird/km/hr] perpendicular to a transect.