Downscaling wind energy resource from mesoscale to microscale model and data assimilating field measurements

The main objective of this research work is to develop and evaluate several coupling methods between operational Numerical Weather Prediction (NWP) model and Computational Fluid Dynamics (CFD) model and data assimilate the field measurements into the CFD model. To address the problem of high spatial variation of the topography on the domain lateral boundaries between NWP and CFD domain boundaries, 3 methods – translation, extrapolation and Cressman interpolation are used to impose the NWP model data on the CFD domain lateral boundaries. Newtonian relaxation data assimilation technique is used to incorporate the field measurement data into the CFD simulations. These techniques are studied in a complex site located in southern France. Comparison of wind profiles between the CFD simulation, measurements and CFD simulation with data assimilation are discussed. This combination of state-of-the-art techniques in NWP, CFD, and field data assimilation will provide the basis of a more accurate wind resource assessment method.


Introduction
The development of wind energy requires precise and well-established methods for wind resource assessment. During the last two decades linear flow models were widely used for wind resource assessment and micro-sitting. But the linear models inaccuracies in predicting the wind speeds in very complex terrain are well known. This led to the use of Computational Fluid Dynamics (CFD), which is capable of modeling the complex flow in fine details around specific geographic features. Numerical Weather Prediction (NWP) or mesoscale models use mathematical models of the atmosphere and oceans to predict the weather by assimilating observation of the current weather conditions to predict the flow characteristics. NWP models are able to predict the hourly wind regime at resolutions of several kilometres, but cannot resolve the wind speed up and turbulence induced by the topography features finer than 1 km. However, CFD has proven successful in capturing flow details at smaller scales. Therefore, the combination NWP and CFD models can result in a better modeling approach for wind energy applications.
Linear model results showed identical flow fields for two opposite wind directions, while non-linear CFD model results are highly asymmetrical and emphasis the use of CFD in complex terrain [1]. The disadvantage of the linear (WAsP) and CFD models over the mesoscale modeling is highlighted in [2] for a complex terrain in western Norway. The authors conclude that the mesoscale simulation with a complete meteorological model showed large wind variation with the site due to the local topographical features. Therefore, the disadvantage of WAsP and CFD models was that mesoscale wind variations were not taken into account in the modeling.
In recent years many different approaches for coupling of mesoscale (NWP) and microscale (CFD) models were carried out. The authors in [3] studied the impact of coupling a CFD (CFD-Urban) model with the Weather Research and Forecasting (WRF) NWP model for contaminant transport and dispersion. The CFD model prediction was significantly improved when using the wind field produced by downscaling WRF output as initial and boundary condition. The key reason is that the turning of lower boundary layer wind and pressure gradient are well represented in the time varying 3D WRF fields. Similarly, [4] coupled the WRF model to a high resolution CFD (CFD ACE/CFD Urban) using the Model Coupling Environmental Library (MCEL). [5] proposed a methodology for using the field measurements as a link between the NWP (ALADIN) model data and CFD (Code_Saturne) model. This methodology previously developed at EDF R&D, showed significant reduction in error in comparison with WAsP simulation. [6] presents an application of a NWP model combined with a CFD (VENTOS) model for improving wind power forecasting. The CFD simulations for every sector allow to establish linear relations between the reference meteorological mast and the turbines of the wind farm. Using the power curve and NWP outputs, they evaluated the wind speed and direction at the reference mast, and they were able to predict the production up to 72 hours ahead. Methodology proposed by [7] provided a high resolution wind mapping of Granada region in Spain. It combines NWP model (SKIRON) and CFD (CFDWind) model with straightforward link that relies on calibrated "geostrophic wind". The authors [8] placed weight on the need to have a valid link between mesoscale modeling and microscale modeling, that is essential for application and verification of modeling results. Different routes were proposed, ranging from direct application of mesoscale data to the recommended and more sophisticated approach involving corrections at mesoscale and microscale model.
Data assimilation is a NWP concept encompassing method for combining observations into numerical models. In [9], the authors successfully implemented a multiscale weather model, which is designed for simulation of weather processes from synoptic (WRF) scale to microscale (LES) with data assimilation. In [10], the process of assimilating mesoscale model data into CFD simulation was done successfully and the preliminary results indicate that its possible to nudge the wind profile towards the WRF while retaining the mass consistency of the CFD model. It uses the spatially varying WRF data as inflow condition for the CFD (Acusolve) model. And it also assimilates vertical wind profiles of the fine-scale WRF data into the CFD model in order to nudge the CFD solution towards the WRF model data.
In this research work, we use operational mesoscale model data as inflow condition for microscale grid and address the problem of topography variation at the border between mesoscale and microscale grid. Three methods (Translation, Extrapolation and Cressman interpolation) are developed to impose inflow conditions on the microscale grid. Moreover, we use Newtonian relaxation data assimilation technique to assimilate the field measurement inside the microscale grid.

