Scientific Papers

A three-dimensional model of terrain-induced updrafts for movement ecology studies | Movement Ecology

Description of Image

In this section, we describe the methodology behind each of the steps used to develop, validate, and apply the improved updraft model. Specifically, “Large-eddy simulation methodology” section describes the LES model and details used to obtain the terrain-resolving turbulent flow field. Next, in “Model development” section we describe the updraft model development, separated in two parts. First, in “Parametric study using synthetic terrains” section, we perform a parametric study using simple gaussian hill geometries to obtain a height adjustment \(f_h\) based on the wind speed at a reference height \(V_\text {ref.}\). In “Model tuning with real terrains” we use LES solutions over real terrains of different levels of complexity to further tune the model and develop adjusting factors based on a measure of terrain channeling and sheltering effects \(f_{Sx}\), and terrain complexity \(f_{tc}\). In “Model validation with field observations” section, we describe model validation using field measurements of vertical velocities over complex terrain. Finally, in “Model illustration with movement model simulations” section we describe the application of a movement model for golden eagles, using both the original and the adjusted model.

Large-eddy simulation methodology

We use the Simulator for Wind Farm Application (SOWFA), developed at the National Renewable Energy Laboratory. SOWFA is an open-source turbulence-resolving LES tool, focused on wind energy applications. Some of the capabilities of SOWFA are its ability to represent full diurnal cycles driven by large-scale mesoscale conditions, different canonical stability states, realistic atmospheric features (e.g., low-level jets, weather fronts, or gravity waves), and features like wind turbines and complex terrain. SOWFA has been extensively used and validated in works within the wind energy community [9, 13, 14, 20, 24, 27, 39]. Its capability of including complex terrain and capturing terrain-generated turbulence has also been investigated [22, 47]. An LES approach is used in SOWFA, meaning that turbulence is resolved in the computational fluid dynamics simulation down to the grid scale, and smaller, sub-grid-scale turbulence is modeled. The reader is referred to the work of Churchfield et al. [10] for more details on SOWFA. The effect of the rough planetary surface is modeled by the Schumann surface stress boundary condition, and the aerodynamic surface roughness is kept constant at 0.15 m.

In the numerical studies conducted, we take a two-step approach in coupling complex terrain to atmospheric turbulence. First, we execute a case where the computational domain has periodic boundary conditions on its lateral boundaries and a flat bottom boundary. This step is known in the literature as a precursor simulation. The goal of the precursor is to spin up the background atmospheric turbulence. When turbulence is developed—identified by the convergence of profiles of turbulence metrics and a resolved energy cascade that follows Kolmogorov’s law [29]—information at the lateral boundary planes is saved. In this step, we assume horizontally homogeneous turbulence by the use of cyclic boundary conditions on the lateral boundaries. Next, a second simulation is set up where the bottom boundary is conformed to a digital elevation model of the terrain geometry and the boundary data saved from the precursor are used as boundary conditions at the inflow planes. The real terrains considered in this work are obtained using NASA’s Shuttle Radar Topography Mission [64] digital elevation model, with resolution of approximately 30 m.

The domain size and grid resolution are responsible for the upper and lower bounds of the scales captured. Due to horizontal homogeneity, the largest scales are of the order of the domain size. The lowest scales resolved are of the order of 4–5 times the grid size [45]. Simulation of atmospheric boundary layers is a multiscale problem and can quickly become computationally expensive if a wide range of scales is needed. Here, we select domain sizes, resolution, and simulation time such that the temporal and spatial scales of interest are properly captured. Note that throughout this work, the LES results are presented in a time-averaged fashion.

Model development

Parametric study using synthetic terrains

We developed the empirical model using a parametric study of idealized Gaussian hill geometries. The conditions are canonical neutrally stratified boundary layers, with a wind speed of 8 m/s at a reference height of 80 m AGL. This reference height was used throughout this study because of the availability of measurements from reanalysis datasets at this height and because it is a typical measurement height on tall meteorological masts where horizontal wind speed data are often collected for wind-energy applications. The numerical setup includes a \(3\times 3\times 1\) km domain with uniform 10-m grid resolution. A capping inversion is enforced by a 10 K stable stratification between 750 and 850 m AGL, effectively limiting the growth of the boundary layer. The domain extents and resolution used are typical values from the atmospheric and wind energy LES community (see [9, 14, 39] and references therein). In this parametric study, we analyze four hills with varying steepness levels, shown in Fig. 1. We investigate the effect of the wind angle \(\alpha -\beta\) by analyzing two angles, 0 and 45 degrees. In the 0 degrees case, hereafter called orthogonal flow, the winds are orthogonal with the ridgetop and, by Eq. (1) the resulting updrafts are the largest possible. In 45 degrees flow, hereafter called oblique flow, the winds approach the slopes at a 45 degree angle. In this study, the hills were designed so that the two of higher steepness are expected to produce separated flow under an 8 m/s wind. The flow separation in such cases induces significantly higher turbulence levels.

