Kunal K Dayal

Gilles Bellon

John E Cater

Michael J Kingan

Rajnish N Sharma

High-resolution mesoscale wind-resource assessment of Fiji using the Weather Research and Forecasting (WRF) model

Keywords: Wind energy, resource assessment, Wind power, WRF, Fiji 1

This study presents a high-resolution mesoscale wind-resource assessment of the small island developing state (SIDS) of Fiji using a 10-year simulation of the Weather Research and Forecasting (WRF) model with convection-permitting resolution. Our analysis evaluates the wind speed and Weibull distributions, diurnal and seasonal wind speed patterns, resource maps of annual and seasonal wind speed, power density, model statistical analysis and interannual wind speed variability. The results reveal that the WRF-model simulated wind resource parameters are in good agreement with observations at 24 existing weather stations. At 55 m above ground, the annual mean wind speed and wind power density varies from 1.5 m/s to 8 m/s and 50 W/m 2 to 300 W/m 2 , respectively, for onshore land areas. Higher wind speeds are observed during austral winter than in austral summer. Forty high wind-resource areas are identified in this study, which were previously unknown. This indicates that there is potential for utility-scale wind power generation at selected locations with wind speed and power density greater than 6.4 m/s and 300 W/m 2 (NREL, Wind Power Class 3). An estimated 1000 MW theoretical potential installed capacity is available for utility-scale wind power applications on Viti Levu and Vanua Levu.

Introduction

Climate change and air pollution have propelled a steady growth of renewable electrical energy generation, and in particular wind energy on the global scale. Over the last decade, the global installed capacity of wind energy grew from 159 GW in 2009 to 651 GW in 2019 [START_REF]Renewables 2020 Global Status Report[END_REF]2]. This is an exponential growth with an estimated doubling time of approximately 4 years. Moreover, work is in progress in many countries with new resource assessments, wind farm construction permits and developments [START_REF]Renewables 2020 Global Status Report[END_REF]2].

Wind-resource assessment at different scales is fundamental for energy-resource modelling and site identification for the development of wind farms to support renewable electricity generation. Numerous methods are available for wind-resource assessment from microscale to mesoscale applications with each having its own benefits and disadvantages [START_REF] Brower | Wind Resource Assessment: A Practical Guide to Developing a Wind Project[END_REF]. For microscale wind-resource assessment of a specific site or location within a domain of 10 -20 km, microscale linear models like the Wind Atlas Analysis and Application Program (WAsP) are used. The model uses ground-based measured wind data for wind resource evaluation and resource mapping at a resolution of 50 m × 50 m or less within the domain around a site [START_REF] Brower | Wind Resource Assessment: A Practical Guide to Developing a Wind Project[END_REF].

For regional wind-resource assessment of a country or a group of islands, numerical models can be used to simulate the wind-resource for the region of study [START_REF] Brower | Wind Resource Assessment: A Practical Guide to Developing a Wind Project[END_REF]. The most common numerical weather prediction (NWP) model used for atmospheric studies, which includes resource mapping, is the Weather Research and Forecasting (WRF) model [START_REF] Powers | The Weather Research and Forecasting Model: Overview, System Efforts, and Future Directions[END_REF]. The usefulness of the WRF model is its ability to simulate the atmospheric processes from the ground level to the stratosphere 20 km above. It provides an approximation of the surface atmospheric parameters like surface winds. It also has the advantage of generating mesoscale wind resource maps at different heights for resource identification, theoretical potential evaluation and for creating high-resolution wind data at km scale at locations whereby measurements are not available.

The WRF model has been used for wind energy resource mapping in Palestine [START_REF] De Meij | Wind energy resource mapping of Palestine[END_REF], assessing wind resource in Greece [START_REF] Giannaros | Performance evaluation of the Weather Research and Forecasting (WRF) model for assessing wind resource in Greece[END_REF], modelling wind resource in complex terrain in Portugal [START_REF] Carvalho | Wind resource modelling in complex terrain using different mesoscale-microscale coupling techniques[END_REF], determining wind power characteristics of Oahu, Hawaii [START_REF] Argüeso | Wind power characteristics of Oahu, Hawaii, Renewable Energy[END_REF], assessing offshore wind resource off the Atlantic Iberian coast and in the Gulf of Thailand [START_REF] Salvação | Wind resource assessment offshore the Atlantic Iberian coast with the WRF model[END_REF][START_REF] Chancham | Offshore wind resource assessment and wind power plant optimization in the Gulf of Thailand[END_REF] and for predicting wind energy resources in Tanzania [START_REF] Enock | Application of WRF mesoscale model for prediction of wind energy resources in Tanzania[END_REF]. In these studies, the model assisted in creating new knowledge about mesoscale wind resources and their characteristics, identified and classified the wind-resource sites and provided useful information for power generation purposes. To date, most downscaling studies for atmospheric parameters have been conducted in the midlatitude regions of the Northern Hemisphere and relatively few have examined semi-arid or tropical locations in the Southern Hemisphere.

Fiji is a small island developing state (SIDS), which is located between the latitudes of 12°S -22°S and longitudes of 177°E -178°W in the tropical Southwest Pacific region. There are 332 islands in Fiji with a total land area of 18,333 km 2 , with 110 islands inhabited. Around 87 % of the total land area is occupied by the two largest islands of Viti Levu and Vanua Levu, with areas of approximately 10,400 km 2 and 5,540 km 2 , respectively. The islands of Fiji are mountainous and of volcanic origin, with maximum altitude of 1300 m above mean sea level.

The overall climate in Fiji is categorized as a tropical marine climate with two seasons, which are classified as austral summer (wet season) extending from November to April and austral winter (dry season) from May to October.

In Fiji, electricity generation is from a mixture of renewable (hydro, biomass and wind) and conventional (industrial diesel and heavy fuel oil) power stations [START_REF] Efl | Energy Fiji Limited[END_REF]. In the year 2019, 1060 GWh of electricity was generated with 58 % (611 GWh) from renewable and 42 % (450 GWh) from conventional power plants. The specific contribution from renewable energybased power plants include: 53 % (559 GWh) from hydropower, 4.6 % (49 GWh) from biomass-based independent power producers and 0.25 % (2.7 GWh) from wind [START_REF] Efl | Energy Fiji Limited[END_REF]. Fiji aims to achieve 100 % electricity generation from renewable energy sources by 2030 [START_REF]Fiji's Intended Nationally Determined Contribution[END_REF].

The extent to which wind electricity generation can contribute to this goal is unknown.

Mapping wind resources on the mesoscale level will assist the SIDS of Fiji in terms of wind resource identification and quantification. Several wind-resource assessment studies have been conducted in Fiji using the Wind Atlas Analysis and Application Program (WAsP) [14 -21] and WindPRO [START_REF] Dayal | Offshore wind resource assessment, site suitability and Technology selection for Bligh Waters Fiji using WindPRO[END_REF]. These studies are confined to micro-scale locations covering small areas surrounding a measuring station in the range of 5 -10 km 2 . Moreover, a recent study in Fiji using 5 -6 years of measured automatic weather station wind data in WAsP identified three potential wind farm sites that can be further investigated for utility-scale wind power applications. The WAsP modelling revealed that each of these three sites could incorporate a minimum of 10 MW using Vergnet 275-kW wind turbines to support the electricity grid network [START_REF] Dayal | Wind resource assessment and energy potential of selected locations in Fiji[END_REF]. To the knowledge of the authors, to date, no validated mesoscale wind-resource assessment has been carried out for the SIDS of Fiji.