Methodology
In this section we describe the site selected for this application, field measurement campaign, inflow condition from the operational mesoscale model and microscale model used.

Description of the site
The area of interest is located in southern France. The site is considered very complex with strong slopes (locally larger than 30°), valleys and forest (more than 80%). The average altitude of site is 700m. The site is located at the centre of the figure 1 (left), with plain in South-East, valley on the North-West and further to the North-West a plateau at altitude of 1100m is seen. Figure 1 shows the elevation map of the site (40x40 km 2 ), the inner simulation domain of (20x20 km 2 ) and the masts locations. The digital elevation map and land use data have resolution of 25m and 50m respectively and are obtained from IGN (Institut Géographique National, France).

Experimental campaign
A one-year field measurement campaign was led by EDF R&D and EDF Energies Nouvelles between June 2007 and June 2008 in order to provide input and validation data. The measurement set-up included 2 sodars, two 50m masts (M and FP) with cup anemometers and vanes, and a 80m mast (M80) with cup and sonic anemometers, vanes, temperature and humidity sensors. A Remtech PA2 sodar was installed besides the 50m M mast in order to provide the vertical profile of wind and turbulence between 100 m and 600 m. A Scintec SFAS sodar, which provided the vertical profile of wind and turbulence between 30 m and 130 m, was located near the M80 mast. Four sonic anemometers installed on the M80 mast measured the three components of the wind averaged over 10 minutes, and their standard deviation. Table 1 shows the list of the measurement locations, instrument device used, height of measurement and period of measurement. The analysis of the sonic and sodar measurement data reveals that the SFAS sodar underestimated the wind speeds by about 12% in average when compared to the sonic anemometers. However, the ratio between sodar and sonic wind speed is highly dependent on the wind direction and is linked to the topography [11]. In the present case, analyses are available at 0h and 12h and forecasts between 1h and 11h and between 13h and 23h with a one hour step. A new forecast is produced every 12 hours, while the longest range of the forecast is 48 hours (www.cnrm.meteo.fr/aladin/). Figure 2 shows the location of the simulation domain (20X20 km 2 ) and surrounding operational mesoscale (ALADIN) grid. Météo-France provided the hourly wind speed and temperature profile for the year 2007, turbulent kinetic energy and dissipation are deduced from the mesoscale profile using theoretical relations. The operational ALADIN runs provide the initial and boundary condition for the CFD modeling.

Microscale model -CFD
The microscale model used for this research work is an open-source CFD code, Code_Saturne developed by EDF R&D and CEREA [12]. It is a general purpose CFD code, which handles complex geometry and physics. The atmospheric module takes into account the large scale meteorological conditions and thermal stratification of the atmosphere. This module is used for wind energy engineering, urban canopy and pollutant dispersion modelling. It is a finite-volume code, robust for application of Reynolds Averaged Navier-Stokes (RANS) and Large Eddy Simulation (LES) models.
Although RANS is less accurate in comparison with LES, RANS is the most commonly used CFD model for simulation involving turbulent flow in industrial and engineering application because it is computationally less expensive. The turbulence in the simulation domain is parameterized by the standard k-ε turbulence closure model as in the previous work [5]. The atmospheric stability is not taken into account in the simulations.
The inlet turbulent kinetic energy, k and dissipation ε, are generated using the frictional velocity, u* deduced from the ALADIN wind profile. A study was conducted to assess the grid resolution for this domain. The model applies no-slip boundary condition on the ground with roughness value deduced from land use file, velocity inflow condition is prescribed on the north and west face for the meteorological situation considered here, outflow condition on the east and south face and symmetry condition at the top. The simulations are run until a stationary state is reached in the domain.

Coupling and data assimilation
In this section, the methods developed for imposing mesoscale model wind profile on microscale grid and assimilating field measurement into the microscale grid are discussed.