Fig. 1
figure 1

The four Gaussian hill geometries used to perform the parametric study. a Elevation; b orographic slope

In order to keep the same incoming flow to both the orthogonal and oblique cases, instead of executing two different wind directions, we opt to manipulate the geometry. Our background wind field was simulated for incoming southwest winds. That means that for our oblique case, the hills shown in Fig. 1 ran from north to south; whereas for the orthogonal flow case, the hills were placed diagonally on the domain, running from the northwest to the southeast. Figure 2 shows the nondimensional time-averaged updraft field obtained from LES at selected heights and estimations from the \(\textsc {BO04}\) model for the orthogonal flow case. In the figure, the top of the hill is exactly on the diagonal of the domain (more clearly visible in Fig. 2i–l), and the arrow notes the direction of the wind. It becomes clear that steep hills induce flow separation and high values of vertical velocity. The imprint of the hill is less evident as the height increases; that is, the updraft field becomes less strong as the height AGL increases, as expected.

Fig. 2
figure 2

Large-eddy simulation solutions of orthogonal flow compared to the \(\textsc {BO04}\) model around hills placed diagonally, with the hilltop placed exactly on the diagonal of the domain shown. Each row corresponds to a height AGL, as indicated by the labels on the left. The black arrow at the bottom row represents the direction of the flow. Each column corresponds to a hill, reproduced at the top row (to scale). Nondimensional vertical velocity from ad LES at 80 m, eh LES at 180 m, and il \(\textsc {BO04}\) model. Note that \(\textsc {BO04}\) model is height-independent

This first study also reveals that positive updraft can be observed on the leeward side of the hill. Asymmetry in the windward and leeward sides of hills is impossible to be represented by the \(\textsc {BO04}\) model. Even with LES, the turbulent flow on the leeward side of steep terrain features can be challenging to model. In the present study, this secondary updraft zone has not been considered.

Analogous LES results to those presented in Fig. 2 were developed for selected AGL values, roughly covering a typical land-based wind turbine rotor-swept zone: 40, 80, 120, and 180 m AGL. For each scenario at each height, the \(\textsc {BO04}\) model results are compared with the LES on a scatter plot. A clear linear relationship is observed for the two less-steep geometries, whereas no clear relationship is found for the steeper hills. We found that the points not following a linear relationship in those cases are related to the separated flow region and strong turbulence. When we consider only the positive vertical velocity values, however, the linear relationship is recovered for all scenarios. Therefore, for each combination of flow angle, hill steepness, and height, we obtain the slope of the linear regression fit of the scatter data. The linear regression fit slope related to each case is shown in Fig. 3. The variation of the slope versus both the height AGL and the cosine of the maximum orographic slope value (see Fig. 1) is presented. Considering that at the limit, a flat terrain produces no orographic-induced vertical velocity, the variation of the slope of the linear regression fit with respect to height can be approximated by a second-order polynomial, while its variation with the cosine of maximum hill slope can be approximated by an exponential relationship.

Fig. 3
figure 3

Slopes of the linear fit on the relationship of the positive vertical velocity for each of the hills and flow direction investigated. Panels a and b show different views of the 3-D data shown in c. The gray surface in c is the surface given by Eq. 2

We fit these points shown in Fig. 3 with a surface that varies quadratically with the height h, and exponentially with the cosine of the slope \(\theta\), in an expression of the form

$$\begin{aligned} f_h = \left( ah^2+ bh + c \right) \cdot \left( d^{-\cos \theta + e} \right) + f \end{aligned}$$