Previous studies on the evaluation of the WRF model for simulating surface winds and the diurnal cycle of wind speed for the SIDS of Fiji showed that the standard 1.33 km by 1.33 km grid resolution in a 20 km -4 km -1.33 km grid nesting set up may not be an appropriate grid size to evaluate the surface winds for Fiji; as it was suspected that the topography was not properly resolved and a grid sensitivity study was recommended [START_REF] Dayal | Evaluating of the WRF model for simulating the surface winds and the diurnal cycle of wind speed for the small island state of Fiji[END_REF]. Following that recommendation, a grid sensitivity study of the WRF model for simulating the surface winds for the SIDS of Fiji was conducted and it was concluded that the 1 km by 1 km grid resolution of the 15 km -5 km -1 km grid configuration was appropriate to evaluate the wind resources for the SIDS of Fiji [START_REF] Dayal | A grid sensitivity study of the Weather Research and Forecasting (WRF) model for simulating the surface winds for the small island state of Fiji[END_REF].

In this paper, we aim to answer two research questions: I.

Can the Weather Research and Forecasting (WRF) model be used to simulate the wind flow around the Fijian Islands in order to get sufficiently accurate information on the wind resources/potential? II. Do the WRF simulations reveal sites that are promising for wind power that have not yet been identified?

This paper is an enhancement of the WRF mesoscale wind resource assessment approach applied to small tropical islands from an initial study [START_REF] Dayal | Evaluating of the WRF model for simulating the surface winds and the diurnal cycle of wind speed for the small island state of Fiji[END_REF] followed by a grid sensitivity study [START_REF] Dayal | A grid sensitivity study of the Weather Research and Forecasting (WRF) model for simulating the surface winds for the small island state of Fiji[END_REF] to finally, the decade-long mesoscale wind resource assessment of Fiji. This study is also an application and evaluation of the mesoscale wind resource modelling using the WRF model.

Therefore, we present the high-resolution mesoscale wind-resource assessment of the SIDS of Fiji, which includes the analysis of the wind speed and the Weibull distributions, diurnal and annual wind speed patterns, resource maps for wind speed and power density. The study also presents the 10-year annual variability of wind speed and the model validation against ground-based observations via a model statistical analysis. The purpose of this study is to test a methodology to assess mesoscale wind-resource on tropical islands and to provide an assessment of the wind-resource in Fiji. Section 2 of this work describes the methodology, the results are presented in section 3, then discussed in section 4 and the conclusions are presented in section 5. A summary of the characteristics of the AWSs from which the wind observations are retrieved is shown in Table 1. The elevation of the AWSs is given in metres above sea level (m a.s.l.). The model terrain elevation of each AWS is the grid cell average taken from the grid cell in which the AWS is located.

Methodology

Study area and wind measuring stations

Model and simulation set-up

The study uses the Weather Research and Forecasting (WRF) model version 3.9.1.1 of the Advanced Research (ARW) solver, a widely used community mesoscale model developed by the National Centre for Atmospheric Research (NCAR) [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF][START_REF] Wang | WRF-ARW User Guide V3[END_REF]. The dynamical downscaling method was used for running simulations, whereby coarse resolution output from an analysis using a General Circulation Model (GCM) are used as initial and boundary conditions to drive a regional numerical model to simulate atmospheric parameters via nesting, considering the local conditions [28 -30].

The initial and boundary conditions for the SIDS of Fiji were obtained from 6-hourly NCEP-FNL (Final) Operational Global Analysis data at grid resolution [START_REF]NCEP FNL Operational Model Global Tropospheric Analyses[END_REF]. The static fields for the topography were obtained from the USGS GMTED2010 dataset. The landwater masks, land use/land cover classification, albedo etc. were obtained and interpolated from the 21-class MODIS and MODIS FPAR. These were taken from the National Center for Atmospheric Research (NCAR) database, at a resolution of 30-arc-seconds [32]. Time varying sea surface temperatures (SSTs) were supplied to the model from NCEP-NOAA at a grid resolution [START_REF]SST Accessed[END_REF].

The input data from NCEP-FNL were downscaled by the WRF model using three nested domains as illustrated in Figure 2. The outermost domain is referred to as d01, the second d02 and the innermost as d03 (which is also the domain shown in Figure 1). The two-way nesting approach was applied for the parent-child grid interaction, where the flow of information goes from the coarser domain to the finer domain, with feedback from the inner domains. All domains are centred at a latitude of 17.73°S and longitude of 177.94°E, which corresponds to the centre of Viti Levu using the Mercator projection as recommended for low latitudes [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF][START_REF] Wang | WRF-ARW User Guide V3[END_REF].

The spatial coverage of each of the three domains are:

, and at a horizontal resolution of , and , respectively. The vertical spacing of the three domains consist of 30 terrain-following vertical coordinates, including eight levels below 1 km, as meteorological parameters of interest are in the lower part of the atmosphere, close to the surface. A timestep of 10 seconds is used in d03 following a 3:1 parent-child grid timestep ratio to ensure numerical stability as recommended in refs. [26 and 27]. The choice of nested grids follows the recommendation of the earlier grid sensitivity study [START_REF] Dayal | Evaluating of the WRF model for simulating the surface winds and the diurnal cycle of wind speed for the small island state of Fiji[END_REF], in keeping with the standard practice of using a 3:1 or 5:1 parent-grid ratio [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF][START_REF] Wang | WRF-ARW User Guide V3[END_REF].

The well-tested tropical suite of the WRF physics parameterization scheme was used for the model physics options [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF][START_REF] Wang | WRF-ARW User Guide V3[END_REF]. The parameterization of deep convection was only used on the larger domains (d01 and d02), as the vertical fluxes due to convection updrafts and downdrafts and the compensating motion outside the clouds cannot be resolved explicitly at grid scales of 5 km and above [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF][START_REF] Wang | WRF-ARW User Guide V3[END_REF]. The remaining physics options are the same for all three domains as shown in Table 2. Simulations for the 10 years were performed in monthly segments on the New Zealand eScience Infrastructure (NeSI) High Performance Computer -Mahuika. The model was initiated from 3 days prior to the start of each month and these 3 days are the model spin-up time and are excluded from the analysis. The model is set to restart every 2 days to avoid model drift as recommended for longer simulations [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF][START_REF] Wang | WRF-ARW User Guide V3[END_REF]. Even though a period of one year has been considered adequate to capture the diurnal and seasonal variations for wind resource assessment [START_REF] Al-Yahyai | Review of the use of numerical weather prediction (NWP) models for wind energy assessment[END_REF], in our case, we use 10 years of simulations to provide a better representation of the diurnal and seasonal variations. The simulation results discussed in the paper are 1-hourly 10-metre-elevation wind speeds retrieved from the simulation on the innermost domain (d03) considering the grid cell containing the AWS as representative of the station.

Model Validation

The validation of a numerical weather prediction (NWP) model can be done according to various methodologies that complement each other [START_REF] Sr | Mesoscale Meteorological Modelling[END_REF]. Five statistical parameters are used for the model evaluation. The Root Mean Square Error (RMSE) is defined as:

in which is the simulated variable at the data point, is the corresponding observed variable and N is the total number of simulation-observation pairs of values.

The mean Bias was used for the evaluation of data tendency. A positive (negative) bias means the simulations overestimate (underestimate) the measured values:

The Standard Deviation Error (STDE) was used to evaluate the dispersion of the error between the observed and the simulated data:

The RMSE, Bias and STDE was used to evaluate the hourly-simulated wind speed against the observed wind speed.

Most emphasis is given to STDE, and this assumption comes from the fact that, even if a simulation has a high RMSE or Bias, if the STDE is low it means that the error is somewhat constant and can be seen as an offset and the simulation physics is correct. If the simulation has a high STDE, the error is more scattered and the simulations have low physical meaning [START_REF] Hong | Dynamical downscaling: Fundamental issues from a NWP point of view and recommendations[END_REF][START_REF] Carvalho | A sensitivity study of the WRF model in wind simulation for an area of high wind energy[END_REF].