Imposing mesoscale model wind profile on microscale grid
The differences in horizontal (Δxy) and vertical resolution (Δz) between mesoscale grid (Δxy ≈ a few km and Δz ≈ some tens of meters) and CFD grid (Δxy ≈ some tens of m and Δz ≈ some meters) are very large, whereas the topography variations are strong at the border of the simulation domain. The relief seen by the mesoscale and microscale codes are different. In particular, there are some unavailable values of the model variables when the relief altitude is higher in the mesoscale grid than in microscale grid. Therefore, imposing the inlet boundary condition on the CFD domain is a difficult task and different strategies need to be implemented. The methods used to impose mesoscale data on microscale grid are: translation, extrapolation and Cressman interpolation (refer to figure 3).
The translation method was developed and implemented in [5]. In translation, the mesoscale wind profiles are translated along the absolute co-ordinates of the microscale terrain. In this method wind speed varies along the height above ground level which is not reliable far from the ground. In extrapolation, the mesoscale wind profiles are extrapolated between the lowest level of the mesoscale profile and the absolute altitude of the local ground level in the microscale grid. Thus, above a certain level, the inlet wind speed obtained in extrapolation method is function of height above sea level.
Translation Extrapolation Cressman interpolation Finally, Cressman interpolation is a method which is widely used in NWP [13]. Using the Cressman interpolation, the wind speed on a boundary face of the microscale grid is calculated as a linear combination of the values provided by the nearby mesoscale data. Figure 3 shows r L and r Z are longitudinal and vertical radius of influence fixed at the beginning of the computation.

Field measurement data assimilation into the microscale model
Data assimilation was originally developed for NWP, in which meteorological observations enter the models forecast and analysis cycles. Nudging, a Newtonian relaxation data assimilation technique is used to incorporate the measurement data into the CFD simulations. It consists on adding to the prognostic equation for the wind components a term that nudges the solution towards the observations [13,14].  , where r is the radial distance from observation point. This type of function used in NWP and it doesn't vary in time since time integration are very small in CFD simulation. u τ is the relaxation time scale. This type of modeling approach is not common in CFD simulation and Code_Saturne is modified to incorporate this additional forcing term.

Case description
The case chosen for initial analysis is typical of the most dominant wind direction (North-Westerly) and is on 07 th December 2007 and the specific time for the CFD simulation is 16:00 CET. The surface wind speed and temperature are 14 m/s and 10 °C, flows from the plateau (North-west) and is perpendicular to the proposed wind farm site.

Grid resolution
The domain used for this study is located at the centre of the figure 1, has dimension of 20x20 km 2 horizontally and is 5.7 km vertically. The masts were located at the centre of the domain. Topography is built using the digital elevation map with resolution of 25 m for the centre (17x17 km 2 ) that is gradually smoothed to 100 m towards the boundaries (1.5 km). A study was conducted to assess the resolution of the grid for this domain. A total of 4 grids were constructed with different horizontal and vertical resolutions (see table 2). The coarse, medium, fine and very fine grids have vertical resolution of 10 m, 5 m, 2 m and 1 m respectively and contain 0.5, 1.2, 2 and 4.7 million hexahedral elements respectively. The horizontal resolution varies between 10 m and 150 m.     5 show the wake formed behind the hill (in front of the mast location M80) for the coarse and very fine grid. Coarse grid and medium clearly under predict the wake size compared to fine and very fine grids. From figures 6 and 7, we can state that the coarse grids over predict the velocity close to the ground and thus are unable to predict correctly the speed-up induced by the topography properly at mast M80 and FP. Grids with vertical resolution below 2m near the ground are able to predict the flow properly and speed-up at the mast location M80 and FP.

Coupling methods
The CFD simulations are carried out for inflow condition computed using 3 methods -translation, extrapolation and Cressman interpolation. Translation and extrapolation can only use a single ALADIN grid profile to impose the inlet boundary condition for all inlet boundary faces. The ALADIN grid profile can be chosen from either the upstream or middle of the computational domain depending on the prevailing wind direction. The drawback of translation and extrapolation method is the difference in the relative altitude of the chosen ALADIN grid profile and the inlet boundary face altitude that varies locally. Cressman interpolation method uses ALADIN profiles in and around the simulation domain to compute the inflow conditions. A sensitivity analysis for different radius in the Cressman interpolation was carried out. Table 3 shows the different radius of the horizontal and vertical radius of influence used in the simulations. Figures 8  and 9 show ALADIN -U and V velocities (red and green colour respectively) and U and V velocities (black colour) computed using the Cressman interpolation for all inlet boundary faces.   The results of CFD simulations using translation, extrapolation and Cressman interpolation are plotted for 2 mast locations along with field measurements. Figures 10 & 11 show the comparison of translation, extrapolation and Cressman interpolation coupling methods at mast M80 and FP location. At M80 mast location, the three consecutive 10 minutes average of cup and sonic anemometer are in good agreement with each other. The sodar measured reduced average velocity in comparison to the cup and sonic anemometer and shows mixed wind profile along the 100 m height. CFD prediction using Code_Saturne showed that Cressman interpolation and extrapolation coupling methods clearly better predict the mean wind speed, compared to translation coupling method. All the three coupling methods are able to capture the shape of wind speed-up. At M mast location, the wind speed is overpredicted by Code_Saturne compared to both cup anemometer and sodar measurements for all three coupling methods. At FP location, the three consecutive 10 minutes averages of cup anemometer show large variation. The CFD models cannot reproduce this variability, as the results of the simulation correspond to a stationary state for the given boundary condition. Cressman interpolation in another case showed also good prediction compared to extrapolation and translation coupling methods.  Figure 10. Comparison of velocity profiles from different coupling methods at M80. Figure 11. Comparison of velocity profiles from different coupling methods at FP.