where \(f_h\) is an adjustment factor to be used on top of the \(\textsc {BO04}\) model. We approximate the surface fit, which yields \(a=0.00004\) m\(^{-2}\), \(b=0.0028\) m\(^{-1}\), \(c=0.8\), \(d=0.35\), \(e=0.095\), and \(f=-0.09\). While this adjustment is mathematically valid for all heights, given the fit is performed at heights within the rotor-swept zone, care should be taken when using it for heights much lower than 30 m and much higher than about 200 m. An illustration of the fitted surface is given in Fig. 3c.

A straightforward way to visualize the adjusted results is to look at their cross sections as opposed to the whole domain as shown previously in Fig. 2. The cases investigated, even though they are fully three-dimensional, can be averaged in the ridge-wise direction and presented as two-dimensional profiles, shown in Fig. 4 for selected heights. The figure compares LES results with the height-adjusted model (Eq. 2) and the original model. Naturally, the height-adjusted results will be a better match on these results since the correction was developed with the same data. It is of interest, however, to note how nonzero updraft regions can extend far upstream of the geometry on the LES results. The extent length seems to be related to the steepness of the hill, where steeper hills have updraft regions extending further upstream.

Fig. 4
figure 4

Curves of nondimensional updraft velocity averaged along the ridge line: ad orthogonal flow; eh oblique flow. While the height adjustment \(W_{0h} = W_0 / f_h\) improves the magnitude with respect to the LES solutions, it does not account for nonzero updrafts upstream of the hill

Taking a closer look at the upstream extent up vertical velocities in each case, we can see that the updraft is not only a function of the slope present directly below the point considered but also contains information from surrounding orographic features. At 40 m AGL, positive updrafts have been observed between 100 and 250 m upstream from the start of the hill (for the Gaussian geometries under investigation we define the start of the hill as 200 m upstream the hilltop). At higher heights, such as 160 m AGL, such distance can be greater than 600 m. At 160 m, the nondimensional vertical velocity values are less than half of those at 40 m. Therefore, for every location, we need to incorporate information about downstream terrain features. A useful terrain metric in this context is that of wind shelter or exposure. Developed by Winstral et al. [69] and Winstral and Marks [68] with focus on snow accumulation and melt, a terrain exposure/sheltering parameter, Sx, was initially formulated based on data from the Green Lakes Valley in Colorado and later validated with data from an exposed site and a sheltered site in the Reynolds Mountain East watershed, located in Idaho. Sx is based on the maximum upwind slopes found between the point of interest and a pie-shaped sector extending upwind of the point in the direction of interest:

$$\begin{aligned} {\widetilde{Sx}}_{{\widetilde{A}},\text {dmax}} \left( x_i,y_i\right) = \max \left[ \tan \left( \frac{ z\left( x_v,y_v\right) – z\left( x_i,y_i\right) }{\sqrt{\left( x_v-x_i\right) ^2 + \left( y_v-y_i\right) ^2}} \right) \right] \end{aligned}$$


where z is the elevation map from a DEM, \({\widetilde{A}}\) is the azimuth of the search direction, \((x_i,y_i)\) are the coordinates of the cell of interest, and \((x_i,y_i)\) are the set of all cell coordinates located along the search vector defined by \((x_i,y_i)\), \({\widetilde{A}}\), and \(\text {dmax}\). The extent of the sector (\(\text {dmax}\)) is a tuning parameter. A negative \({\widetilde{Sx}}\) value indicates an exposed location, whereas a positive value indicates a sheltered location. As recommended by the authors of the model, we average \({\widetilde{Sx}}\) across an upwind window of directions, resulting in a measure that is more robust to both natural and systematic deviations from the recorded data, Sx. The mean maximum upwind slope parameter is then given by

$$\begin{aligned} Sx_{A,\text {dmax}} \left. \left( x_i,y_i\right) \right| _{A_1}^{A_2} = \frac{1}{n_v} \sum _{A=A_1}^{A_2} {\widetilde{Sx}}_{A,\text {dmax}} \left( x_i,y_i\right) \end{aligned}$$


where \(A_1\) and \(A_2\) define the outer limits of the upwind window, A bisects \(A_1\) and \(A_2\), and \(n_v\) is the number of search vectors in the window defined by \(A_1\) and \(A_2\).