The surface wind speed distribution is well approximated by a Weibull distribution [START_REF] Sathyajith | Wind Energy: Fundamentals, Resource Analysis and Economics[END_REF][START_REF] Manwell | Wind Energy Explained: Theory, Design and Application[END_REF] and this distribution is therefore used to evaluate the model-simulated wind distributions against observed wind distributions. The Weibull distribution is a two-parameter function, mathematically represented by a probability density function (PDF), while the cumulative distribution function is represented by . These functions are defined as: and where is the reference wind speed at the height of measurement, the non-dimensional shape parameter, and the scale parameter whose dimension coincides with that of .

Note that . The scale parameter, , indicates on average how windy the site is, whilst the shape parameter, , indicates how 'peaked' the distribution is. The model and the measured Weibull and parameters are used for comparison.

The Pearson correlation coefficient was used to evaluate the diurnal cycle of wind speed and the Weibull probability density functions. The Pearson correlation coefficient is a statistical measure that quantifies the strength of the relationship between the relative variation of two variables; in our case, the simulated and the measured yearly averaged diurnal values and the Weibull parameter values.

Resource maps

The 10-year model simulation results from domain 3 (d03) were used to generate the annual mean wind speed, seasonal mean wind speed and wind power density resource maps at 55 m elevation, which is equivalent to the typical wind turbine hub height.

This was done to assist in identifying potential wind resource sites for utility-scale wind power applications.

Interannual wind speed variability

Interannual wind speed anomalies were calculated by subtracting the measured long-term mean wind speed from the model annual mean wind speed for individual stations for 10years. These were plotted to study the variation of the wind speeds from year to year.

Results

This section presents the results in two parts: the model validation and the resource maps.

Model Validation

The 10-metre-elevation results for the parameters used to validate the model against groundbased measurements, include: the model statistical analysis, wind speed distribution and Weibull PDFs, and the diurnal and annual wind speed patterns. The model results are compared with locations where the measurements are available from AWSs.

Model Statistical Analysis

Tables 345present the model statistical analysis for 3 years for the leeward side, windward side and the outer islands AWSs, respectively. These also present the Since the mean STDE for all AWSs is low, which is below the border threshold of ± 2.50 m/s, this indicates that the model adequately captures the physical processes [35, 36 and 37].

The 10-metre-elevation model errors reported are assumed to be the equivalent at higher elevations of interest, for instance at the hub-height of 55 m.

Wind Speed Distribution and Weibull PDFs

The wind speed distribution provides site-specific information in terms of the Weibull and parameters. Tables 6 -8 present the wind speed distribution parameter analysis for 3 years for the leeward side, windward side and the outer islands, respectively. Higher Weibull (> 4 m/s) and (> 2) parameters are observed and modelled at Momi, Rakiraki, Yaqara, Udu, Koro Island, Levuka, Vunisea, Viwa and Yasawa. These sites are considered windy and have little variation from the mean wind speed, which indicates the presence of consistent, regular wind resources. The comparison of the Weibull PDFs provides an opportunity to visualise and compare the measured and the modelled wind speed and its frequency of occurrence. Figures 3456present the comparison of the 3-year modelled and the measured Weibull PDFs together with the 10years modelled long-term average (WRF LA) for the leeward side, windward side, outer islands and for all the AWSs, respectively. These figures also list the Weibull and parameters for each PDF. The comparison between the modelled and observed Weibull PDFs show that the model agrees moderately well with the observations (with an average correlation of 0.88 for 24

AWSs). Nonetheless the model tends to record lower frequency of lower wind speeds of 1 -2 m/s on the leeward and the windward sides of the larger islands, and 1 -3 m/s for the outer islands compared with the observations. The opposite is seen for wind speeds from 3 -6 m/s for the leeward side, 3 -8 m/s for the windward side, and 5 -10 m/s for the outer islands.

Similar model behaviour is also observed for all the 24 AWSs shown in Figure 6. This indicates that the model tends to overestimate lower wind speeds and underestimate medium to higher wind speeds. The measured Weibull and parameters are slightly lower in magnitude than modelled. The mean Weibull PDF correlation across all stations is 0.88, with 0.89 for the leeward side, 0.86 for the windward side and 0.87 for the outer islands indicating that the results are in good agreement. The 3-year and the 10-year modelled Weibull PDFs have little variation in terms of the Weibull and parameters and the Weibull PDFs correlation is almost 1, which shows that the 3-year period is representative of the 10-year period used in the simulations.

Diurnal Wind Speed Patterns

The diurnal wind speed pattern is the variation of wind speed over a 24-hour period. For a tropical island, it is a representation of the land-sea breeze circulation and provides useful information about the strength of the wind speed during different times of the day. the AWSs, the leeward side AWSs, the windward side AWSs and the outer islands' AWSs, respectively. Note that the diurnal cycles of individual stations are not included in this paper, but rather an average across all stations, stations on the leeward side, windward side and the outer islands is used. The model simulated diurnal cycle for 3 years and 10 years is generally in good agreement with the observations. The overall diurnal cycle for all the 24 AWSs is similar to the diurnal cycles on the leeward and the windward sides. Comparison of the 3 year and the 10-year modelled diurnal cycles shows little variation, with a correlation of almost 1 (Figure 7). The model is able to capture the diurnal cycle patterns for both sides of the larger islands, although it systematically overestimates the wind speed. For the leeward and the windward sides of the larger islands, a strong diurnal cycle is observed and modelled, which can be seen in Figures 8 and9. For the outer islands, which are relatively smaller, a weaker diurnal cycle is observed (Figure 10), which is expected for small tropical islands as ventilation prevents the genesis of the land-sea breeze circulation [START_REF] Ulrich | Superenhancement of Precipitation at the Center of Tropical Islands[END_REF]. This small observed diurnal cycle is not correctly simulated by the WRF model, which only simulates a small early-morning increase in wind speed. This is probably caused by the lack of representation of small-scale topography in the model, which can reduce the efficiency of the ventilation process. The model performance is better on the leeward-side stations with a smaller estimate of the mean diurnal wind speed in comparison with the windward-side and the outer-islands stations. The model-simulated seasonal cycle of wind speed is in good agreement with the observations for the study period. The model successfully captures the seasonal cycle of wind speed for either sides of the larger islands and for the outer islands, although modelled wind speed is overestimated in comparison with the measurements. Again, this overestimation is smaller for the leeward side of the larger islands than for the windward side and the outer islands. A similar model behaviour is observed for the seasonal cycles for all the AWSs. As per the station plots, slightly higher wind speeds are observed during austral winter (May -October), during which the southeast trade winds are more persistent than during austral summer (November to April), during which the wind climatology is mostly dominated by synoptic systems. Little variation can be observed between the 3-year and the 10-year modelled seasonal cycle.

Seasonal Cycle of Wind Speed

Resource Maps

Long-term-mean Wind Speed Map at 55 m

The long-term-mean wind speed map represents the spatial distribution of the average wind speed over the simulation. These maps at higher heights such as at the hub-height of a wind turbine are useful for identifying potential wind resource sites that can be used for utilityscale wind power applications and for siting wind turbines for wind farms. m/s to 8 m/s. The dominant wind direction is south-easterly and the blocking effect of the islands can be seen in their wakes as well as in the upwind flow; this blocking effect means that the flow is diverted around the islands, accelerating the wind on the sides of the islands, particularly in the channel between the two main islands. Following the NREL Wind Power Classification [41] using the 50-meter-elevation reference, numerous potential wind resource locations on shore can be identified with wind speeds greater than 6.4 m/s, which corresponds to Wind Power Class 3.

Seasonal Wind Speed Maps at 55 m

The seasonal wind speed map represents the seasonal spatial distribution of the average hourly wind speed over a season in a year.

Long-term-mean Wind Power Density Map at 55 m

The long-term-mean wind power density at the hub-height of a wind turbine is also useful in identifying potential wind resource sites for utility-scale wind power applications. Figure 18 presents the 55-metre-elevation long-term-mean wind power density map calculated using 10 years of data. The plotted values are the wind power density area-averaged over a grid cell using a constant air density of 1.225 kg/m 3 . Figure 18. Long-term-mean 55-metre-elevation wind power density map of Fiji.

