Estimation of LiDAR error over complex terrain covered with forest using numerical tools

– Since a few years, a new wind measurement instrument has been competing with standard cup anemometers: the LiDAR. The performances of this instrument over complex terrain are still a matter of debate and this is mainly due to the ﬂow homogeneity assumption made by the instrument. In this work, the error caused by this hypothesis was evaluated with the help of OpenFOAM 1.7, MeteoDyn WT 4.0 and WAsP Engineering for a LiDAR deployed on a complex site covered with dense forest. The assessment of the CFD model ﬁrstly revealed the signiﬁcant impact of both the location and nature of the inlet boundary condition. Despite the presence of terrain complexity within a radius of 340 m around the remote sensor, an averaged error of less than 3% was observed, suggesting that the LiDAR is only aﬀected by topographic variations in the immediate vicinity of the scanned volume.


Introduction
Currently, cup anemometers are the standard in the wind industry due their robustness and available documentation regarding their uncertainty [1].However, they must be mounted on meteorological masts, which represent important investments.Masts are vulnerable to harsh climate and their relocalization is complicated.For these reasons, another type of wind measuring instrument has recently come into play: remote sensing devices.Mostly based on Doppler effect, these devices are more commonly known as SODARs (Sonic Detection and a Corresponding author: e.jeannotte@gmail.comRanging) and LiDARs (Light Detection and Ranging).These instruments emit whether sound or light in the atmosphere and measure the frequency shift of the reflected waves in order to calculate the wind speed.SO-DARs and LiDARs are very interesting options for the industry.For instance, they can measure wind speeds at high altitudes (up to 200 m) and this for many heights simultaneously.Over relatively flat terrains, remote sensing devices showed good performances when compared to cup anemometers [2].Despite these advantages, there is one main drawback associated to them.While retrieving wind speed, wind flow within the sampled volume is assumed to be homogeneous, which entails serious error over complex terrain [3,4].Among proposed solutions [3], it is Article published by EDP Sciences suggested to model the error with the help of numerical models such as WAsP Engineering.In the present study, this approach was revisited using the open source CFD toolbox OpenFOAM.
The main objective of this work was to evaluate the LiDAR error over complex terrain covered with forest with the aid of OpenFOAM in order to compare the results to experimental data and results obtained with the commercial softwares MeteoDyn and WAsP Engineering.However, it was considered appropriate to firstly proceed to an evaluation of the model.

Mathematical model
First, it is important to present the assumptions applied within the context of this study.The flow was considered incompressible, neutrally stratified and steady.Coriolis force was neglected since only the surface layer was considered.In order to consider the effect of forest, the same approach as [5] was implemented, which consists of modeling the forest as a porous medium characterised by a drag coefficient (C d ) and a leaf area distribution (α).Source terms depending on these two variables were added to the RANS equations of motions and turbulence to simulate pressure loss and increase of turbulence respectively.It was chosen to bring these modifications to the standard kturbulence model first for its simplicity, and second for its already proven capability of modeling atmospheric flow [5,6].

Evaluation of the model over complex terrain
The LiDAR used for the present study is located in Anse-à-Valleau (AAV), Gaspe Peninsula on a complex site almost entirely covered with dense forest.Terrain complexity was analysed using IEC 61400-12-2 draft standard considering topography up to a range of 340 m around the remote sensor (see Fig. 5), which corresponds to 8 times the radius of the scanned disc when the instrument shoots at 80 m.It was classified as being moderately (Class 2) to highly (Class 4) complex.
In order to evaluate the performance of the model, data from a measurement campaign undertaken from the 14/12/2007 to 29/12/2009 were used.More precisely, a flow simulation was run for all 10 • directions where experimental data were available, i.e. [270 • -350 • ].The calculated ratio of horizontal velocity between two cup anemometers was compared to experimental measurements for each wind direction.Locations of the masts (AAV3 & AAV6) on which the anemometers are mounted are shown is Figure 2.

Mesh
A mesh aligned with the flow was constructed for each wind direction.A domain sensitivity analysis was performed on wind direction 270 • to determine the optimal size of the domain (see Fig. 2 for dimensions).The background mesh (Δx and Δy) measures ∼35 m.Local refinement was done within radii of 300 m and 150 m away from each mast, dividing the background mesh by a factor of 2 in both directions.In order to correctly capture the gradients near the ground, a vertical expansion ratio of 1.045 was used.This results in cells ∼1 m to ∼160 m high, from lower to upper boundary respectively.
A GCI [7] was conducted using the x-component of velocity at AAV6.The sizes of the considered meshes were 1 525 997, 3 237 951 and 7 153 126 cells.The grid refinement factor, as recommended by [7], was kept ≈1.3.The test revealed a discretization error of 3.86% for the mesh used in the present study.