The Sx quantity is used differently in this work. Based on its original definition, it gives information about upstream features, whereas our interest is primarily on downstream features. For that, we flip the wind direction that goes into the Sx calculation so that \(A = (\text {wdir}+180)\%360\) with wdir being the wind direction given in degrees and the % operator denoting rest of division (modulo operator). The consequence for that is that we have a positive angle on the windward side and negative on the leeward side. We use a 30-degree window for all Sx derivations and \({\widetilde{Sx}}\) calculated in 5 degree increments. Given the range of values observed previously, \(\text {dmax}\) is selected to be 500 m.

Model tuning with real terrains

Next, we develop additional empirical corrections to account for flow effects over real terrain. Real terrain rarely contains isolated features like the idealized hills analyzed above, and thus complex flow patterns are expected to occur, including some wind channeling effects. We emphasize that in the development of the adjustment factors, we do not attempt to model the complex turbulent flow on the leeward side of the terrain; rather, our focus is on the upwind side.

We selected two regions containing different types of terrain features, shown in Fig. 5. The first is relatively mildly sloping unstructured terrain near the Top of the World wind farm in Wyoming (referred to as Wyoming region). In this region, flow separation would not be expected except under high wind conditions at a few isolated locations. The second region is within the Valley and Ridge physiography of the Appalachian Mountains in central Pennsylvania, known for its long linear ridges (referred to as Pennsylvania region). The domain in this case includes a steep ridge similar to the idealized one above and represents a more extreme terrain type. In this region, strong separated flow is expected to occur. The choice of these sites is primarily due to the contrast in terrain, but also the Wyoming region has a high density of golden eagles and has had documented fatalities at several wind farms, whereas the Valley and Ridge region of Pennsylvania is a known migratory path for golden eagles.

Fig. 5
figure 5

Regions considered for model tuning. a, d 10 m digital elevation map, b, e corresponding slope, and c, f aspect. ac Pennsylvania region; df Wyoming region. Note the different scales

In this part of the study, the numerical domain extents were increased to \(5\times 5\times 2\) km, with a uniform 10 m grid resolution. The taller domain allows the boundary layer to deform according to the underlying terrain without introducing numerical artifacts at upper boundary. The choice of wind speeds and wind directions to be simulated as canonical neutral boundary layers followed typical conditions observed at the site of interest during the years of 2017 and 2018. Historical data were obtained through the Wind Integration National Dataset (WIND) Toolkit [16]. Wind roses of the region reveal that the wind is predominantly from the west during the summer and spring seasons. For the analysis, two wind speeds were chosen. One of 8 m/s representing the low wind speeds observed during the day at those seasons, which corresponds to about 78% of the winds encountered at the site, and another wind speed of 15 m/s, representing an extreme. A wind speed of 15 m/s was only observed during 4% of the time at the region, but represents a challenging, yet realistic scenario that is important to account for in the model.

As mentioned in “Parametric study using synthetic terrains” section, we use the modified Sx parameter to account for downstream features. An adjustment for downstream features based on the tangent of the exposure quantity Sx, considering the flipped wind direction, is given as

$$\begin{aligned} f_{Sx} = 1 + \tan \left[ Sx\left( A=(\text {wdir}+180)\%360, \,\text {dmax}=500 \text { m} \right) \right] \end{aligned}$$


where % denotes the modulo operator, and the wind direction is given using typical wind direction convention where 0 degrees represents wind coming from the north. Note the flipping of wind direction gives us the correct sign of the factor \(f_{Sx}\) that allows straightforward multiplication of all the individual factors.

Upon exploration of the LES results, it becomes clear that the impact of small terrain features is less defined at higher heights. That is, the updraft becomes less sensitive to isolated terrain features and instead resembles an integration of the effects of nearby features. Based on observations of the LES results at different heights, we experimented with a direct convolution of terrain metrics with a Gaussian function. Using the slope and aspect quantities blurred by Gaussian kernels resulted in blurred versions of the \(\textsc {BO04}\) model updraft field, which more closely matched that obtained using LES. We then devise another aspect of our improved model, which is to use metrics convoluted with Gaussian functions. The original updraft coefficient given by Eq. (1) is then computed using the filtered metrics. The Gaussian blurring kernels have an empirically determined standard deviation \(\sigma\) that varies linearly with height in the form of

$$\begin{aligned} \sigma = \min (0.8h + 16, \, 300) \end{aligned}$$