The 10-year average wind power density varies from 50 W/m 2 to 300 W/m 2 at various locations depending on its topography and geographic location for the onshore land area of the major and the smaller islands. The offshore wind power density varies from 100 W/m 2 to 300 W/m 2 from near-shore areas to open sea areas, respectively.

Interannual Wind Speed Variability

The interannual wind speed variability refers to the magnitude of fluctuation in the annual mean wind speed on a year-to-year basis at a site. It is an important consideration for the assessment of the long-term wind resource and for the confirmation of the consistency of the wind speed over a time period of multiple years. To evaluate the interannual wind speed variability, the variance of the wind speed anomaly for the model is plotted and the interannual wind speed anomalies for the measured and the modelled wind speed are plotted separately.

Figure 19 presents the 10-year mean variance of the 10-metre wind speed anomaly for the model for the 24 AWSs. The variance of the model anomaly varies with a minimum of 0.03 m/s at Nadi to a maximum of 2.02 m/s at Viwa with an average of 0.68 m/s for the 24 AWSs.

Eighteen out of twenty-four stations have the variance of the model anomaly less than 1 m/s. The trends in the interannual wind speed anomalies for the model and the measurements on the leeward side, windward side and the outer islands show that there is little variation (< 0.8 m/s) in the measured and the modelled wind speed over the 3 -10 year period. Since there is little variability in the wind speed, energy production would not change significantly from year to year.

Systematic Bias in Wind Power Density

An evaluation of the systematic bias in wind power density at the 24 AWSs for 3 years is presented in Table 9. This shows the overestimation of the model wind speed. Wind power density is calculated using the measured and the modelled wind speed measurements and the difference is calculated to determine the systematic bias in the 10-metre wind power density at the 24 AWSs. The wind power density (WPD) is calculated using the standard equation ( , where is the air density and is the wind speed) with site-specific air density ranging from 1.084 kg/m 3 to 1.168 kg/m 3 . The site-specific air density has been estimated using the WAsP Air Density Calculator which is a function of altitude (elevation in metres above sea level) and mean air temperature at the same altitude. A lapse rate of 6.5 K/km and a sea level pressure of 1013.25 hPa are assumed and the estimate is based on a standard atmosphere. WAsP software uses the barometric reference data from CFSR reanalysis database [START_REF] Mortensen | Wind Atlas Analysis and Application Program: WAsP 12 Help Facility[END_REF].

The systematic bias computed for wind power density is assumed to be the error in the energy production as it is a function of power density and the other parameters considered in the energy calculation equation are constant. The systematic bias for wind power density is assumed to be constant with height, as the non-linear relationship of wind power density and wind speed has been considered in the systematic bias at 10-metre elevation. A lower systematic bias can be expected at higher heights for the wind power density in a similar manner to the wind speed bias.

The mean systematic bias in the wind power density at the 24 AWSs is 9.68 W/m 2 , which is an overestimation of 33.8 %. The systematic bias is less than 15 W/m 2 at 13 AWSs, 15 -30 W/m 2 at 6 AWSs and greater than 30 W/m 2 at 5 AWS. Even though it is assumed that the wind speed bias is constant with height, there is a chance that this bias will be smaller at higher elevations as there will be little surface friction at the hub-height of 55-metres in comparison with that at an elevation of 10-metres.

Discussion

The modelling exercise using 10 years of data presents an opportunity to quantify and validate the wind-resource for the SIDS of Fiji for the first time at a mesoscale level. The diameter (160 m) spacing across and downstream using a 6 x 6 array size [43, 44 and 45], the total available potential wind resource area of ( ) corresponds to an estimated 1000 MW theoretical potential installed capacity. To confirm this, micrositing with a higher resolution model is required to make a decision considering the terrain complexities and the errors associated with the mesoscale modelling.

For the smaller outer islands there is also potential for wind power development for rural electrification as the wind resource is available. However, such a development would have to be complimentary to other power systems such as diesel power plants which currently provide electricity in these locations. Outer islands like Gau, Ovalau, Taveuni and Kadavu can accommodate a few MW scale wind farms as there are isolated mini-grid electricity networks, which are mostly supported by diesel power plants and a micro-scale hydropower plant in the case of Taveuni.

Other work has identified outer island sites of Gau [START_REF] Singh | A Pre-feasibility Study of Wind Resources in Vadravadra, Gau Island, Fiji[END_REF], Kadavu [START_REF] Sharma | Wind energy resource assessment for the Fiji Islands: Kadavu Island and Suva Peninsula[END_REF] and Taveuni [START_REF] Kumar | Wind characteristics and energy potentials at Wainiyaku Taveuni, Fiji[END_REF] as potential sites for rural electrification, which the current mesoscale work also identifies.

Comparing the results of this study with a recent WAsP microscale wind resource study undertaken by the authors [START_REF] Dayal | Wind resource assessment and energy potential of selected locations in Fiji[END_REF], the WRF model is able to successfully identify the three reported potential utility-scale wind farm sites of Rakiraki, Nabouwalu and Udu. Looking at the seasonal cycle of wind speed, the measured and the modelled wind speed seasonal variability is small (< 1 m/s) over a period of 3 -10 years at all sites where measurements are The WRF model has a tendency to overestimate lower wind speeds and underestimate higher wind speeds [START_REF] Al-Yahyai | Review of the use of numerical weather prediction (NWP) models for wind energy assessment[END_REF]. The overestimation of wind speed can be observed at the following station locations where the measured wind speed is low: Keiyasi, Labasa, Seaqaqa, Wainikoro, Nausori, Saqani and RKS on the larger islands (% Bias > 49 %). The underestimation of wind speed can be observed at Rakiraki, Udu and Yasawa whereby the measured wind speeds are medium to high. Similar model behaviour is reported for studies done for Greece [START_REF] Giannaros | Performance evaluation of the Weather Research and Forecasting (WRF) model for assessing wind resource in Greece[END_REF] and Hawaii [START_REF] Argüeso | Wind power characteristics of Oahu, Hawaii, Renewable Energy[END_REF].

The results generated by numerical weather prediction models are only approximate due to simplifications of the topography and the physical processes which they use in order to make the problem tractable [START_REF] Carvalho | Wind resource modelling in complex terrain using different mesoscale-microscale coupling techniques[END_REF]. There is a substantial difference in terms of the wind speed measurements available at 10-metre-elevation from AWSs, which are essentially measurements made at a single point, whereas the wind speeds calculated using the model are the grid cell average which corresponds to an area of . The grid cells are relatively large with an assumed topography and roughness that can be very different from the real conditions [START_REF] Larsén | The selective dynamical downscaling method for extreme-wind atlases[END_REF]. Also, the measurements are made at an elevation of 10-metres, where structures in the close vicinity might influence the measurements. Such effects are not included in the WRF mesoscale model.

There is also a possibility that the stations are out of calibration, as for these sites, the model is able to capture the diurnal cycle pattern, which is similar to the measured and the long-term average diurnal cycles but there is overestimation of wind speed.

The terrain representation in the model can also contribute to errors. In Table 1 there are differences in the model and the real elevation at almost all station locations of between -62.4 m (-8 %) and 54.4 m (170.5 %). For most stations where the topography is less resolved in the model than the real topography, there is an overestimation of wind speed. These include the 6 AWSs mentioned earlier where there is a positive larger difference between the model and the measurements indicating overestimation of wind speed.

Significant sources of error in any numerical weather prediction model (such as the WRF model presented in this paper) include: the initial and the boundary conditions (in our case we are using coarse resolution output results from GCMs, which is almost the same for the entire domain as the resolution is around ); resolution of the domain (in our case it is an average of ); terrain and vegetation characteristics (in our case we used average values for each grid cell); nudging and assimilation techniques; and the complexity of the terrain effects and their representation in the model (in our case we have seen differences) [START_REF] Awan | Parameterization-Induced Error Characteristics of MM5 and WRF Operated in Climate Mode over the Alpine Region: An Ensemble-Based Analysis[END_REF]. The choice of the physical schemes can also contribute to errors in the model, as these schemes are based on many assumptions, and these assumptions may fail, or give an inadequate response to certain synoptic forcing, limiting their application.