Data assimilation
In this section, the implementation and testing of assimilation procedure is carried out for the same case that occurred on 07 th December 2007 at 16:00 UTC. Cressman interpolation is used for imposing the inlet boundary condition for the CFD simulation with data assimilation. The idea is to assimilate one of the mast and use the other two masts are used for validation purpose. The location of M80 mast in between M and FP masts is an ideal choice for field measurement assimilation. For this initial analysis the sonic anemometer measurements are assimilated into the CFD simulation.
The observed velocities obs u in Eq. (4) have to be interpolated in the vicinity of M80 mast location.
The volumetric Cressman interpolation is used to interpolate U and V components of the observed velocity. Sensitivity analysis is carried out to determine the longitudinal radius of influence (r L ) and vertical radius of influence (r Z ) of the volumetric Cressman interpolation. Sonic anemometer measurements are available at 10, 25, 45 and 78 m heights. The volumetric Cressman interpolation is carried out with longitudinal radius (r L ) of 400, 600 and 800 m and vertical radius (r Z ) of 100, 150 and 200 m. Then, the Cressman spreading function, from Eq. (4) is computed. Figures 12 and  13 show the horizontal and vertical cross-section of the Cressman spreading function applied at the M80 mast location respectively. The Cressman spreading function is imposed at 50 m above ground, which corresponds to the centre of the vertical range covered by sonic measurements. The value of the Cressman spreading function, W (x, y, z, t) = 1at 50 m above ground and it drops exponentially outwards to W (x, y, z, t) = 0 . The purpose is to nudge M80 sonic anemometer and monitor the influence of nudging at the M and FP mast location. The radius of influence of volumetric Cressman interpolation and Cressman spreading function are chosen in such way that both masts (M and FP) are inside the assimilation zone. To quantify the value of u τ , sensitivity study was conducted with following values: 35, 50, 60 and 100 s for sonic measurement assimilation.    In the figures 14 and 15, CFD represents the simulation with Cressman interpolation at inlet faces and without assimilation, and CFD+nudging represents the simulation with Cressman interpolation with field measurement assimilation. From the mean velocity profiles (figure 14) at M80 location, CFD+nudging showed good improvement in predicting wind speed using nudging. At M mast location, CFD+nudging improved the prediction compared to the CFD but still shows an overestimation compared to cup anemometer and long-range sodar measurements. At FP mast location CFD+nudging improved the prediction compared to CFD. CFD simulation with assimilation predictions is in good agreement at M80 and FP mast, at M location the results over predict the wind speed by 11% when compared with cup anemometer and by 30% when compared to sodar at 200 m. Two additional case studies for both Northwesterly and Southeasterly wind direction were carried out and it showed that CFD simulation with assimilation at M80 location are able to correct the CFD simulation without assimilation towards the measurement data.

Conclusion
A new methodology has been developed for improving the prediction of annual average wind speed over a complex terrain using coupled mesoscale and microscale model using field measurements assimilation into the micro scale (CFD) model. To better use all the operational mesoscale model data on the microscale CFD grid, a new coupling method, based on Cressman interpolation of all the vertical profiles surrounding the site, was developed and implemented. This coupling method is applicable for any type of terrain but it provides benefits especially in complex terrain. Sensitivity analyses were carried out for: grid size, determination of the optimum radius of influence for the Cressman interpolation at inlet boundary and testing different coupling methods. To improve further the CFD model prediction a data assimilation technique is used. It uses field measurements data from the site to modify the CFD solution towards the measurements (nudging). Nudging technique was implemented into the CFD model and sensitivity analysis were carried out to optimise the radius of influence of the volumetric Cressman interpolation, nudging coefficient and Cressman spreading function. The results with the comparison of annual averaged wind speed was carried out for the year 2007 using measurements, WAsP results, methodology proposed by [5] and current methodology using CFD simulations without assimilation and CFD simulations with assimilation will be published.