where h is the height AGL in meters. The minimum guarantees a maximum kernel size at much higher heights. Considering the Wyoming region, the blurring step is illustrated in Fig. 6, showing the dimensional vertical velocity field obtained using LES at different heights, the original model, their counterpart Sx– and height-adjusted, and the proposed blurring approach. As noted earlier, the original model is sensitive to small-scale terrain features, resulting not only in localized overestimation and underestimation of the updrafts, but also on a field that is too well defined with respect to the underlying terrain and not realistic.

Fig. 6
figure 6

Illustration of the intermediate steps of the model, showing the dimensional vertical velocity field around the Wyoming region (Fig. 5d–f). a, e LES; b, f original \(\textsc {BO04}\); c, g proposed model with Sx and height adjustment; d, h proposed model with height and blur adjustment. Each row illustrates one height AGL, indicated on the left

The blurred adjustment results in a flow field that more closely matches the LES model, although it generally underestimates the values from the LES. Therefore, a final adjustment was developed related to terrain complexity. The Pennsylvania region, with its significant spread in elevation-induced flow separation and high levels of turbulence as well as vertical velocities that exceeded 1.5 m/s at the top of the ridge. The Wyoming region with mildly complex terrain, on the other hand, has mostly attached flows. In order to account for the steep hills where the vertical velocity is significantly higher, we introduce a terrain complexity factor. Many different measures have been used to describe the complexity of peaks, valleys, and ridges—for example, variance of elevations, autocorrelation of elevation, relief, contour density, rugosity, curvature, slope change, aspect change, etc. A good discussion on the topic can be found in the work of Huaxing [23]. For the purposes of this work, we are interested in a metric that is simple and is localized. By localized, we mean that in a large region with heterogeneous terrain, complex regions will be assigned a high value for its local features, while flat regions will be assigned a low value. We derive a simple metric that is based on the mean elevation surrounding each DEM cell and further scaled by the relief in the same local region. A terrain complexity factor is given by

$$\begin{aligned} \text {tc} = \frac{ \hspace{0.83328pt}\overline{\hspace{-0.83328pt}\langle z \rangle \hspace{-0.83328pt}}\hspace{0.83328pt}- \langle z \rangle _\text {min} }{ \langle z\rangle _\text {max} – \langle z\rangle _\text {min} } \end{aligned}$$


where the angled bracket symbols refer to quantities within a \(500 \times 500\) m area, centered around the cell of interest, and the overbar represents a mean quantity. The expression on the denominator is the relief of each region. The terrain complexity is further multiplied by another parameter, \(\Lambda\), which is a linear function of the height AGL, \(\Lambda =h/40\). The normalizing value of 40 m was obtained from regression analysis comparing LES data to the model, and represents a mean value across all cases investigated. Finally, the terrain complexity factor is given as the following:

$$\begin{aligned} f_{tc} = 1 + \Lambda \cdot \text {tc}. \end{aligned}$$


The \(\Lambda\) parameter makes the terrain complexity factor increase with height. It essentially gives a larger weight for high heights without compromising lower heights. This is necessary to counterbalance the second-order polynomial that grows fast as height increases.

Finally, we compound the adjustment factors discussed to come up with an improved model for the orographic updrafts, \(w_{0i}\). The final adjusting factor is

$$\begin{aligned} F = \frac{f_{Sx} \cdot f_{tc}}{f_h}, \end{aligned}$$


and the model reads

$$\begin{aligned} W_{0i} = F \cdot W’_0, \qquad w_{0i} = V_\text {ref.} W_{0i} \end{aligned}$$


where \(V_\text {ref.}\) is the horizontal wind speed at a reference height, 80 m, and the \(W’_0\) is the original coefficient (Eq. (1)) computed using slope and aspect blurred by the Gaussian kernel with standard deviation given by Eq. (6). As mentioned, the idea is that horizontal wind speeds at this reference height can be obtained either by field experiments (where 80 m is a height where instruments are typically deployed) or by analysis models, thus removing the need to obtain the wind speeds at the height of interest.

Statistical comparison of the resulting model with LES results over the Wyoming and Pennsylvania terrains is included in “Model comparison to LES using real terrains” section, computed using the differences between modeled updraft and LES updraft in regions where the modeled vertical velocities are positive.