In terms of the applicability and value, the finest reanalysis data available globally is ERA5

that is available at

Conclusions

This study presents a high-resolution mesoscale wind-resource assessment of the small island developing state (SIDS) of Fiji using a 10-year simulation of the Weather Research and Forecasting (WRF) model with convection-permitting resolution. Our analysis evaluates the wind speed and Weibull distributions, diurnal and annual wind speed patterns, model statistical analysis, resource maps of annual and seasonal wind speed, power density and interannual wind speed variability and anomalies.

The results revealed that the WRF model simulated wind-resource parameters are in good agreement with the observations. The model can be reliably used to simulate the wind flow around the Fijian Islands in order to get reliable information on the wind resource. At 55 m above ground, the annual mean wind speed and wind power density varies from 1.5 m/s to 8 m/s and 50 W/m 2 to 300 W/m 2 , respectively, for onshore land areas. Higher wind speeds are observed during austral winter than in austral summer. Forty high wind-resource areas are identified in this study, which were previously unknown. This indicates that there is potential for utility-scale wind power generation at selected locations with wind speed and power density greater than 6.4 m/s and 300 W/m 2 (NREL Wind Power Class 3 at an elevation of 50 m).

As the first mesoscale wind resource modelling exercise for the SIDS of Fiji for the recent decade (2009 -2018), our study supports the idea that the SIDS in the Pacific can utilize their wind-resource to ease its dependence on diesel-based power generation. It also provides a methodology for mesoscale wind-resource assessment for the SIDS in the Southwest Pacific together with the identification and verification of potential wind resource sites.

In future work, it is recommended that the WRF model results are coupled with the linear wind flow model WAsP. Such an approach could improve the accuracy of the wind speed prediction by taking into account high-resolution micro-scale features of the topography and vegetation characteristics. This could also assist in determining the accuracy of the wind speed prediction as well as evaluating the wind energy potential for the identified locations from the high-resolution mesoscale wind resource mapping exercise.

Acknowledgement

This study is part of a PhD research project, for which the lead author is the recipient of scholarship funding from the New Zealand Ministry of Foreign Affairs and Trade (NZMFAT). The authors gratefully acknowledge the Fiji Meteorological Services for providing the ground-based wind data measurements from AWSs in Fiji and the New Zealand eScience Infrastructure (NeSI) for the use of the high-performance computing facilities in this research.

Figure 1

 1 Figure1shows the model topographical height contours at a resolution of 30-arc-seconds for the study area, which covers the larger islands and the smaller islands in the surrounding waters and is marked with the locations of 24 automatic weather stations (AWSs) of Fiji Meteorological Services. Time series and uninterrupted measured wind data is available for 3 years (2016 -2018) from the 24 AWSs at an elevation of 10 m above ground level. This data is used for the model validation exercise for the surface winds.

Figure 1 .

 1 Figure 1. Topographical map of the study area (Fiji) with the locations of the AWSs.

Figure 2 .

 2 Figure 2. Simulation domains.

Figure 3 .Figure 4 .

 34 Figure 3. Comparison of the Weibull PDFs for the Leeward Side Stations.

Figure 5 .

 5 Figure 5. Comparison of the Weibull PDFs for the Outer Islands Stations.

Figures 7 -

 7 10 present a comparison of the modelled (WRF) and the measured (M) diurnal wind speed patterns for 3 years together with the 10-year modelled long-term (WRF LA) average for all

Figure 7 .

 7 Figure 7. Diurnal wind speed pattern for all 24 AWSs.

Figure 8 .

 8 Figure 8. Diurnal wind speed pattern for the Leeward Side AWSs.

Figure 9 .Figure 10 .

 910 Figure 9. Diurnal wind speed pattern for the Windward Side AWSs.

Figures 11 -Figure 11 .

 1111  present the comparison of the modelled (WRF) and the measured (M) seasonal cycle of wind speed for the 3-year period together with the 10-year modelled longterm (WRF LA) average for the leeward side, windward side, outer islands and for all the AWSs, respectively.

Figure 12 .

 12 Figure 12. Seasonal cycle of wind speed for the Windward Side AWSs.

Figure 13 .

 13 Figure 13. Seasonal cycle of wind speed for the Outer Islands AWSs.

Figure 14 .

 14 Figure 14. Seasonal cycle of wind speed for all the 24 AWSs.

  Figure 15 presents the 55-metre-elevation long-term-mean wind speed map for the SIDS of Fiji. Ten years of data was used to calculate the long-term mean. The 55-metre-elevation has been chosen because a typical wind turbine hub height, which is suitable for power generation in the SIDS of Fiji, has a hub height of 55 m. Note that wind turbines with a hub-height of 55 m are employed at the existing Butoni wind farm in Fiji[START_REF] Efl | Energy Fiji Limited[END_REF].

Figure 15 .

 15 Figure 15. Long-term-mean 55-metre-elevation wind speed map of Fiji.

  Figures 16 and 17 present the 55-m-elevation mean seasonal wind speed map of Fiji for the two seasons of austral summer and austral winter calculated using 10 years of data.

Figure 16 .

 16 Figure 16. Seasonal 55-metre-elevation wind speed map of Fiji for austral summer.

Figure 17 .

 17 Figure 17. Seasonal 55-metre-elevation wind speed map for austral winter.

Figure 19 .

 19 Figure 19. Variance of wind speed anomaly for 10 years for 24 AWSs.

  leeward side, windward side and the outer islands at 10 m elevation, respectively.

Figure 20 .

 20 Figure 20. Measured interannual wind speed anomalies for the leeward side AWSs.

Figure 21 .Figure 22 .

 2122 Figure 21. Measured interannual wind speed anomalies for the windward side AWSs.

Figure 23 .

 23 Figure 23. Model interannual wind speed anomalies for the leeward side AWSs.

Figure 24 .

 24 Figure 24. Model interannual wind speed anomalies for the windward side AWSs.

Figure 25 .

 25 Figure 25. Model interannual wind speed anomalies for the Outer Islands AWSs.

  overall model performance is in good agreement with the ground-based measurements. The WRF mesoscale NWP model can successfully reproduce the observed wind climate on the ground. High wind resource areas are identified at higher altitude inland areas of the larger islands of Viti Levu and Vanua Levu, nearshore land areas on the eastern and the western sides of the larger islands and the nearby onshore areas adjacent to the channel between the two larger islands.Mesoscale wind resource maps of the long-term-mean wind speed and the long-term-mean wind power density at 55-meter-elevation can be compared with the NREL Wind Power Classification[41]. This reveals forty potential wind resource locations with wind speed and power density greater than 6.4 m/s and 300 W/m 2 (these values corresponds to Wind Power Class 3 for the 50 metre classification, which is considered feasible for utility-scale wind power applications). For the two main islands of Viti Levu and Vanua Levu where the national electricity grid is available, there are more than ( grid cells) of areas on each island where the wind speed and the power density are greater than the Wind Power Class 3 thresholds. Using the rule-of-thumb that of spatial area is required to site a 10 MW wind farm with wind turbines having 5-rotor

  available. The annual wind speed variability is low because Fiji is situated in the Southeast Trade winds zone, which are steady large-scale winds. The interannual variability is small as well (~0.8 m/s). The small natural variability of wind speed guarantees regular energy production. The El-Nino, La-Nina and the ENSO events are the main modes of interannual variability observed on the ground. These are also captured by the model and are the cause of crests and the troughs of the interannual wind speed anomalies in Figures 20 -25. According to the NASA Earth Observatory [46], El-Nino events occurred between 2009 -2010, 2014 -2016 and 2018 -2019, La Nina events occurred between 2008 -2009, 2010 -2012 and 2017 -2018 and ENSO events occurred between 2014 -2016.

  31 km x 31 km grid resolution every hour since 1979, this research created high-resolution wind resource data at 1 km x 1 km resolution for the recent decade (2009-2018) for the whole of Fiji, even in locations where ground-based measurements are not available. This data has been used for potential wind resource site identification and to evaluate the theoretical potential of wind energy for Fiji. It can also be used to evaluate the technical potential by coupling it with a microscale model. The advantage of the microscale model (WAsP) is that it uses the wind atlas and generalisation methodology, which removes all the terrain influence from the model data going upwards to the free stream wind speed (geostrophic wind speed) and works downwards taking into account high-resolution terrain data to create a generalised wind climate which can be used at a predicted site to evaluate the wind resources with greater accuracy then the input data from the model [41].