Boundary conditions
For lateral boundaries, a zero gradient was used for all variables.Ground conditions were prescribed according to [8], who state that a shear stress based on local flow conditions must be applied.Forest creates an unkown displacement height at the lower boundary which makes impossible the direct prescription of velocity (U ) and turbulent quantities (k, ) at the top boundary.In this work, shear stress was fixed through the gradient of velocity by selecting the desired friction velocity (u * ).TKE was assigned a zero gradient.The vertical gradient of dissipation rate was expressed in terms of the turbulent viscosity (ν t ) and the friction velocity, which eliminates all dependence on the vertical position and thus the displacement height.This set of top boundary conditions (Eq. ( 1)), developed by [9], was validated over flat terrain and results agree with [5].
At inlet, neutrally stratified logarithmic wind profiles were imposed for U , k and .These profiles are an approximation for lack of experimental results over complex terrain.In order to limit the effects of such an approach, an artificially smooth transitional zone was added between the inlet and the real topography.That zone starts with 1km of flat terrain followed by 1 km of terrain that slowly introduces topography.The friction velocity was tuned for each direction to match the average experimental velocity at AAV6.The roughness length, neglecting the displacement height, was deduced from fully developed profiles obtained on flat terrain covered with forest.This value (∼2.6 m) was used for directions where the wind comes solely from the land.For directions where the inlet was located partially on the ground and the water, an interpolation was made between the previously mentionned value and the typical roughness length of water (∼0.0002m).At outlet, a zero-gradient was fixed for all quantities except for pressure for which a value of zero was imposed.  of blackspruce forest was chosen and the drag coefficient was varied (0.15, 0.015, 0.0015) to determine the optimal value.The effect of constant and variable forest height was also examined.For the first scenario, a constant forest height of ten meters was used.For the second case, forest height was determined with the help of maps distributed by the "Ministères des ressources naturelles du Québec".In order not to model the effects of forest twice (through the roughness length and the porous media approach), a roughness length of 0.05 m was applied everywhere forest was present.

Meteodyn parameters
The commercial software Meteodyn WT 4.0 was used as a benchmark for the current experiment.Forest is also treated as a porous media and turbulence is modeled with a one equation closure model based on turbulent length scale.Efforts were made to use approximately the same parameters as the ones used in OpenFOAM when possible.As opposed to OpenFOAM, it is not possible to place the boundaries independently.Based on the dimension from AAV3 to inlet found in the domain sensitivity analysis (135 nodes × 35 m ≈ 5 km) and the distance accounting for smooth introduction of topography (2 km), distance between the center of the domain and lateral boundaries was set to 7 km.To refine a particular location, one must define a result point in the domain, a desired minimum horizontal cell dimension and a horizontal expansion factor.Then, the software starts from this location with the mentionned cell dimension and contructs the mesh up to the boundaries following the given horizontal expansion.A coarse vertical expansion factor was used due to limited computational power available.Mesh parameters used in MeteoDyn are summarized in Table 1.

Results
Figure 1 shows the ratio of horizontal velocity between masts AAV3 and AAV6.The blue forks indicate one standard deviation of the ratio for each direction.
For directions 270 • to 290 • , both sofwares perform very similarly and are at the limit of one standard deviation below the average ratio.For these directions, the wind comes almost entirely from the land.The logarithmic profiles used in these cases, even with a smoothing procedure, are not perfectly suited for forested and complex terrain.Then, for 300 • to 350 • , results considerably improve, especially for MeteoDyn.This seems to be due to the inlet boundary moving more and more in the water, which makes the inflow boundary condition for both models more appropriate.Meteodyn is better reproducing experimental data for 290 • to 320 • .It was found that Meteodyn uses topographic and roughness data up to 2 km outside the computational domain to define its boundary conditions.Consequently, for directions 290 • to 310 • , Meteodyn would then use a wind profile coming mostly from the water, which would explain the better results.This was further investigated in OpenFOAM by pushing back the inlet of case 320 • so that it falls in the water.Results improved as seen in Figure 1.
The variable forest height approach does not seem to significantly improve the results when compared to the constant height approach.The results presented here for a constant forest height (black markers) were obtained using a mean height of 10 m.Note that the distribution of forest on the site of AVV is quite uniform, which may explain the similar performances of both scenarios.However, using accurate forest height data would probably be crucial on a site with highly heterogeneous forest distribution.
It is difficult to comment about the effect of the forest drag coefficient when the wind comes from the land (270 • to 290 • ) for both masts.Results suggest that regardless of the drag coefficient used, the wind deceleration is proportional for both masts, thus yielding similar horizontal velocity ratios.This is no longer valid as the wind comes more and more from the water.Results vary slightly depending on the drag coefficient.Surprisingly, not much difference can be observed between C d = 0.15 and C d = 0.015.For wind directions 330 • to 350 • , where we are most confident about the inlet boundary condition as no assumptions were made regarding z 0 at inlet, C d = 0.0015 is in very good agreement with experimental data and Meteodyn, especially for the last direction.Thus, it can be concluded that the tuning process of the drag coefficient yielded a value of C d = 0.0015.