To briefly summarize the previous sections, the model development process is comprised of the following steps:

  • Replace the multiplying wind speed \(w_0/W_0\) of the original \(\textsc {BO04}\) model by a reference wind speed obtained at a fixed, reference height, called \(V_\text {ref.}\);

  • Use a smoothed version of the digital elevation model (by means of a blurring kernel) to compute the aspect and slope quantities, yielding a quantity we term \(W’_0\);

  • Apply empirical adjustments \(f_i\) as multiplying factors on top of a non-dimensional \(W’_0\) quantity to determine the dimensional vertical velocity \(W_{0i}\).

A summary of inputs and outputs of the model is given in Fig. 7.

Fig. 7
figure 7

Inputs and outputs for the proposed model. Details regarding the methods behind the proposed model block are given in “Model development” section

Model validation with field observations

As discussed in the previous sections, the model was developed from LES results, and while the SOWFA code used in this work has also been used and validated in the wind energy and atmospheric modeling communities (see references given in “Large-eddy simulation methodology” section), vertical velocity is not a quantity often investigated in detail. Therefore we sought to compare the model results with field data. Our particular interest is in measured vertical velocities at heights spanning a few hundred meters AGL, corresponding to the rotor-swept zone of typical land-based wind turbines.

Several candidate field studies were considered, including Bolund hill, Askervein hill, Perdigao mountain, and Alaiz mountain. Bolund hill is a small peninsular feature that is 12 m high, 130 m long, and 75 m wide located near the city of Roskilde in Denmark. Its geometrical shape induces complex three-dimensional flow. The experimental campaign [3], however, contains masts that are only 10 m tall. Askervein hill [60], located in Scotland, has a nearly Gaussian shape with mild slopes. The campaign, performed in 1985, has masts as tall as 50 m, but most of the data are taken using significantly shorter masts. Most importantly, however, no vertical component of the velocity is given, only statistics of such component. Perdigao mountain was the site of a thorough field campaign, conducted in 2018 [18]. Perdigao is a double-ridge case, located in Portugal, selected for its challenging interaction of one ridge’s wake on the other. The campaign had tens of masts spread across the region, and three of them contains vertical velocity data up to 100 m AGL. We looked at several months worth of data and noted very strong diurnal cycle characteristics imprinted on the vertical velocity time-history. That means that the region was subject to strong convective effects during daytime and stable boundary layers at nighttime. Both the model developed in this work and the original one are not designed for non-neutral conditions, so they are not applicable to the conditions observed during the campaign.

For model validation we selected the Alaiz mountain, Spain field campaign. The Alaiz region is presented in Fig. 8. It is a mountain where flows from every direction are subject to complex orographic features. In particular, southerly flows experience complex mountainous regions with small ridges before a steep incline, and northerly flow is subject to a steep incline. A DEM of the Alaiz mountain region with 2-m resolution was provided by authors of the experiment.

Fig. 8
figure 8

Overview of the Alaiz mountain region. a General surroundings with the MP5 mast location shown by a marker. Zoom around the MP5 mast, b elevation; c slope

We selected hilltop reference mast MP5 that contained instruments and measurements of the vertical velocity at heights up to 118 m AGL was selected (see Fig. 8). Three-dimensional wind speed data collected by sonic anemometers were available between the months of July 2017 and July 2019. We focused on the period spanning July to Dec 2017. The three components of the wind speed are available in 10-min interval means, including the maximum, minimum, and standard deviation of several variables within each 10-min interval. For this validation exercise, we focus on the comparison of vertical velocity component only.

No strong convective updrafts were observed in the data, suggesting the orographic lift conditions for which the model was developed. We assume that any thermally induced vertical disturbances present are averaged out in the 10 min windows. Buoyant disturbances are naturally present in the atmospheric boundary layer due to the heat flux from the ground surface. Given vegetation and other ground irregularities present at the site, localized temperature differences are generally random. Large pockets of positive and negative vertical velocity develop. When coupled with a non-zero horizontal wind speed, the disturbances vary in time and space in an unpredictable manner. Our assumption is rooted in numerical experiments of unstable boundary layers and holds true for canonical unstable simulations with uniform heat flux for averaging periods as short as 5 min. Further analysis of the temporal window of the available data was performed, indicating thermal updrafts had no impact on the data.