Table 1 .

 1 Details and position of AWSs.

		Latitude	Longitude	Elevation	WRF
	Station Name	(° S)	(° E)	(m a.s.l.)	Elevation (m a.s.l.)
	Keiyasi	-17.8795	177.7552	89.8	63.7
	Koro Island	-17.3450	179.4184	108.8	80.1
	Korolevu	-18.2129	177.7304	25.7	41.5
	Labasa	-16.4333	179.4000	8.5	13.1
	Levuka	-17.6403	178.8186	41.5	43.1
	Lomaivuna	-17.8714	178.3601	122.1	128.5
	Momi	-17.8952	177.2668	43.8	25.9
	Nadarivatu	-17.5676	177.9632	824.1	761.8
	Nadi	-17.7599	177.4448	20.7	14.9
	Nausori	-18.0464	178.5591	5.7	4.9
	Rakiraki	-17.3404	178.2214	8.1	3.4
	Rarawai	-17.5564	177.6814	9.3	17.0
	RKS Lovoni	-17.7260	178.5551	17.1	14.1
	Saqani	-16.4749	179.7089	30.0	43.2
	Seaqaqa	-16.4758	179.1578	101.8	109.5
	Sigatoka	-18.1422	177.5039	6.7	12.0
	Suva	-18.1475	178.4536	6.6	15.6
	Tokotoko	-18.2186	178.1700	4.9	6.2
	Udu	-16.1411	-179.9947	43.7	0.3
	Viwa	-17.1494	176.9117	2.0	1.0
	Vunisea	-19.0469	178.1654	31.9	86.3
	Wainikoro	-16.3044	179.5586	15.1	15.3
	Yaqara	-17.4330	177.9774	20.0	11.2
	Yasawa	-16.6984	177.5746	40.0	5.1

Table 2

 2

		. Physics Parameterization
	Parameterization Options
	Long-wave radiation	RRTMG
	Short-wave radiation	RRTMG
	Microphysics	WRF Single-Moment 6-class
	Surface Layer	MM5 Similarity
	Planetary Boundary Layer	Yonsei University
	Cumulus	New Tiedtke
	Land Surface Model	Noah

  3-year averaged wind speed diurnal cycle correlation between the model and the measurements for individual station locations. The mean bias for 3 years across all 24 AWSs is 0.62 m/s, which is a 20.2 % overestimate in comparison with the measured wind speed. The mean RMSE is 2.22 m/s, the mean STDE is 2.00 m/s and the mean diurnal cycle correlation is 0.81 for the locations tested.

Table 3 .

 3 Wind statistics from the leeward side AWSs(2016 -2018)

		Measured	WRF	Bias		RMSE	STDE	Correlation
	Station	Mean (m/s)	Mean (m/s)	(m/s) Bias %	(m/s)	(m/s)	(%)
	Keiyasi	1.15	1.86	0.71	61.9	1.51	1.34	0.99
	Labasa	1.57	2.94	1.37	87.3	1.99	1.44	0.91
	Momi	4.14	3.20	-0.94	-22.6	2.23	2.03	0.70
	Nadarivatu	2.81	3.54	0.73	25.9	2.00	1.87	0.58
	Nadi	2.77	2.79	0.02	0.7	1.68	1.68	0.89
	Rakiraki	6.10	5.53	-0.57	-9.4	2.30	2.23	0.95
	Rarawai	1.77	1.61	-0.17	-9.4	1.38	1.37	0.99
	Seaqaqa	1.36	2.20	0.83	61.2	1.78	1.57	0.97
	Sigatoka	1.89	2.41	0.53	27.8	1.38	1.28	0.98
	Wainikoro	1.65	3.49	1.84	111.1	2.61	1.85	0.94
	Yaqara	3.98	3.93	-0.05	-1.3	1.73	1.73	0.86
	Average	2.65	3.05	0.39	30.3	1.87	1.67	0.89

Table 4 .

 4 Wind statistics from the windward side AWSs(2016 -2018)

		Measured	WRF	Bias	Bias	RMSE	STDE	Correlation
	Station	Mean (m/s)	Mean (m/s)	(m/s)	%	(m/s)	(m/s)	(%)
	Korolevu	3.30	3.88	0.58	17.5	1.85	1.76	0.94
	Lomaivuna	1.98	2.61	0.63	31.8	1.65	1.52	0.97
	Nausori	1.94	3.08	1.14	58.9	2.59	2.32	0.90
	RKS Lovoni	2.06	3.95	1.89	91.8	2.67	1.89	0.10
	Saqani	2.80	4.17	1.36	48.7	4.15	3.92	0.47
	Suva	2.36	3.24	0.88	37.4	2.22	2.04	0.93
	Tokotoko	2.57	3.18	0.61	23.6	2.18	2.09	0.95
	Udu	5.36	4.80	-0.56	-10.5	1.96	1.88	0.72
	Average	2.80	3.61	0.82	37.4	2.41	2.18	0.75
		Table 5. Wind statistics from the outer islands' AWSs (2016 -2018)
		Measured	WRF	Bias		RMSE	STDE	Correlation
	Station	Mean (m/s)	Mean (m/s)	(m/s) Bias %	(m/s)	(m/s)	(%)
	Koro Island	3.94	4.47	0.53	13.4	2.00	1.93	0.69
	Levuka	3.82	5.15	1.33	34.7	2.79	2.46	0.85

Table 6 .

 6 Weibull Wind Distribution parameters for leeward side AWSs(2016 -2018)

		Measure	WRF			Measure	WRF
		d	Weibul			d	Weibul			Weibull
		Weibull	l	Weibul		Weibull	l	Weibul	Bias	Correlatio
	Station	(m/s)	(m/s)	l Bias Bias %			l Bias	%	n
	Keiyasi	1.17	2.02	0.85	72.5	1.17	1.21	0.05	4.0	0.95
	Labasa	1.70	3.35	1.65	97.5	1.35	2.69	1.34	99.5	0.60
	Momi	4.56	3.63	-0.93	-20.3	1.99	2.31	0.31	15.8	0.95
	Nadarivat
	u	3.15	4.07	0.92	29.5	1.71	2.17	0.46	26.6	0.90
	Nadi	3.03	3.17	0.14	6.0	1.65	2.33	0.68	41.3	0.95
	Rakiraki	7.11	6.21	-0.90	-12.6	2.44	2.42	-0.02	-0.7	0.97
	Rarawai	1.91	1.79	-0.12	-6.3	1.40	1.42	0.02	1.8	1.00
	Seaqaqa	1.48	2.44	0.96	64.8	1.23	1.42	0.19	15.8	0.94
	Sigatoka	2.11	2.69	0.58	27.6	1.58	1.47	-0.11	-7.1	0.98
	Wainikoro	1.88	3.99	2.11	112.6	1.40	2.28	0.89	63.7	0.59
	Yaqara	4.47	4.32	-0.15	-3.2	2.24	1.96	-0.28	-12.7	0.99
	Average	2.96	3.43	0.46	33.5	1.65	1.97	0.32	22.55	0.89

