Wind Resource Mapping in Nepal INTERIM MESOSCALE WIND MODELLING REPORT May 2016 This report was prepared by DTU Wind, in association with 3E and WindPower Nepal under contract to The World Bank. It is one of several outputs from the wind Renewable Energy Resource Mapping and Geospatial Planning Nepal [Project ID: P150328]. This activity is funded and supported by the Energy Sector Management Assistance Program (ESMAP), a multi-donor trust fund administered by The World Bank, under a global initiative on Renewable Energy Resource Mapping. Further details on the initiative can be obtained from the ESMAP website. This document is an interim output from the above-mentioned project. Users are strongly advised to exercise caution when utilizing the information and data contained, as this has not been subject to full peer review. The final, validated, peer reviewed output from this project will be the Nepal Wind Atlas, which will be published once the project is completed. Copyright © 2016 International Bank for Reconstruction and Development / THE WORLD BANK Washington DC 20433 Telephone: +1-202-473-1000 Internet: www.worldbank.org This work is a product of the consultants listed, and not of World Bank staff. The findings, interpretations, and conclusions expressed in this work do not necessarily reflect the views of The World Bank, its Board of Executive Directors, or the governments they represent. The World Bank does not guarantee the accuracy of the data included in this work and accept no responsibility for any consequence of their use. The boundaries, colors, denominations, and other information shown on any map in this work do not imply any judgment on the part of The World Bank concerning the legal status of any territory or the endorsement or acceptance of such boundaries. The material in this work is subject to copyright. Because The World Bank encourages dissemination of its knowledge, this work may be reproduced, in whole or in part, for non-commercial purposes as long as full attribution to this work is given. Any queries on rights and licenses, including subsidiary rights, should be addressed to World Bank Publications, The World Bank Group, 1818 H Street NW, Washington, DC 20433, USA; fax: +1-202-522-2625; e-mail: pubrights@worldbank.org. Furthermore, the ESMAP Program Manager would appreciate receiving a copy of the publication that uses this publication for its source sent in care of the address above, or to esmap@worldbank.org. Interim mesoscale wind modelling report for Nepal Jake Badger1 , Andrea N. Hahmann1 , Patrick J. H. Volker1 , Jens Carsten Hansen1 1 Department of Wind Energy, Technical University of Denmark (DTU), Risø Campus, Denmark May 25, 2016 2 METHOD Abstract This document reports on the methods used in Phase 1 of The World Bank wind mapping project for Nepal. The interim mesoscale modelling results were calculated from the output of simulations using the Weather, Research and Forecasting (WRF) model. We document the method used to run the mesoscale simulations and to gen- eralize the WRF model wind climatologies. In addition to the data for Nepal, maps for Bhutan are shown in an appendix. 1 Introduction The conventional method used to produce estimates of wind resource over large areas or regions, such as on a national scale, is to analyze wind measurements made at a number of sites around the region, as in for example the European Wind Atlas (Troen and Petersen, 1989). In order for this method to work well, there needs to be a good spatial coverage of high-quality data. This criterion is sometimes difficult to satisfy and therefore other methods are required that typically give good indications of the geographical distribution of the wind resource, and as such will be very useful for decision making and planning of feasibility studies. Numerical wind atlas methodologies have been devised to solve the issue of insufficient wind measurements. The latest methodology developed at at DTU Wind Energy uses the Weather Research and Forecasting (WRF) model in a dynamical downscaling mode to produce mesoscale analysis. It is this method that is employed in this study and described in this report. The method has recently been documented in Hahmann et al. (2014b) and verified against tall masts in the North and Baltic Sea. This report is structured as follows: Sections 2 and 3 describe the general method and the specific modelling setup of the WRF modelling systems used in the generation of the Nepal phase 1 output. In Section 4 the results are presented, including some examples of local wind climates. Finally, Section 5 presents some conclusions. 2 Method Numerical wind atlas methodologies have been devised to solve the issue of insufficient wind measurements. Two methodologies have been developed and used at DTU Wind Energy. The first methodology is the KAMM/WAsP method developed at Risø National Laboratory. It has been used extensively for a number of national projects. The origins of the method are described in Frank and Landberg (1997) and further details of the downscaling method developed are found in Badger et al. (2014). The KAMM/WAsP methodology has since been upgraded to use a newer and more sophisticated mesoscale model, namely the Weather Research and Forecasting (WRF) model. The wind atlas method used in this study was calculated by carrying out a large number of 10 days mesoscale model simulations using the WRF model to cover a multiyear period. The output from the WRF simulations is analysed in a number of ways. For example, investigation of the dynamic variation of wind speeds as a function of time of day and month of year. Specific meteorological phenomena in the model output relevant to wind energy can be investigated, and an understanding of the important meteorological phenomena is sought. To use the simulation data for wind resource assessment, the data must be post processed. This post 1 3 MODELLING processing includes calculating statistics from a very large dataset and the generalization of the wind climatologies. Wind climate estimates derived from mesoscale modelling and measurements can be compared in a proper way by the use of the generalization of the wind climatologies. Without the generalization step no verification is possible, because the surface description within the model does not agree with reality, and therefore modelled winds will not agree with measured winds, except perhaps in extremely simple terrain or over water far from coasts. The generalization method has been used extensively in a number of wind resource assess- ment studies, particularly within the KAMM/WAsP method. The WRF wind atlas method with generalization and validation was first carried out within the Wind Atlas for South Africa project (WASA, 2014) and described in Hahmann et al. (2014a). For more details on the generalization method see Appendix A. 3 Modelling The WRF Model (Skamarock et al., 2008) is a mesoscale numerical weather prediction system designed to serve both operational forecasting and atmospheric research needs. The simu- lations used to generate the interim wind modelling results utilize the Advanced Research WRF (ARW-WRF) version 3.5.1 model released on 23 September 2013. The WRF modelling system is in the public domain and is freely available for community use. It is designed to be a flexible, state-of-the-art atmospheric simulation system that is portable and efficient on available parallel computing platforms. The WRF model is used worldwide for a variety of applications, from real-time weather forecasting, regional climate modelling, to simulating small-scale thunderstorms. Although designed primarily for weather forecasting applications, ease of use and quality has brought the WRF model to be the model of choice for downscaling in wind energy applications. This model was used in wind-related studies concerning: wind shear in the North Sea (Pe˜ na and Hahmann, 2012) and over Denmark (Draxl et al., 2014), organized convection in the North Sea (Vincent et al., 2012), low-level jets in the central USA (Storm en et al., 2009), wind climate over complex terrain (Horvath et al., 2012), gravity waves (Lars´ et al., 2012), extreme winds (Lars´en et al., 2013), among many others. 3.1 Model setup The simulations for the interim wind modelling were calculated on a grid with horizontal spacing of 45 km × 45 km (outer domain, D1, with 116 × 94 grid points), 15 km × 15 km (first nested domain, D2, with 181 × 121 grid points) and 3 km × 3 km (second nest, D3, with 476 × 206 grid points). Maps of the model domains are displayed in Fig. 1. The surface roughness length for innermost domain, D3, is given in Fig. 2. In the vertical the model was configured with 50 levels with model top at 20 hPa. This is a special model configuration adapted to the occurrence of more deep convective activity in this region. The lowest 10 of these levels are within 1000 m of the surface and the first level is located at approximately 11 m AGL. Table 1 lists the details of the model configuration, including the model parametrizations used in the simulations. The actual namelist used in the simulations is presented in Appendix B. 2 3.1 Model setup 3 MODELLING Figure 1 – WRF model domains configuration and terrain elevation (m). Top left: 45 km × 45 km domain (D1), top right: 15 km x 15 km (D2) and bottom: 3 km × 3 km (D3). The inner lines show the position of D2 and D3 in D1 and D2, respectively. The colour scale indicates the terrain height. 3 3.1 Model setup 3 MODELLING Figure 2 – WRF model domain D3 surface roughness length. The horizontal grid spacing is 3 km × 3 km. The colour bar to the bottom right indicates the values of surface roughness length. 4 3.1 Model setup 3 MODELLING Table 1 – Summary of model and system setup and physical parameterizations used for the WRF simulations. Model setup: WRF (ARW) Version 3.5.1. Mother domain (D1; 116 × 94 grid points) with 45 km grid spacing; 2 nested domains: D2 (181 × 121 grid points) using 15 km and D3 (476× 206 grid points) with 3 km horizontal grid spacing on a Mercator projection (see Fig. 1). 50 vertical levels with model top at 20 hPa; 10 of these levels are placed within 1000 m of the surface; The first 7 levels are located approximately at: 11, 33, 55, 77, 100, 120 and 142 m. MODIS (2001–2010) land-cover classification of the International Geosphere-Biosphere Programme. Simulation setup: Initial, boundary conditions, and fields for grid nudging come from the The ERA Interim (ERA-I) at 0.7◦ × 0.7◦ resolution. Runs are started (cold start) at 00:00 UTC every 10 days and are integrated for 11 days, the first 24 hours of each simulation are disregarded. Sea surface temperature (SST) from Optimum Interpolation Sea Surface Temperature (OISST) at 0.25◦ × 0.25◦ resolution (Reynolds et al., 2010) and are updated daily. Model output: hourly (lowest 18 vertical levels) for D3. Time step in most simulations: approx. 135 seconds. One-way nested domains; 5 grid point nudging zone. Grid nudging on D1 only and above level 10; nudging coefficient 0.0003 s−1 for wind, temper- ature and specific humidity. No nudging in the PBL for temperature and specific humidity. Physical parameterizations: Precipitation: WRF Single-Moment 6-class scheme (option 6), Tiedtke scheme cumulus pa- rameterization (option 6) turned off on D3. Radiation: RRTMG scheme for longwave (option 4) and shortwave (option 4) radiation, which includes the MCICA method of random cloud overlap. PBL and land surface: Mellor-Yamada-Janjic (Eta) TKE scheme (Janjic, 1994) (option 2), Eta Similarity Scheme (option 2) surface-layer scheme, and Noah Land Surface Model (option 2). The surface roughness is kept constant at their winter value. Diffusion: Simple diffusion (option 1); 2D deformation (option 4); 6th order positive definite numerical diffusion (option 2); rates of 0.06, 0.08, and 0.1 for D1, D2, and D3, respectively; vertical damping. Positive definite advection of moisture and scalars. 5 3.2 Data processing 3 MODELLING Most choices in the model setup are fairly standard and used by other modelling groups. The only special setting for wind energy applications is the use of a constant surface roughness length, thus disabling the annual cycle available in the WRF model. This choice is consis- tent with the generalization procedure discussed in section 2 and Appendix A. A few other parameterization settings are updated for equatorial conditions compared to other wind atlas simulations: more vertical levels and raised model top, more sophisticated microphysics and convective scheme and updated radiation parameterizations. Figure 3 – WRF model simulation schematic showing how the simulation period is covered by a succession of overlapping 11 day simulations. The first day of the simulations, which overlaps with the last day of the previous simulation, is for model spin-up and is not used in subsequent analysis. The final simulation covered the 10-year period 2004/08/4– 2013/08/14, and was run in a series of 11-day long overlapping simulations, with the output from the first day of each simulation being discarded, see Fig. 3. This method is based on the assumptions described in Hahmann et al. (2010) and Hahmann et al. (2014b). The simulation used grid nudging that continuously relaxes the model solution towards the gridded reanalysis but this was done only on the outer domain and above the boundary layer (level 10 from the surface) to allow for the mesoscale processes near the surface to develop freely. Because the simulations were re-initialized every 10 days, the runs are independent of each other and can be integrated in parallel reducing the total time needed to complete a multi-year climatology. The grid nudging and 10-days reinitialization keeps the model solution from drifting from the observed large- scale atmospheric patterns, while the relatively long simulations guarantee that the mesoscale flow is fully in equilibrium with the mesoscale characteristic of the terrain. 3.2 Data processing Wind speeds and directions are derived from the WRF model output, which represents 10- minutes or hourly instantaneous values. For evaluating the model wind speed climatology, the zonal and meridional wind components on their original staggered Arakawa-C grid were interpolated to the coordinates of the mass grid. The interpolated wind components were then used to compute the wind speed. For a given height, e.g., 100 m, wind speeds are interpolated 6 4 RESULTS between neighboring model levels using logarithmic interpolation in height. It was found that this interpolation procedure preserves more of the original features in the model wind profile compared to other schemes (e.g., linear or polynomial interpolation of the wind components). The various data processing steps are shown in Fig. 4. Static Data Height 10min Data pro- WRF Raw 60min (WindsT) Time-series Generalization .lib reduction Interpolation cessing Raw 10min (Winds) Figure 4 – Schematic representation of the data processing used to create the wind climate files that compose the WRF-based NWA. For each model grid point inside Nepal in domain D3 time-series for the entire period for the wind speed, wind direction at 5 heights, and 1/L were generated. The generation of the time-series is a rather time consuming process because the WRF output files are stored for every three hours for the whole domain. The generation of time-series requires that for every grid-point in the considered region all files for the whole period have to be accessed. 4 Results In this section the results in the form of the annual mean wind climate are presented based on the 10 years of simulation, covering the period 2004/08/4 to 2013/08/14 inclusive. First the simulated winds are presented. These represent the annual mean wind speed and power density at 100 m a.g.l. directly from the modelling, see Figs. 5 and 6. Therefore, the winds in these maps reflect the orography and surface roughness length as they are represented in the model rather than the real orography and roughness length. Please note for the power density calculation the air density is provided by the mesoscale model simulation. Figs. 7 and 8 show the generalized winds. These represent the annual mean wind speed and power density at 100 m a.g.l. for standardized condition of flat terrain with surface roughness length of 10 cm everywhere. The winds in these maps reflect the variation of the winds due to all influences other than the microscale orography and surface roughness change. Whereas, the generalized power density, which uses a constant air density is the variation of power density due to variation of the wind speed distribution alone. An example of generalized wind climate file data is given in Fig. 9. Figure 10 shows the location of the generalized wind climate files. One generalized wind climate file is created for every WRF model grid point inside Nepal. These files can be used in the WAsP software to calculate the predicted wind climate accounting for highly detailed microscale orography and surface roughness change effects for a particular site of interest. Next the generalized winds are presented. These represent the annual mean wind speed and power density at 100 m a.g.l. for standardized condition of flat terrain with surface roughness length of 10 cm everywhere, as shown in Figs. 7 and 8. Now the winds in these maps reflect 7 4 RESULTS Figure 5 – Mean annual simulated wind speed for Nepal at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The colour scale indicates the wind speed in m s−1 . 8 4 RESULTS Figure 6 – Mean annual simulated wind power density for Nepal at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The colour scale indicates the wind power density in W m−2 . Note: for the power density calculaton the air density is from the mesoscale model simulation. 9 4 RESULTS Figure 7 – Mean annual generalized wind speed for Nepal at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The standard conditions are flat terrain with uniform surface roughness length (10 cm). The colour scale indicates the wind speed in m s−1 . 10 4 RESULTS Figure 8 – Mean annual generalized wind power density for Nepal at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The standard conditions are flat terrain with uniform surface roughness length (10 cm). The colour scale indicates the wind power density in W m−2 . Note: for the power density calculation only the air density is constant at 1.25 kg/m3 . 11 4 RESULTS Figure 9 – Example of the data contained within a generalized wind climate file data. This data can be used in the WAsP software to make predictions of the wind resources at a specific site of interest accounting for the microscale effects due to orography and surface roughness changes. The location is in eastern Nepal. 12 4 RESULTS Figure 10 – Top: The location of the generalized wind climate data for the whole of Nepal shown in Google Earth. Bottom: A detail of generalized wind climate data coverage including how a user of the data can find out about the data filename using Google Earth. 13 4 RESULTS the variation of the winds due to all influences other than the microscale orography and surface roughness change. Please note for the power density calculation only the air density is constant at 1.25 kg/m3 , so that variation of power density is due to variation of the wind speed distribution alone. The figures show only a small part of the information contained in the generalizated wind climate. An example of generalized wind climate file data is given in Fig. 9. Figure 10 shows the location of the generalized wind climate files that are created for every WRF model grid point inside Nepal. These files can be used in the WAsP software to calculate the predicted wind climate accounting for highly detailed microscale orography and surface roughness change effects for a particular site of interest. Results for Bhutan are given in Appendix C. The 100 m above surface simulated and generalized mean wind speed and wind power density are given. The maps presented for Bhutan in Figs. 12, 13, 14, 15 correspond to the Figs. 5, 6, 7, 8 shown in this section for Nepal. 14 5 CONCLUSIONS 5 Conclusions This report has described the Phase 1 mesoscale wind modelling for Nepal. The simulation methodology, the configuration of the WRF model and the generalization method have been reported. The results of the wind modelling are presented in the form of simulated and generalized wind maps, and in the form of generalized wind climate data files. In addition to the result for Nepal, the results for Bhutan are given as this was also covered, intentionally, by the calculation domain. The measurement data is essential to the validation work required in Phase 3. Suggestions for regions for the measurement masts are given in the report titled Candidate Site Identifi- cation Report. Through the measurements programme a better understanding of the wind energy relevant meteorology of the country will be gained, an improved configuration of the modelling system will be developed and tested, and an uncertainty estimate of the final wind atlas can be determined. 15 REFERENCES REFERENCES References Badger, J., H. Frank, A. N. Hahmann, and G. Giebel, 2014: Wind-climate estimation based on mesoscale and microscale modeling: Statistical-dynamical downscaling for wind energy applications. J. Appl. Meteor. Climatol., 53, 1901–1919. na, and G. Giebel, 2014: Evaluating winds and vertical wind Draxl, C., A. N. Hahmann, A. Pe˜ shear from WRF model forecasts using seven PBL schemes. Wind Energy, 17, 39–55. Frank, H. and L. Landberg, 1997: Modelling the wind climate of Ireland. Bound.-Layer Me- teor., 85 (3), 359–378, doi:{10.1023/A:1000552601288}. Hahmann, A., C. Lennard, J. Badger, C. Vincent, M. Kelly, P. Volker, B. Argent, and J. Refslund, 2014a: Mesoscale modeling for the Wind Atlas of South Africa (WASA) project. Tech. rep., http://orbit.dtu.dk/services/downloadRegister/102673293/ Mesoscale_modeling.pdf, DTU Wind Energy. Hahmann, A. N., D. Rostkier-Edelstein, T. T. Warner, F. Vandenberghe, Y. Liu, R. Babarsky, and S. P. Swerdlin, 2010: A Reanalysis System for the Generation of Mesoscale Climatogra- phies. J. Appl. Meteor. Clim., 49 (5), 954–972, doi:{10.1175/2009JAMC2351.1}. Hahmann, A. N., C. L. Vincent, A. Pe˜ na, J. Lange, and C. B. Hasager, 2014b: Wind cli- mate estimation using WRF model output: method and model sensitivities over the sea. International Journal of Climatology, doi:10.1002/joc.4217. Horvath, K., D. Koracin, R. Vellore, J. Jiang, and R. Belu, 2012: Sub-kilometer dynamical downscaling of near-surface winds in complex terrain using WRF and MM5 mesoscale models. J. Geophys. Res., 117, D11 111, doi:DOI10.1029/2012JD017432. Janjic, Z. I., 1994: The step-mountain eta coordinate model: Further developments of the convection, viscous sublayer, and turbulence closure schemes. Monthly Weather Review, 122 (5), 927–945, doi:10.1175/1520-0493(1994)122 0927:TSMECM 2.0.CO;2. Kelly, M. and I. Troen, 2014: Probabilistic stability and “tall” wind profiles: theory and method for use in wind resource assessment. Wind Energy, in press. en, X. G., J. Badger, A. N. Hahmann, and N. G. Mortensen, 2013: The selective dynamical Lars´ downscaling method for extreme-wind atlases. Wind Energy, 16, 1167–1182, doi:10.1002/ we.1544. en, X. G., S. Larsen, and A. N. Hahmann, 2012: Origin of the waves in A case-study Lars´ of mesoscale spectra of wind and temperature, observed and simulated’: Lee waves from the Norwegian mountains. Q. J. R. Meteorolog. Soc., 138 (662, Part A), 274–279, doi: {10.1002/qj.916}. na, A. and A. N. Hahmann, 2012: Atmospheric stability and turbulence fluxes at Horns Rev Pe˜ — An intercomparison of sonic, bulk and WRF model data. Wind Energy, 15, 717–731, doi:DOI:10.1002/we.500. Reynolds, R. W., C. L. Gentemann, and G. K. Corlett, 2010: Evaluation of aatsr and tmi satellite sst data. J. Climate, 23 (1), 152–165, doi:DOI10.1175/2009JCLI3252.1. 16 REFERENCES REFERENCES Skamarock, W. C., et al., 2008: A Description of the Advanced Research WRF Version 3. Tech. Rep. NCAR/TN–475+STR, National Center for Atmospheric Research. Storm, B., J. Dudhia, S. Basu, A. Swift, and I. Giammanco, 2009: Evaluation of the weather research and forecasting model on forecasting low-level jets: implications for wind energy. Wind Energy, 12 (1), 81–90. Troen, I. and E. L. Petersen, 1989: European Wind Atlas. Published for the Commission of the European Communities, Directorate-General for Science, Research, and Development, Brussels, Belgium by Risø National Laboratory. Tuller, S. E. and A. C. Brett, 1984: The characteristics of wind velocity that favor the fitting of a Weibull distribution in wind-speed analysis. J. Appl. Meteor. Clim., 23 (1), 124–134, doi:10.1175/1520-0450(1984)0232.0.CO. Vincent, C. L., A. N. Hahmann, and M. C. Kelly, 2012: Idealized mesoscale model simulations of open cellular convection over the sea. Bound.-Layer Meteor., 142 (1), 103–121, doi: DOI10.1007/s10546-011-9664-7. WASA, 2014: The Wind Atlas for South Africa. [Online], http://wasa.info.org. 17 A DETAILED DESCRIPTION OF GENERALIZATION A Detailed description of generalization A.1 Basic generalization equations The generalization of WRF model winds is an extension of the KAMM/WAsP generalization method described in Badger et al. (2014). In the first step, the time series of wind speed and direction are corrected for orography and roughness change, which are a function of wind direction and height. Given a time series of wind speed, u = u(z, t), and wind direction, φ = φ(z, t), which are functions of height and time, intermediate values, u ˆ, are given ˆ and φ by u ˆ= u (1) (1 + δAo )(1 + δAr ) ˆ = φ − δφo , φ (2) where δAo , δφo and δAr and are generalization factors for orography in wind speed and direction and roughness change, respectively. From the time series of corrected wind speed and direction ”wind classes” are determined. The binning is based on wind direction sectors, wind speed and surface stability according to the Obukhov length as described in section A.2. From the binning, mean values of wind speed, u, and wind direction, φ and typical Obuhov length L˜ , together with the frequency of occurrence, F , of each bin are determined. For simplicity, we will drop the over-bar from the equations that follow, but it is understood that they are applied to the mean values of each bin and not the individual time series values. From the corrected wind speed value we obtain an intermediary friction velocity, u ˆ∗ ˆ κu ˆ∗ = u (3) ln[(z/z ˜ )] ˆ0 ) + ψ (z/L where zˆ0 is the downstream surface roughness length and ψ is a stability correction function that adjust the logarithmic wind profile due to non-neutral stability conditions and κ is the von K´arm´ an constant. The stability correction uses the relationship: −31.58[1 − exp(−0.19z/L)] if x ≥ 0 ψ (z/L) = (4) 2 log[0.5(1 + x)] + log[0.5(1 + x2 )] − 2 tan−1 (x) + 1.5746 if x < 0 where x = (1 − 19z/L). We use this function with a typical value of the Obukhov length from each wind class bin (see table 2). This procedure avoids using the similarity theory on wind profiles that lie outside the bounds of validity of the theory and that sometimes occur in the WRF simulations. In the next step, we use the geostrophic drag law, which is used for neutral conditions to determine nominal geostrophic wind speeds, G ˆ , and wind directions, αG , are calculated, using the intermediate friction velocity and wind direction: 2 ˆ∗ ˆ=u G ln uˆ∗ −A + B2, (5) κ fzˆ0 ˆG = − sin−1 B u φ ˆ∗ , (6) κG ˆ 18 A.2 Sectorization A DETAILED DESCRIPTION OF GENERALIZATION where A = 1.8 and B = 5.4 are two empirical parameters and f is the Coriolis parameter, and ˆG is the angle between the near-surface winds and the geostrophic wind. Near the equator, φ where f can become too large or undefined, it is reset to its value at a latitude of 10◦ . ˆ∗G , for a standard roughness length z0,std , To obtain a new generalized friction velocity, u Equation 5 is reversed by an iterative method, 2 ˆ∗G ˆ=u G ln ˆ∗G u −A + B2, (7) κ f z0,std Finally, the generalized wind speed, uG , is obtained by using the logarithmic wind profile law ˆ∗G u z uG = ln . (8) κ z0,std A.2 Sectorization Table 2 – Stability ranges and typical values used in the generalization procedure. Stability class Obukhov length Typical Obukhov value range (m) ˜ (m) L Very unstable -50 < L < -100 -75 Unstable -100 < L < -200 -150 Near unstable -200 < L < -500 -350 Neutral L < -500; L > 500 10000 Near stable 200 < L < 500 350 Stable 50 < L < 200 125 Very stable 10 < L < 50 30 To apply the generalization procedure to the WRF-model output, winds from the mesoscale model simulations are binned according to wind speed (usually in 2.5 m s−1 bins), wind direction (usually 48 sectors of 7.5◦ width) and seven stability class based on the Obukhov length that is also an output from the WRF simulation. The ranges for the stability classes are listed in Table 2 together with the “typical” length used in the generalization. The procedure is carried out for each model grid point independently. In practice, time series of wind speed and direction at the desired vertical levels and 1/L are extracted from the model output files. The generalization procedure is then carried out on each time series file. A.3 Weibull distribution fit The frequency distribution of the horizontal wind speed can often be reasonably well described by the Weibull distribution function (Tuller and Brett, 1984): kw −1 k kw u u F (u) = exp − , (9) Aw Aw Aw 19 A.3 Weibull distribution fit A DETAILED DESCRIPTION OF GENERALIZATION where F (u) is the frequency of occurrence of the wind speed u. In the Weibull distribution the scale parameter Aw has wind speed units and is proportional to the average wind speed calculated from the entire distribution. The shape parameter k (≥1) describes the skewness of the distribution function. For typical wind speed distributions, the kw -parameter has values in the range of 2 to 3. From the values of Aw and kw , the mean wind speed U ( m s−1 ) and mean power density E (W m−2 ) in the wind can be calculated from: 1 U = Aw Γ 1 + (10) kw 1 3 3 E = ρAw · Γ 1 + (11) 2 kw where ρ is the mean density of the air and Γ is the gamma function. We use the moment fitting method as used in the Wind Atlas Analysis and Application Program (WAsP) for estimating the Weibull parameters. The method is described in detail in Troen and Petersen (1989). Basically this method estimates Aw and kw to fit the power density in the time series instead of the mean wind speed. The Weibull fit is done for the ensemble of wind speeds in each wind direction bin (usually 12 direction sectors) for each standard height (usually 5 heights: 10, 25, 50, 100 and 200 m) and standard roughness lengths (usually 5 roughness: 0.0002 (water), 0.03, 0.1, 0.4, 1.5 m). The 25 Weibull fits for each wind direction sector use the method described above. This sector-wise transformation of Weibull wind statistics—i.e. transforming the Weibull Aw and kw parameters to a number of reference heights over flat land having given reference roughnesses—uses not only the geostrophic drag law, but also a perturbation of the drag law, with the latter part including a climatological stability treatment. The transformation and stability calculation is consistent with that implemented in WAsP and outlined in Troen and Petersen (1989), with further details given in Kelly and Troen (2014). The transformation is accomplished via perturbation of both the mean wind and expected long-term variance of wind speed, such that both Weibull-Aw and kw are affected. When purely neutral conditions (zero stability effects) are presumed for the wind statistics to be transformed, there is still a perturbation introduced, associated with the generalized (reference) conditions in the wind atlas. This perturbation uses the default stability parameter values found in WAsP; it is negated upon subsequent application of the generalized wind from a given reference height and roughness to a site with identical height and surface roughness, using WAsP with its default settings. The climatological stability treatment in the generalization depends on the unperturbed Weibull parameters and effective surface roughness (Troen and Petersen, 1989), as well as the mesoscale output heights and wind atlas reference heights (though the latter disappears upon application of wind atlas data via WAsP). Figure 11 shows the structure of the resulting WAsP ”lib” file. It is structured as Weibull Aw ’s and kw ’s for each sector, height and standard roughness length. The first row contains information about the geographical location of the wind climate represented in the lib-file. The second row lists the number of roughness classes (5), heights (3), and sectors (12), respectively. In the third and fourth row, the actual roughness (m) and heights (m) are listed. Below these header lines, a succession of frequencies of wind direction (1 line), values of Weibull-Aw (1 line) and Weibull-kw (1 line) for each roughness class and height are printed 20 A.3 Weibull distribution fit A DETAILED DESCRIPTION OF GENERALIZATION Figure 11 – Contents of WAsP generalized wind climate file. This climate is for a location close to for each sector (12 sectors per line). This type of file can be used and displayed (Figure 9) in WAsP. 21 B WRF NAMELIST B WRF namelist &time_control auxinput4_interval = 360,360,360, input_from_file = ’’,’’,’’, end_minute = 0,0,0, history_interval = 1440,720,60, start_day = ’’,’’,’’, debug_level = 0, start_second = 0,0,0, io_form_auxinput4 = 2, start_month = ’’,’’,’’, io_form_input = 2, end_year = ’’,’’,’’, end_hour = 0,0,0, end_day = ’’,’’,’’, start_year = ’’,’’,’’, start_hour = 0,0,0, ignore_iofields_warning = ’’, start_minute = 0,0,0, interval_seconds = 21600, restart_interval = 100000, frames_per_outfile = 1, 1, 3, iofields_filename = ’WAfields.txt’,’WAfields.txt’,’WAfields.txt’, io_form_restart = 2, restart = ’’, end_second = 0,0,0, auxinput4_inname = ’wrflowinp_d’, io_form_boundary = 2, end_month = ’’,’’,’’, io_form_history = 2, io_form_auxinput2 = 2, / &domains parent_time_step_ratio = 1,3,5, num_metgrid_levels = 33, time_step_fract_den = 1, e_vert = 41,41,41 grid_id = 1,2,3, i_parent_start = 1, 29, 43, j_parent_start = 1, 25, 41, s_sn = 1,1,1, s_we = 1,1,1, e_we = 116, 181, 476, e_sn = 94, 121, 206, 22 B WRF NAMELIST smooth_option = 2, time_step = ’’, feedback = 0, time_step_fract_num = 0, parent_id = 1,1,2, parent_grid_ratio = 1,3,5, num_metgrid_soil_levels = 4, dx = 45000.0,15000.0,3000.0 dy = 45000.0,15000.0,3000.0 p_top_requested = 5000, max_dom = 4, eta_levels = ’’, / &physics mp_physics = 4, 4, 4 ra_lw_physics = 1, 1, 1 ra_sw_physics = 1, 1, 1 radt = 10, 10, 10 sf_sfclay_physics = 2, 2, 2 sf_surface_physics = 2, 2, 2 bl_pbl_physics = 2, 2, 2 bldt = 0, 0, 0 cu_physics = 1, 1, 0 cudt = 5, 5, 5 fractional_seaice = 1, seaice_threshold = 0., isfflx = 1, ifsnow = 0, icloud = 1, surface_input_source = 1, num_land_cat = 21, num_soil_layers = 4, sst_update = 1, maxiens = 1, maxens = 3, maxens2 = 3, maxens3 = 16, ensdim = 144, / &fdda grid_fdda = 1, 0, 0 gfdda_inname = "wrffdda_d", gfdda_end_h = 300, 0, 0 gfdda_interval_m = 360, 0, 0 fgdt = 0, 0, 0 if_no_pbl_nudging_uv = 0, 0, 0 23 B WRF NAMELIST if_no_pbl_nudging_t = 1, 0, 0 if_no_pbl_nudging_q = 1, 0, 0 if_zfac_uv = 1, 0, 0 k_zfac_uv = 10, 0, 0 if_zfac_t = 1, 0, 0 k_zfac_t = 10, 0, 0 if_zfac_q = 1, 0, 0 k_zfac_q = 10, 0, 0 guv = 0.0003, 0.000075, 0.000075, gt = 0.0003, 0.000075, 0.000075, gq = 0.0003, 0.000075, 0.000075, if_ramping = 0, dtramp_min = 60.0, io_form_gfdda = 2, / &dynamics h_sca_adv_order = 5,5,5, diff_6th_factor = 0.12,0.12,0.12, zdamp = 5000.0,5000.0,5000.0, rk_ord = 3, damp_opt = 0, non_hydrostatic = ’’,’’,’’, km_opt = 4, moist_adv_opt = 1,1,1, v_mom_adv_order = 3,3,3, w_damping = 1, diff_opt = 0, h_mom_adv_order = 5,5,5, time_step_sound = 6,6,6, scalar_adv_opt = 1,1,1, v_sca_adv_order = 3,3,3, khdif = 0,0,0, diff_6th_opt = 2,2,2, kvdif = 0,0,0, dampcoef = 0.01,0.01,0.01, / &bdy_control nested = ’’,’’,’’, relax_zone = 4, specified = ’’,’’,’’, spec_bdy_width = 5, spec_zone = 1, / &grib2 / 24 B WRF NAMELIST &namelist_quilt nio_tasks_per_group = 0, nio_groups = 1, / 25 C RESULTS FOR BHUTAN C Results for Bhutan Additionally to Nepal, the most inner model domain with a 3 km horizontal grid-spacing covered the entire country Bhutan (Fig. 1). In Figs. 12 and 13 the annual mean wind speed and power density at 100 m a.g.l. is plotted for Bhutan. The annual mean generalized wind speed and power density at 100 m a.g.l. for standardized condition of flat terrain with surface roughness length of 10 cm everywhere are shown in Figs. 14 and 15. The figures result from the same simulations as described in Sect. 3. Figure 12 – Mean annual simulated wind speed for Bhutan at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The colour scale indicates the wind speed in m s−1 . 26 C RESULTS FOR BHUTAN Figure 13 – Mean annual simulated wind power density for Bhutan at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The colour scale indicates the wind power density in W m−2 . Note: for the power density calculaton the air density is from the mesoscale model simulation. 27 C RESULTS FOR BHUTAN Figure 14 – Mean annual generalized wind speed for Bhutan at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The standard conditions are flat terrain with uniform surface roughness length (10 cm). The colour scale indicates the wind speed in m s−1 . 28 C RESULTS FOR BHUTAN Figure 15 – Mean annual generalized wind power density for Bhutan at 100 m a.g.l. from WRF simulation at 3 km × 3 km grid spacing for the period 2004/08/4 to 2013/08/14 inclusive. The standard conditions are flat terrain with uniform surface roughness length (10 cm). The colour scale indicates the wind power density in W m−2 . Note: for the power density calculation only the air density is constant at 1.25 kg/m3 . 29