Comparison between the field data and both the original and proposed models were done for every 10-min-mean datapoint. For each datapoint, we obtain the wind direction and the horizontal wind direction, which are directly used in both the \(\textsc {BO04}\) and the proposed model. Time-series of the vertical velocity is then obtained for both models and simple statistics metrics are obtained. Some times with missing field data indicate that the vertical velocity was not available at those times, but horizontal velocity was, and thus the modeled values could be determined. For details on the instruments and data processing and quality, the reader is referred to the experiment’s original references [37, 52]. Results of the validation with field observations are presented in “Model validation with field data” section.

Model illustration with movement model simulations

In the final part of the study, our goal was to determine the extent to which simulated movement tracks of a soaring bird are dependent on the updraft model chosen, and whether the proposed updraft model yields more realistic results at higher flight altitudes where orographic lift weakens. We use two sites with different terrains and wind flow regimes and estimate the updraft using the \(\textsc {BO04}\) model and the proposed model at two different heights. The flight tracks are simulated with a heuristic individual-based movement model (IBMM) within the Stochastic Soaring Raptor Simulator framework [50].

In the IBMM, the simulated bird reduces its energy expenditure by selecting a movement path based solely on local updraft availability and intensity above a threshold value. This is a purposely simplistic approach that does not consider the possibility that the bird makes movement decisions based on the overall structure of the landscape/updraft field within view; rather, it assumes that the bird responds rapidly to updraft conditions experienced during flight in its immediate vicinity [32]. The model has only two parameters: a specified direction of motion (\(\phi\)) and the species-dependent threshold updraft velocity (\(w_\text {thr}\)) that will sustain soaring/gliding flight. In this example we use 0.85 m/s for \(w_\text {thr}\), based on a minimum sink speed calculation [43] using data from Lish et al. [35] for golden eagles. The model track simulation goes as follows: the individual moves in constant-length steps (here, we use 30 m) across the domain from a specified starting location. At each step it assesses updraft availability and intensity along a 60 degree arc one step ahead of its current position (at \(\phi\), \(\phi \pm 15\) deg, \(\phi \pm 30\) deg). The modeled updraft velocity field is bilinearly interpolated to these locations. If the interpolated updraft velocity at one or more of these five locations exceeds \(w_\text {thr}\), the individual moves to the location of maximum updraft. If not, it selects its next location randomly (i.e., the movement path becomes a directed random walk). The model is run for 1000 simulated tracks across the model domain from various starting positions along a model boundary, and a presence map is generated using a Gaussian smoothing filter over the 1000 tracks.

The two sites selected in this example are in central Pennsylvania within the Valley and Ridge physiography of the Appalachian Mountains, and near Altamont Pass in the Diablo Range of California. The central Pennsylvania site, a different site than one previously used to develop the model, is dominated topographically by distinct and steep northeast-to-southwest trending ridges, whereas the Altamont Pass region is less well-structured, with a mix of mildly sloping and steeper terrain, including many small-scale features that modify the windfield. Using the WIND Toolkit, mean wind conditions are obtained for the fall season (September–November) for the Appalachian Mountains region, and for the spring season (March–May) for the Altamont Pass region. We obtain the mean wind conditions at two different heights: 80 m and 180 m above ground level; 80 m is the reference height of the proposed model, and 180 m is the highest height where atmospheric data are available through the WIND Toolkit. The seasons were selected because they represent migratory seasons with significant movement activity in each respective region. An analysis of WIND Toolkit wind roses during daytime hours for each season was performed with the goal of selecting both typical and atypical (although realistic) conditions for each site. Wind roses for each site at each height are shown in Fig. 9.

Fig. 9
figure 9

Wind roses for daytime hours at two different heights for a fall season at the Appalachian mountains in central Pennsylvania; b spring season at the Altamont Pass region in California

A single wind direction was selected for each site, corresponding to the prevailing west-northwest direction for the Appalachian region, and west-southwest for the Altamont Pass region. For both locations we obtained the mean and the standard deviation of the wind speed at the prevailing wind direction. We devised three wind conditions, based on the mean value plus or minus one standard deviation. We refer to these synthetic conditions as “low,” “moderate,” and “high” wind conditions, for each height and for each location. The final values for each condition are shown in Table 1. The updrafts were computed for each wind condition from Table 1 using both \(\textsc {BO04}\) and the proposed model. The IBMM was run for each case, with the updraft model (\(\textsc {BO04}\) vs. proposed) as the only difference between the model runs. Results are presented in “Movement model simulations” section.

Table 1 Typical wind conditions for each site used for the movement model simulation

Description of Image

Source link