Table 7 .

 7 Weibull Wind Distribution parameters for windward side AWSs(2016 -2018)

		Measure	WRF			Measure	WRF
		d	Weibul			d	Weibul			Weibull
		Weibull	l	Weibul	Bias	Weibull	l	Weibul		Correlatio
	Station	(m/s)	(m/s)	l Bias	%			l Bias Bias %	n
	Korolevu	3.65	4.31	0.67	18.3	1.94	2.14	0.20	10.3	0.96

Table 8 .

 8 Weibull Wind Distribution parameters for outer islands' AWSs(2016 -2018)

		Measured	WRF			Measured	WRF			Weibull
		Weibull	Weibull	Weibull	Bias	Weibull	Weibull Weibull	Bias	PDF
	Station	(m/s)	(m/s)	Bias	%			Bias	%	Correlation
	Koro	4.47	5.02	0.55	12.4	2.06	2.57	0.50	24.7	0.95
	Levuka	4.30	5.78	1.48	35.2	1.75	2.85	1.10	63.2	0.76
	Vunisea	4.66	5.10	0.44	9.6	2.02	2.03	0.01	0.5	0.99
	Viwa	4.20	6.59	2.39	57.2	1.63	2.51	0.88	54.2	0.67
	Yasawa	7.19	6.64	-0.56	-7.6	2.61	2.88	0.27	10.3	0.98
	Average	4.97	5.83	0.86	21.4	2.01	2.57	0.55	30.6	0.87

Table 9 .

 9 Systematic Bias for Wind Power Density (WPD) for all AWSs at 10 m

	Station	Measured Mean (m/s)	WRF Mean (m/s)	D3 = 1 km Bias(m/s)	Measured WPD (W/m 2 )	WRF WPD (W/m 2 )	Bias (W/m 2 )
	Keiyasi	1.15	1.86	0.71	0.88	3.73	2.85
	Koro Is	3.94	4.47	0.53	35.38	51.67	16.29
	Korolevu	3.30	3.88	0.58	20.95	34.05	13.10
	Labasa	1.57	2.94	1.37	2.26	14.84	12.58
	Levuka	3.82	5.15	1.33	32.44	79.50	47.05
	Lomaivuna	1.98	2.61	0.63	4.49	10.28	5.79
	Momi	4.14	3.20	-0.94	41.30	19.07	-22.23
	Nadarivatu	2.81	3.54	0.73	12.03	24.04	12.02
	Nadi	2.77	2.79	0.02	12.39	12.66	0.27
	Nausori	1.94	3.08	1.14	4.26	17.06	12.80
	Rakiraki	6.10	5.53	-0.57	132.56	98.76	-33.80
	Rarawai	1.77	1.61	-0.16	3.24	2.44	-0.80
	RKS	2.06	3.95	1.89	5.10	35.96	30.86
	Saqani	2.80	4.17	1.37	12.80	42.27	29.48
	Seaqaqa	1.36	2.20	0.84	1.46	6.17	4.71
	Sigatoka	1.89	2.41	0.52	3.94	8.17	4.23
	Suva	2.36	3.24	0.88	7.68	19.86	12.19
	Tokotoko	2.57	3.18	0.61	9.91	18.78	8.87
	Udu	5.36	4.80	-0.56	89.62	64.36	-25.26
	Viwa	3.84	5.92	2.08	33.07	121.17	88.10
	Vunisea	4.02	4.58	0.56	37.84	55.96	18.12
	Wainikoro	1.65	3.49	1.84	2.62	24.80	22.18
	Yaqara	3.98	3.93	-0.05	36.79	35.42	-1.37
	Yasawa	6.27	5.87	-0.40	143.46	117.72	-25.74
	Average	3.06	3.68	0.62	28.60	38.28	9.68


David Tschumperlé

LIC-BASED REGULARIZATION OF MULTI-VALUED IMAGES

In this paper, a general multi-valued image regularization method based on LIC's (Line Integral Convolutions [4]) is proposed. From the investigation of recent approaches based on multi-valued diffusion PDE's, we show how a regularization process is naturally decomposed, first as the estimation of its underlying smoothing geometry, and then, as the application of a locally and spatially oriented smoothing. Performing this last part using LIC's significatively improves the overall regularization process both in visual quality and processing time. We illustrate three different applications of our general regularization framework : Color image denoising, inpainting and magnification.

INTRODUCTION