Modelling LiDAR error
The LiDAR is located roughly 80 m away from the meteorological mast AAV6.Therefore, simulations performed in previous section were used in this section.On top of the error caused by the LiDAR's homogeneity assumption, a difference between LiDAR and cup  measurements is expected caused by terrain separation between both instruments (see Fig. 4).
The method proposed by [3] was used to evaluate LiDAR error.The cartesian velocity components (U x , U y and U z ) obtained with the CFD at locations P1, P2, P3, P4 on Figure 3 are projected onto the lines of sight of the LiDAR resulting in a set of 4 radial velocities.Using equations shown in step 3, the horizontal velocity at the center of the scanned disc can be deduced assuming homogeneous flow.

Results
First, both terrain separation and LiDAR error were analysed separately in Figure 6.Then, Figure 7 compares the combination of both with experimental data.The latter was recorded between 17/02/2011 and 28/03/2011 where the remote sensor and the anemometer were measuring at ∼80 m.
Since there is no major obstacle in the surroundings of the mast and the LiDAR as shown in Figure 5, the small deviations (0−1.5%)from a ratio of one were expected.Both softwares show quite good agrement especially for the interval [300 • −320 • ].Very little deviations from the unity ratio due to LiDAR error are observed.Again, both softwares are within 1% agreement except for [310 • −320 • ] where MeteoDyn is a bit more conservative.Finally, both graphs were combined in order to make a comparison with experimental results.Results obtained with the WAsP Engineering script developed by [3] are also included.The size of the domain is 3 km × 3 km and is centered on LiDAR location.Horizontal resolution was set 20 m, which is the same used by [3].OpenFoam and Meteodyn tend to slightly overestimate the error whereas the opposite is observed with WAsP Engineering.Despite a small difference between modeled and experimental error for direction 320 • , all softwares are well within uncertainty.
Terrain classification was made using topography up to 340 m away from the LiDAR.When looking at Figure 5, there is indeed considerable topographic changes up to this radius.However, when considering only the actual area covered by the disc when the sensor scans at 80 m, which is a circle of diameter ≈85 m, terrain shows less than 5 m variation in elevation and less than 5 • slopes.Therefore, it suggests that the flow seen by the instrument might not be affected by topography as far as ≈340 m away.

Conclusion
First, an atmospheric flow model capable of handling both complex terrain and variable height forest was implemented in the open source CFD package Open-FOAM.Performances of the model were assessed with the help of MeteoDyn WT 4.0 and experimental measurements.The effect of variable forest height as opposed to constant forest height was investigated.Since the distribution of forest is quite homogeneous on the site studied, no significant difference was observed between both approaches.A tuning process was used to determine the most appropriate value of drag coefficient to be used as no experimental data were available for AAV.Among 3 tested values, 0.0015 yielded best results.Results also showed the significant impact of both the location and the nature of the inlet boundary condition.An improvement of the results was noticed as the inlet moves in the water due to the validity of the logarithmic profiles imposed.It was further validated by moving the inlet of case 320 • completely in the water.
Secondly, the flow model developed was used to evalute LiDAR error following the procedure developed by [3].Even though terrain was classified as moderately to highly complex, resulting LiDAR error is very low, which questions both the validaty of the terrain classification and the numerical simulations.The first aspect might deserve more attention since CFD and experimental data agree with each other.Results suggest that the size considered for the terrain classification is of high importance.The LiDAR does not seem to be affected by topographic changes as far as 340 m away from its location, but more by topography very close to the scanned volume.

Fig. 1 .
Fig. 1.Horizontal velocity ratios at 60 m above ground level of AAV3/AAV6 per direction obtained with OpenFoam and Meteodyn accompanied with experimental data.

Fig. 4 .
Fig. 4. Locations of CFD defined points of interest.The prime ( ) means "CFD calculated".The points Mast1' and LiDAR occupy exactly the same location but are calculated differently.Mast2' is the CFD calculated point corresponding to AAV6.