Obtaining regularized versions of noisy or scratched data has always been a desirable goal in the fields of computer vision and image processing. It is useful, either to restore images corrupted by noise (which is the most direct application of image regularization) or -more indirectly -as a pre-processing step that eases further analysis of the considered data. Since the pioneering work of Perona-Malik [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], anisotropic diffusion PDE's (Partial Differential Equations) raised a strong interest for this purpose : such equations have the ability to smooth data in a nonlinear way, allowing the preservation of important image features (contours, corners or other discontinuities). Thus, many diffusion PDE's have been proposed so far for the restoration of scalar and multi-valued images ( [START_REF] Aubert | Mathematical Problems in Image Processing: PDE's and the Calculus of Variations[END_REF][START_REF] Kimmel | Images as embedded maps and minimal surfaces: movies, color, texture, and volumetric medical images[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Sapiro | Geometric Partial Differential Equations and Image Analysis[END_REF][START_REF] Sapiro | Anisotropic diffusion of multi-valued images with applications to color filtering[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] and references therein). Despite this wide range of existing formalisms, all methods have something in common : they locally smooth the image along one or several directions that are different at each image point. Typically, the principal smoothing direction is often parallel to the contours within the image, resulting then in an anisotropic regularization process. Thus, defining a correct smoothing behavior is one of the key point of a good regularization algorithm. Recently, authors of [START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] proposed general frameworks able to design a specific regularization process from a given underlying local smoothing geometry : one first retrieves the geometry of the structures inside the image (generally by the computation of the so-called structure tensor field). Then, a local geometry of the desired smoothing is defined by the mean of a second field of diffusion tensors (depending on the first one). Finally, one step of the smoothing process itself is performed through a specific diffusion PDE. In this paper, we first review these efficient regularization methods. Then, keeping in mind the idea of separating the regularization process from the design of its smoothing geometry, we propose a new LIC-based framework that regularizes a multi-valued image according to a defined local smoothing geometry. Our method has two main advantages compared to classical PDE implementations : In one hand, it better preserves the orientations of small image structures. In the other hand, it runs up to three times faster. Finally, we illustrate the effectiveness of our generic LICbased regularization method, with results on color image restoration, inpainting and non-linear magnification, among all possible applications.

PDE-BASED REGULARIZATION

Let us consider a corrupted multi-valued image I : Ω → R n (n = 3 for color images) defined on a domain Ω ⊂ R 2 . We denote by I i : Ω → R, the particular channel i of the image I : ∀X ∈ Ω, I (X) = I 1(X) I 2(X) ... I n(X)

T .

Regularization algorithms in [START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] consist in computing firstly the smoothed structure tensor field G σ = G * G σ , where G σ is a 2D gaussian kernel (with a variance σ) and G : Ω → P(2) is the field of the symmetric and positive matrices defined as : ∀X, G (X) = n i=1 ∇I i(X) ∇I T i(X) . As noticed in [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF], G σ(X) is a good estimator of the local multi-valued geometry of I at X : its spectral elements give at the same time the vector (color) variations (eigenvalues λ 1 , λ 2 of G σ ) and the orientation (edges) of the local structures (eigenvectors u⊥v of G σ ) of I at each point X ∈ Ω. The convolution of G by G σ allows the estimation of a more coherent multi-valued geometry. Starting from this simple geometric measure G σ of the image, authors of [START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] design then a tensor field T : Ω → P(2) which defines the desired local smoothing behavior of the regularization process :

∀X, T (X) = f 1 (λ 1 , λ 2 ) uu T + f 2 (λ 1 , λ 2 ) vv T (1)
where v corresponds to the principal eigenvector of T (X) .

Basically, f 1 and f 2 define the strengths of the desired smoothing along corresponding directions u and v, at X. The regularization process itself is then done by evolving a diffusion PDE, either ∂Ii ∂t = div(T∇I i ) [START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] or ∂Ii ∂t = trace(TH i ) [START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF] where H i stands for the hessian matrix of I i(X) . Comparisons between these two equations is discussed in [START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF] and is out of the scope of this paper. Nevertheless, we remind the key idea that the smoothing is performed after the design of a local smoothing geometry defined by a field T of diffusion tensors. Note also that T might be updated during the process (after one or several PDE iterations). Applying such diffusion equations is very time consuming : I is regularized little by little since PDE-based processing is very local, even with recent and fast implementations [START_REF] Weickert | Parallel Implementations of AOS Schemes : A Fast Way of Nonlinear Diffusion Filtering[END_REF].

LIC-BASED REGULARIZATION

In order to speed up the smoothing process itself, we rather propose to replace the application of diffusion PDE's by a specific LIC-based method. Nevertheless, the first step (i.e. the computation (1) of the tensor field T defining the local smoothing geometry) is preserved.

LIC and single direction smoothing

LIC's (Line Integral Convolutions) have been first introduced in [START_REF] Cabral | Imaging vector fields using line integral convolution[END_REF] as a technique to create a textured representation of a vector field w : Ω → R 2 . The general idea, originally expressed in a discrete form, was to smooth an image I noise -containing only noise -by averaging the pixel values along the integral lines of w for all points X ∈ Ω. A continuous formulation of a LIC is then :

I LIC (X) = 1 N σ -σ f (a) I noise (C (X,a) ) da (2)
where f : R → R is an even function (decreasing on R + , typically a 1D gaussian) and C : Ω × R → Ω defines the integral curve of w starting from X and parameterized by a ∈ R, such that :

C (X,0) = X ∂C ∂a (X, a) = w(C (X,a) )
The normalization factor is N = σ -σ f (a) da. The σ parameter controls the smoothing length (size of the neighborhood used on the integral line by the LIC filter). Intuitively, eq.( 2) actually smoothes the image along the vector field w. This vector flow visualization problem has been also recently tackled by PDE-based methods in [START_REF] Becker | PDE methods in flow simulation post processing[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF]. It is basically based on the regularization methods described in section 2, with a field T of diffusion tensors defined as : ∀X ∈ Ω, T (X) = w (X) w T (X) .

Note that these smoothing methods are not strictly equivalent, but they are based on the same rough idea : smoothing an image along given directions defined by a vector field w.

We investigate the theoretical link between diffusion PDE's and LIC-based methods in [START_REF]Tschumperlé Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's[END_REF].

Multi-directional smoothing

The similarity between these two techniques performing single direction smoothing naturally suggests to replace the smoothing phase of general regularization PDE algorithms by a LIC-filtering process. Naturally, this implies that we have to extend the LIC principle (that considers only one single direction) to deal with general tensor fields T : Ω → P(2) instead of vector fields w. We propose to do so by averaging multiples LIC's starting from the same point X, but following different paths.

Let us denote by U θ the vector U θ = (cos θ sin θ) T . Then, the vector w θ (X) = T (X) U θ verifies :

• If T (X) is isotropic then w θ (X) = αU θ . • If T (X) is anisotropic and directed along U θ , then w θ (X) αU θ . • If T (X) is anisotropic and orthogonal to U θ , then w θ (X)
0. That can be understood as follows : the more U θ represents a part of T, the higher will be the norm w θ (X) . For a given θ, we propose to compute a LIC-filtering starting from X and directed by the vector field w θ , then we average all these "atomic" LIC filters for all θ ∈ [0, π]. This range is sufficient to reach the entire plane since LIC's are computed backward and forward along the integral lines C θ . Finally, our regularization equation, computing an anisotropically smoothed version I regul of an image I noisy is : ∀X,

I regul (X) = 1 N π 0 σ w θ (X) -σ w θ (X) f (a) I noisy (C θ (X,a) ) da dθ
(3) where the normalization factor is N = f (a) da dθ,

w θ (X) = T (X) U (θ) and f (a) = exp - a 2 2σ
σ is a user-defined parameter setting the global strength of the smoothing filter. C θ : Ω×R → Ω represents the integral line of the field w θ starting from X and parameterized by a :

C θ (X,0) = X ∂C θ ∂a (X, a) = w θ (C θ (X,a) ) = T(C θ (X,a) ) U (θ)
It is easy to verify that :

• The equation (3) respects the extremum principle, since only pixel averaging is performed.

• For points X ∈ Ω where T (X) is isotropic, the smoothing is also isotropic (averaging along all directions of the plane), since ∀θ, w θ (X) is constant. • For points X where T (X) is highly anisotropic, the smoothing is done only along the principal eigenvector of T (X) . These properties are also verified by most of the proposed anisotropic regularization PDE's, and ensure that the smoothing process is indeed driven by the tensor field T.

• The algorithme is very stable, even with very large time step dt, since only pixel averaging is performed. Moreover, the rotational invariance of the scheme is ensured from a numerical point of view.

APPLICATIONS

The Fig. 1 illustrates different applications of our LIC-based regularization algorithm on color images. Processing time is mentionned on the figure captions. Experiments have been done on a Sun Sparc 1.5 Ghz. Results are obtained quite fastly compared to PDE techniques. Actually, one often needs only one iteration to achieve the regularization since the LIC-based scheme is stable by construction. For our experiments, we defined T with eq.( 1) and

f 1/2 (λ 1 , λ 2 ) = 1 ( +λ1+λ2) p 1/2 (4)
where = 0.01, p 1 = 0.5 and p 2 = 2 (for denoising purpose) or p 2 = 10 (for inpainting and magnification purposes). The discretization is done with classical finite differences to compute the smoothed structure tensor field G σ (updated at each iteration), and discrete sums to approximate the two integrals of (3). The discretization of θ is quite large (dθ = 30 o ) while a is more precisely approximated (da = 0.2). We applied eq.( 3) for :

• Color image denoising (Fig. 1a,b,c) : Only one iteration of our algorithm has been necessary to obtain the results. Note how the noise is removed while the variously oriented image structures are well preserved (rotational invariance property of the scheme).

• Color image inpainting (Fig. 1d) : For this application, we applied eq.( 3) only on the points inside a user-defined mask (a checkerboard mask in our case). 100 iterations have been needed to get the results.

• Non-linear interpolation for magnification (Fig. 1e) :

These two examples shows how our LIC-based regularization technique (3) can perform super-resolution. A linear interpolation of the small image is regularized while keeping the know points unchanged (similar to inpainting with a very sparse grid mask). You can get more results and the algorithm executable at : http://www.greyc.ensicaen.fr/˜dtschump/greycstoration/

CONCLUSION

We have proposed a very generic and efficient regularization algorithm acting on multi-valued images. It gathers the best of classical filtering techniques used in the field of diffusion PDE's and LIC's. This allows us to propose a new regularization framework that improve result quality and processing time. Our LIC-based method is simple to implement and potentially handles a very wide range of image processing applications.

  (a) Removing grain noise from a 512 × 512 color image. (detail, 1 iteration, processing time : 19 seconds) (b) Removing pattern noise from a 555 × 367 color image. (detail, 1 iteration, processing time: 11 seconds) (c) Removing JPG artefacts from a 300 × 300 color image. (1 iteration, processing time : 13 seconds) (d) 290 × 290 color image inpainting. (100 iterations, processing time : 1 minute 26 s) (e) Non-linear image magnification (1 iteration, processing time : respectively 20 seconds) From left to right : original color image, nearest-neighbor interpolation, bicubic interpolation, our LIC-based interpolation (3).

Fig. 1 .

 1 Fig. 1. Results of our LIC-filtering method applied to various image processing issues.