Laurent Menut
email: menut@lmd.polytechnique.fr

Modeling of mineral dust emissions with a Weibull wind speed distribution including subgrid scale orography variance

The modeling of mineral dust emissons requires a fine knowledge of the wind speed close to the surface. In regional and global models, Weibull distributions are often used to better represent the subgrid scale variability of the wind speed. This distribution mainly depends on a k parameter, itself currently parameterized as a function of the wind speed value. In this study, we propose to add the potential impact of the orography variance in the wind speed distribution by changing the k parameter value. Academic test cases are designed to estimate the parameters of the scheme. A realistic test case is performed over a large domain ecompassing the northern part of Africa and Europe and for the period of 1st January to 1st May 2012. The simulations's results are compared to PM 10 surface concentations and AERONET Aerosol Optical Depth and Aerosol Size Distribution. We show that with the orography variance, the simulations results are closr to the ones without, showing this additional variability is not the main driver of possible errors in mineral dust modeling.

Introduction

For many geophysical studies, the wind speed close to the surface is a key parameter. This includes the analysis and modeling of mineral dust emissions. Primarily, the understanding of these emissions was performed using in situ or wind tunnel experiments, [START_REF] Gillette | Modeling dust emissions caused by wind erosion[END_REF][START_REF] Shao | Effect of saltation bombardment on the entrainment of dust by wind[END_REF]. This means that the deduced parameterizations were first mainly representative of a local spatial area. Using these experimental studies, Dust Production Models (DPM) were designed and later implemented in regional and global transport models. The transition between the local scale and the large scale modeled fluxes is ensured by using tuning parameters. For example, [START_REF] Tegen | Modeling of mineral dust in the atmosphere: Sources, transport, and optical thickness[END_REF] uses a simple relationship for the vertical flux estimation, including a scaling factor, [START_REF] Marticorena | Modeling the atmospheric dust cycle: 2 Simulations of Saharan dust sources[END_REF] linked the horizontal F h (saltation) to the vertical F v (sandblasting) flux using constant α factors as α = F v /F h (with a dependency on the soil characteristics only). Other well-known and used DPMs such as [START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distribution in source areas[END_REF], [START_REF] Shao | A model for mineral dust erosion[END_REF] or Kok et al. (2012), among others, are more realistic and calculate the vertical dust flux using more parameters. But even if they are more complex, they are also using tuning factors. Once the DPM scheme is included in transport models, simulations are most often validated with derived budget but not direct measurements of emissions (as this is not possible to measure emissions out of a wind tunnel). The most employed proxy is the Aerosol Optical depth (AOD) derived from satellite data or photometers from the AERONET network, [START_REF] Holben | An emerging groundbased aerosol climatology: Aerosol Optical Depth from AERONET[END_REF]). More recently, other derived parameters are used such as SEVIRI to estimate the hourly frequency of emissions, [START_REF] Schepanski | A new Saharan dust source activation frequency map derived from MSG-SEVIRI IR-channels[END_REF]). These AOD measurements are the results of vertically integrated concentrations, converted using uncertain optical properties of aerosols. It is therefore difficult to directly use these data to better understand or constrain the surface emissions process. This explains that the DPM may vary a lot between models, leading to a large variability of the mineral dust emissions estimates, [START_REF] Huneeus | Global dust model intercomparison in AeroCom phase I[END_REF]Cuevas et al. 2015).

To fill the gap between locally estimated emissions and a regional scale model, recent studies explored the way to describe the DPM's input parameters with distributions in place of single mean values [START_REF] Grini | Roles of saltation, sandblasting, and wind speed variability on mineral dust aerosol size distribution during the Puerto Rican Dust Experiment (PRIDE)[END_REF][START_REF] Pryor | Empirical downscaling of wind speed probability distributions[END_REF][START_REF] Menut | Sensitivity of hourly Saharan dust emissions to NCEP and ECMWF modelled wind speed[END_REF]). This corresponds to the way to represent variables also with their subgrid scale variability. This variability may be taken into account for many parameters in the DPM: recently, Menut et al. (2013b) use spatially high resolution databases of roughness lengths, soil and surface parameters to represent the existing variability in a model grid cell representing several tens of squared kilometers. This is also essential for the wind speed used in the emissions calculation. The mean wind speed used in a model often represents averaged values both in the temporal and spatial dimensions. But, as the mineral dust emissions in a threshold problem, the use of a single mean value may induce large modeling errors. For example, for particular soil conditions, emissions of dust may appear only up to 8.5 m s -1 : if the "real" mean wind speed is 8 m s -1 , there is no emission at all. But if the "modeled" wind speed is 9 m s -1 , emissions are calculated. This model error is possible since the wind speed model uncertainty is ≈ 1 m s 1 in regional models (see [START_REF] Gómez-Navarro | Sensitivity of the WRF model to PBL parametrisations and nesting techniques: evaluation of wind storms over complex terrain[END_REF] among others). To avoid these possible large errors due to the conjuction of wind speed uncertainty and threshold calculation using this wind speed, DPM use wind speed distributions such as the Weibull distribution, [START_REF] Weibull | A statistical distribution function of wide applicability[END_REF]. These distributions represents the fact that the mean wind speed may be biased or uncertain, but is also a way to represent the fact that during one hour and for a large grid cell, the instantaneous wind speed may vary a lot.

In this study, the main goal is to improve the subgrid scale variability of the mean wind speed by adding additional information in the k parameter estimation. This additional information is the orography variance of a model grid cell. The choice to add this specific variability was done because this information is available but not used in the model and because it seems logical that the surface wind speed is not spatially constant over surfaces of tens of squared kilometers. Then, in a model, when using a "mean" wind speed, it is clear that this value does not reflect completely all possible values observed in a grid cell where the orography may vary.

By adding a relationship between k and the orography, we add realism to the simulation and we also expect to reduce the uncertainty linked to the model resolution (the more realistic subgrid scale variability is taken into account the less the results depend on the horizontal resolution). Using the WRF and CHIMERE models, simulations are done with and without the modification of k. Using comparisons to measurements, we quantify if the modeled results are improved.

First, the data and tools used in this study are presented with the observations in Section 2 and the models in Section 3. The Weibull distribution is presented in Section 4 and the proposed modified formulation of the shape parameter k in Section 5. Results with an academic test case are presented in Section 6. Then, in Section 7, the impact of this additional term is evaluated on a realistic test case representing the simulation of the period January-April 2012 over a large domain encompassing the northern Africa (including Sahara and Sahel) and the Mediterranean region. Conclusions are presented in Section 8.

The observations

In this study, an academic and a real test case are presented. For the real test case, several observations are used: (i) the surface concentrations of PM 10 as measured by the Sahelian Dust transect (SDT), (ii) the AERONET photometers measurements for the Aerosol Optical Depth (AOD), the Angstrom exponent (Ae) and the Aerosol Size Distribution (ASD). Finally, the statistical scores used to quantify the results are presented in this section.

a. Surface concentrations of PM 10

For the surface concentrations of PM 10 , we use the Sahelian Dust Transect (SDT) measurements. These measurements are performed in Banizoumbou (Niger), Cinzana (Mali), M'Bour and Bambey (Senegal) as described in [START_REF] Marticorena | Temporal variability of mineral dust concentrations over West Africa: analyses of a pluriannual monitoring from the AMMA Sahelian Dust Transect[END_REF]. The concentrations are measured using a Tapered Element Oscillating Microbalance with a PM 10 inlet. These microbalances are at the same place that the photometers of the AERONET network.

b. Aerosol Optical Depth

The stations are selected to be representative of different locations, then different proximities to mineral dust sources. Banizoumbou and Cinzana are located in the center of western Africa and also correspond to stations of the SDT. They are considered as locations close to the sources. Dakar is located close to the coast in the western side of west Africa and Izana (Tenerife Island) is located in the Canary Islands. Dakar and Izana are more representative of stations 'under mineral dust plumes' after long-range transport.

The aerosol optical properties are compared between observations and model using the AERONET (AErosol RObotic NETwork) measurements [START_REF] Holben | An emerging groundbased aerosol climatology: Aerosol Optical Depth from AERONET[END_REF]. First, the comparison is done using the Aerosol Optical Depth (AOD) and for a wavelength of λ =550nm. Second, the comparison is performed using the Aerosol Size Distribution (ASD) estimated after inversion of the photometers data as described in [START_REF] Dubovik | A flexible inversion algorithm for retrieval of aerosol optical properties from Sun and sky radiance measurements[END_REF]. The level 2 data are used. For each AERONET station used in this study and listed in Table 1, the inversion algorithm provides volume particle size distribution for 15 bins, logarithmically distributed for radius between 0.05 to 15 µm.

c. Statistical scores

The statistics are calculated for correlation (R t ), Root Mean Square Error (RMSE) and bias. The correlation used in this study is the Pearson's correlation. Each correlation provides specific information on the quality of the simulation. The temporal correlation, noted R t , is estimated station by station and using daily averaged data in order to have homogeneous comparisons between all variables. This correlation is directly related to the variability from day to day, for each station. O t,i and M t,i represent the observed and modeled values, respectively, at time t and for the station i, for a total of T days and I stations. The mean time averaged value X i is:

X i = 1 T T ∑ t=1 X t,i (1)
The temporal correlation R t,i for each station i is calculated as:

R t,i = ∑ T t=1 (M t,i -M i ) (O t,i -O i ) ∑ T t=1 (M t,i -M i ) 2 ∑ T t=1 (O t,i -O i ) 2
(2)

The mean temporal correlation, R t , used in this study is thus:

R t = 1 I I ∑ i=1 R t,i (3)
with I the total number of stations. The normalized Root Mean Square Error is expressed as:

RMSE = 1 T 1 I T ∑ t=1 I ∑ i=1 O t,i -M t,i O t,i 2 (4)
for all stations i and all times t.

The models

For the results calculations, several models are used. First, we describe the Dust Production Model. Second, and for the real test case, we describe the Weather and Research Forecasting (WRF) model which calculates the meteorological fields and the CHIMERE chemistry-transport model which calculates the concentrations of mineral dust in the troposphere.

a. The Dust Production Model

The mineral dust flux is estimated following two steps. First, the saltation flux is calculated with the scheme of [START_REF] White | Encyclopedia of Fluid Mechanics[END_REF]:

F h (D p ) = ρ air g u 3 * 1 - u T * u * 1 + u T * u * 2 (5)
where u * is the "saltation" friction velocity, i.e the friction velocity calculated using the roughness length, z 0 . The roughness length is estimated using the global 6 km horizontal resolution 'Global Aeolian Roughness Lengths from ASCAT and PARASOL' dataset, [START_REF] Prigent | Comparison of satellite microwave backscattering (ASCAT) and visible/near-infrared reflectances (PARASOL) for the estimation of aeolian aerodynamic roughness length in arid and semi-arid regions[END_REF]. g is the acceleration of gravity, g = 9.81ms -2 , ρ air is the air density (depending on the meteorological model).

u T * is the threshold friction velocity depending on the soil particle diameter size D p and z 0 and is calculated following the parameterization of [START_REF] Shao | A simple expression for wind erosion threshold friction velocity[END_REF]:

u T * (D p ) = a n ρ p gD p ρ air + γ ρ air D p (6)
with the constant parameters a n = 0.0123 and γ=300 kg.m -2 . The particle density, ρ p =2.65 10 3 kg.m -3 is chosen to be representative of quartz grains clay minerals. The saltation flux is estimated only when u * > u T * (D p ) for a given soil particle diameter D p . One saltation flux is calculated for each grain size of the soil distribution, then all fluxes are integrated.

Second, the sandblasting flux is calculated following the Alfaro and Gomes (2001) scheme, modified in [START_REF] Menut | Sensitivity of the sandblasting fluxes calculations to the soil size distribution accuracy[END_REF]. The calculation is based on the partitioning between the dust cohesion energy and the kinetic energy of each saltating agregate. The emitted mineral dust is splitted into three aerosols modes (fine, coarse and big). These modes are represented using log-normal distributions with diameters d 1 =1.5 10 -6 m, d 2 =6.7 10 -6 m and d 3 =14.2 10 -6 m and their associated standard deviation, respectively σ 1 =1.7, σ 2 =1.6 and σ 3 =1.5. For each mode, a kinetic energy is defined and the flux is calculated as:

F v,m,i (D p ) = N class ∑ k=1 π 6 ρ p β p i (D p,k )d 3 m,i e i dF h (D p,k ) (7)
where N class is the number of intervals discretizing the soil size distribution in the range [D min p : D max p ] and d m,i the mean mass diameter. β is a constant as β =16300 and p i are factors depending on the soil size, after [START_REF] Alfaro | Modeling mineral aerosol production by wind erosion: Emission intensities and aerosol size distribution in source areas[END_REF]. One specificity of this scheme is that, in addition to the mean flux value, the emission size distribution depends on the wind speed. The higher the wind speed, the finer the emitted particles.

The surface and soil databases are described in Menut et al. (2013b). The soil characteristics are from the STATSGO-FAO database, [START_REF] Nachtergaele | Harmonized World Soil Database, FAO, Rome, Italy and IIASA[END_REF]), a global dataset with a native resolution of 30' × 30'. For the surface characteristics, the NCAR USGS database is used, with the same resolution than the soil database, [START_REF] Loveland | Development of a global land cover characteristics database and IGBP DISCover from 1 km AVHRR data[END_REF]. For each model grid cell, we thus have the relative percentage of all soil and surface characteristics. For each of these percentage, a sandblasting flux is calculated then cumulated to have the total sandblasting flux in each model grid cell.

Complementary to the soil and surface characteristics, an erobility factor is also necessary to estimate the relative part of the grid cell which could emit mineral dust. A global database was derived using MODIS satellite data to estimate this factor, as described in [START_REF] Beegum | Simulating Aerosols over Arabian Peninsula with CHIMERE: Sensitivity to soil, surface parameters and anthropogenic emission inventories[END_REF]. A map of erodibility is presented in Figure 1. This Figure also shows the orography variance (in m 2 ) provided by the WRF model database. This variance varies between 0 to 15000 m 2 on the domain. The most important variability is observed over mountainous areas, when the resolved topography corresponds to an average of very differents altitudes above sea level.

b. The WRF and CHIMERE model

The WRF and CHIMERE models are used to calculate the concentrations of aerosol in the troposphere. The system is an 'off-line' coupling model meaning that the meteorological fields are first calculated with WRF, then these meteorological fields are used by CHIMERE for the transport.

WRF is used in its 3.6.1 version, [START_REF] Skamarock | A Description of the Advanced Research WRF Version 2[END_REF]). This regional modeling uses NCEP global meteorological fields using a spectral nudging to ensure that the large-scale circulation is well taken into account. The complete set-up for this model is presented in [START_REF] Menut | Observations and regional modeling of aerosol optical properties, speciation and size distribution over Northern Africa and western Europe[END_REF].

CHIMERE is used in its version described in Menut et al. (2013a) and including the last updates presented in [START_REF] Mailler | CHIMERE-2017: from urban to hemispheric chemistry-transport modeling[END_REF]. For this study, the gaseous and aerosol chemistry is not used. The only aerosol taken into account is mineral dust: emissions, transport, mixing and wet and dry deposition. For all other parameters, we also used exactly the same configuration as in [START_REF] Menut | Observations and regional modeling of aerosol optical properties, speciation and size distribution over Northern Africa and western Europe[END_REF]. The modeled AOD is calculated by FastJX for several wavelengths, [START_REF] Wild | Fast-J: Accurate Simulation of In-and Below-Cloud Photolysis in Tropospheric Chemical Models[END_REF]. The boundary conditions for mineral dust are those of the climatology provided by GOCART, [START_REF] Ginoux | Sources and distributions of dust aerosols simulated with the GOCART model[END_REF]) and described in (Menut et al. 2013a).

In order to quantify the impact of our changes during a long period, the simulations ranges from 1st January to 1st May 2012. These four months correspond to the dry period in western Africa, when major mineral dust events occur. The modeled domain has an horizontal resolution of 60 × 60km and represents a large domain encompassing the northern part of Africa and the Mediterranean sea, as presented in Figure 1. This domain is strictly the same for the WRF and the CHIMERE models. The horizontal resolution was selected for different reasons: (i) this study is about subgrid scale variability: with this resolution of tens of squared kilometers, the variability is non negligible, (ii) for mineral dust studies, it is important to have a large modelled domain, encompassing all possible sources and long-range transport. This simulations ranging 4 months and over a large domain, this resolution make the computational cost manageable. Note that this resolution was used in previous studies, such as [START_REF] Menut | Observations and regional modeling of aerosol optical properties, speciation and size distribution over Northern Africa and western Europe[END_REF], [START_REF] Briant | Aerosol-radiation interaction modelling using online coupling between the WRF 3.7.1 meteorological model and the CHIMERE 2016 chemistry-transport model, through the OASIS3-MCT coupler[END_REF], [START_REF] Menut | An alternative way to evaluate chemistrytransport model variability[END_REF] and shows satisfactorily results compared to observations.

The main goal of this study is to change the wind speed. This wind speed is calculated by WRF and then used by CHIMERE for several processes: the transport, the calculation of the friction velocity (for the dry deposition) and the mineral dust emissions. The use of the Weibull distribution is made only in CHIMERE and only for the mineral dust emissions. For all other processes, the mean wind speed remains unchanged and as it is provided by WRF.

The Weibull distribution

This section presents the main concept of the Weibull distribution, the parameters to define and the way to estimate them.

a. The Weibull distribution formulation

The Weibull distribution is expressed as the following probability density function, following [START_REF] Babiker | A statistical model for horizontal mass flux of erodible soil[END_REF]:

p(|U sg |) = k λ |U sg | λ k-1 exp - |U sg | λ k (8)
where |U sg | is the subgrid-scale wind speed value, p(|U sg |) is its probability density function, k is a dimensionless shape parameter and λ (m s -1 ) is a scale parameter related to the mean of the distribution (in general the mean wind speed value, |U|).

The subgrid-scale wind speed |U sg | is estimated from the grid-scale mean wind speed |U| as:

|U sg | = 2i N w |U| (9)
where i ranges from 1 to N w , the number of subgrid-scale wind speed values defined to estimate accurately the distribution.

The λ value is estimated with the k parameter, which enables to reduce the Weibull distribution uncertainty to the k parameter only. The expression of λ is:

λ = U Γ 1 + 1 k (10)
with Γ the gamma function. The main problem for an accurate calculation of the Weibull distribution is thus focused on the accuracy of the k parameter.

b. Current estimate of k

The sensitivity to k is a well-known problem. Numerous studies were dedicated to its study and we present here only a few examples. Using measurements made at one site, [START_REF] Babiker | A statistical model for horizontal mass flux of erodible soil[END_REF] made a sensitivity study to quantify the impact of the mineral dust emissions fluxes to several k values ranging from 1 to 3. Also using measurements (30 years and 1432 differents measurements sites in US), [START_REF] Gillette | Modeling dust emissions caused by wind erosion[END_REF] estimated k values between 0.5 and 3.86. [START_REF] Christofferson | A simple estimator of the shape factor of the two-parameter Weibull distribution[END_REF] proposed a methodology to estimate the k value from numerous wind speed observations. [START_REF] Pérez | Weibull wind speed distribution: Numerical considerations and use with sodar data[END_REF]) used sodar data to retrieve wind speed values close to the surface: in this case, several k values are tested and compared to the observations, in a range from 1 to 4. [START_REF] Kelly | Weibull-k revisited: "Tall" profiles and height variation of wind statistics[END_REF] studied the potential variability of k and presented results with values ranging between 1.5 and 3. [START_REF] Gryning | Weibull Wind-Speed Distribution Parameters Derived from a Combination ofWind-Lidar and Tall-Mast Measurements Over Land[END_REF] presented a similar technique, but using mast and lidar measurements and concluded that k varies between 1.8 and 3.1 (depending on the altitude and the sites localization). [START_REF] He | Probability distributions of land surface wind speeds over North America[END_REF]) used wind data from 720 stations (from 1979 to 1999) and showed that the k parameter could depend on the atmospheric stability: the Weibull distribution is narrower during the night, due to intermittent turbulence. Statistically, they showed that the k value could depend on the wind speed value with a larger distribution for weak winds (high wind variability) and narrower for strong winds (well established wind speed regimes). The first studies showed that a constant value of k is not realistic and this parameter has to be variable in space and time.

More specifically for the mineral dust emissions problem, some other studies were already dedicated to the best as possible estimate of k. [START_REF] Menut | Sensitivity of hourly Saharan dust emissions to NCEP and ECMWF modelled wind speed[END_REF] used a constant value of k, limiting the variability. IIt was done to ensure a certain stability of the calculation in case of operational forecast calculations (i.e no day to day tuning of the model). This approach has to be replaced in models in order to be more realistic. [START_REF] Ridley | Toward resolution-independent emissions in global models: Impacts on the seasonal and spatial distribution of dust[END_REF]) made two simulations for the same domain and the same period: one with an high resolution to explicitely calculate the wind speed variability. The second one, with a coarsest resolution, to test the deduced k parameters used with a Weibull distribution. The approach is realistic but valid mainly for one studied case. A change in the period of the domain requires to replay the two simulations, fine and coarse.

A relationship between k, the mean wind speed |U| and its variance σ U was presented by [START_REF] Justus | Methods for estimating wind speed frequency distributions[END_REF] (hereafter called J78). This scheme is the only we found to express the k value as a function of a meteorological parameter. This is why it was selected for this study. It is expressed as:

k = U σ U 1.086 (11)
A direct relationhip linking k and σ U is used in [START_REF] Grini | Roles of saltation, sandblasting, and wind speed variability on mineral dust aerosol size distribution during the Puerto Rican Dust Experiment (PRIDE)[END_REF] as:

k = 0.94 U (12)
This expression does not take into account subgrid scale variability such as the orography. It considers that the Weibull shape mainly depends on the wind speed magnitude. This latter relation is used in many studies such as [START_REF] Su | Numerical simulations of Asian dust storms using a coupled climate-aerosol microphysical model[END_REF] with the CAM3 model applied over China, (Zhang et al. 2016) with CAM5 and at the global scale (2 o resolution) and [START_REF] Tegen | Modelling soil dust aerosol in the Bodélé depression during the BoDEx campaign[END_REF] for the comparison to observed wind values in Bodele. Note that in these studies, they are using a dust production model for which the size distribution is not dependent on the wind speed (unlike the model used in this study). • For |U|=8 m s -1 , the distribution ranges from 0 to 15 m s -1 . The J78 shape parameter value is k=2.66 • For |U|=14 m s -1 , the distribution ranges from 0 to 28 m s -1 . The J78 shape parameter value is k=3.52

Having a k value evolving with the mean wind speed value, it is interesting to notice that, for |U|=8 m s -1 , the J78 distribution has the largest width, when for |U|=14 m s -1 , the distribution has values between the distributions calculated with k=3 and 4.

Adding the orography variability in k a. The f k factor

The main goal of this study is to estimate a k value taken into account the orography subgrid variability. This quantity represents a factor directly affecting the surface wind speed and is currently not taken into account in the models using this wind speed. For that, we thus propose to add a variability term to the parameterization of J78. The new k depending on the orography is thus expressed as:

k = f k U σ U 1.086 (13)
where f k is a factor, using a simple function, respecting the following rules:

• First, we have to ensure that k will evolve around its mean value and within an acceptable range. We thus define two constants, f k,min and f k,max , designed to act as limiters. They represent the minimum and maximum possible factor around the mean k value. • We consider that over flat terrain, the wind speed module would have a lower variability that over complex terrain with a changing orography. The value of f k =1 represents a "mean" orography variance in the selected range of possible values of orography. For f k < 1, we have k decreasing, representing a more peakly wind speed distribution and thus an orography variance less important than the "mean" value: over flat terrains, the wind speed is more spatially homogeneous, because there is no orography. • The change is designed to act only over erodible surfaces since this is only applied for mineral dust emissions. It is thus necessary to find the maximum of orography variance until where the mineral dust emissions occur. The orography variance may be very large over mountainous areas but these areas are non erodible (being mainly constituted of rocks) and are thus not taken into account. These constraints lead us to define a simple formulation such as:

f k = f k,min + ( f k,max -f k,min ) ×     1 - c 1 + a exp -b σ z σ z,max     (14)
where: σ z is the orography variance of the model cell where the dust flux is calculated. σ z,max represents the maximum of orography variance for which we want to reduce the k value. The a, b and c values are just constants dedicated to modify the shape of the function. Here, we want to have a smooth function and we selected a=20, b=10, c=1. Some sensitivity tests were done and results were not sensitive to a change in these values. A schematic representation of f k is presented in Figure 3. For this study, we selected a f k decrease/increase of ±20%, leading to define f k,min =0.8 and f k,max =1.2. These values were arbitrarily selected but represent "reasonable" variability of the k parameter, regarding the numerous studies having fitted this parameter using wind speed measurements.

Academic test case modeling

In order to understand the variability of k as a function of the mean wind speed and various values of orography variances, a simple academic test case is performed. For this test case, the DPM of Alfaro and Gomes ( 2001) is used alone, in a 0D context (meaning that there is no time or space, but only a calculation using a varying academic wind speed value). In order to quantify the maximum of the impact of the f k factor, the Weibull distribution is calculated with three differents hypothesis and results are displayed in Figure 4: • f k =1 corresponding to the J78 scheme without any change • f k with σ z =0.01 × σ z,max corresponding to a low orography variability. In this case, we consider that the orography variance corresponds to 1% of the maximum orography variance. • f k with σ z =σ z,max corresponding to a maximum of orography variability (100%)

Results are first presented in Figure 4 (top) for the evolution of k depending on the mean wind speed value |U|. For f k =1, noted [J78], the k value regularly increases when |U| increases. For the two other cases, the increase has a similar trend, but the absolute values of k differ. For example, and for |U|=10 m s -1 , we have k=3 for the J78 scheme, but k ≈ 3.5 and 2.5 when the orography is low and high, respectively (Sigma=1% and Sigma=100%).

For the same three expressions of k, results are presented as Weibull distribution in Figure 4 (middle). The distribution with J78 without any change is inbetween the two others distribution when f k is applied. sponding integrated values are presented in Table 2. The flux calculated with the most important variability (then low k and Sigma=100%) show the highest values. In Table 2, this correponds to F v =60.16 kg m -1 day -1 when F v =53.91 kg m -1 day -1 when there is no variability due to orography. For a higher value of k, the distribution is narrower and this leads to less important fluxes:

F v =39.19 kg m -1 day -1 .
Emissions are diagnosed for three main modes as defined previously. Here, they are expressed in radius, to be later compared with AERONET measurements. The highest impact on the fluxes appear for the fine mode of the distribution, r p =0.75 µm. This mode corresponds to the optically efficient part of the size distribution. The change of k has thus an impact on the vertical flux but also on the deduced Aerosol Optical depth (AOD). As presented in Table 2, with no variability, AOD=0.92. If the orography variability is taken into account, AOD=0.54 and 1.10 for an orography variance of 1 and 100%, respectively. This huge impact, both on emissions fluxes and on AOD, shows the high sensitivity of the determination of the k parameter for mineral dust model studies.

Real test case

The two simulations, without and with the orography impact on the k parameter value, are called hereafter "J78" and "J78+Var oro", respectively. The simulations containing only mineral dust as aerosol, the results are compared to measurements performed in locations where mineral dust dominates the aerosol composition: the AERONET stations located in Western Africa and the surface PM 10 measurements of the Sahelian Transect.

a. Determination of the k parameter over the domain

The first step for the real test case is to calculate the parameters involved in equation 14. Among all parameters, it is necessary to select a correct value for σ z,max . To do this, we calculate the distribution of the SSO variance. This is calculated for the whole domain and for the erodible surfaces only. The result is displayed in Figure 5. The SSO variability between 0 and 500 m 2 represents occurence between 80 and 90%. Between 500 and 1000m 2 , it reduces to less than 10% for both cases. There is more low variabilty for the whole domain, probably due to the fact that ocean model cells are taken into account. The values higher than 1000 m 2 represent a few percents of the distribution. As we want to act mainly on the range where the variability is the most frequent over the erodible areas, we select σ z,max =1000 m 2 . The results of the k parameter calculation using these criteria is presented in Figure 6. Values are estimated for cells with an erodibility of 10% at least. For the other cells, the k parameter is defined in the model as described by J78, but is not used for mineral dust calculation, the cells being considered as non erodible. For a wind speed of 10 m s -1 , J78 scheme gives as result k=2.97 (as presented in Table 2). Using the orography variance, the k values are varying between 2.2 and 3.8. The lowest values are obtained where the orography is high: for example to the north of the Bodele region and in the Hoggar. The highest values are defined where the orography is low and these areas represents the majority of the erodible cells for the studied domain.

Globally, we probably will have mineral dust emissions less important than using only the J78 parameterization. But, we expect to have a different spatial variability of the emissions since the change in the k parameter depends on varying orography over the whole domain.

b. Results for surface concentrations of PM 10

The PM 10 surface concentrations are compared between the model and the observations and with an hourly time-step. Results are first presented as time series in Figure 7 and for the sites of Banizoumbou, Cinzana and Dakar. The observations, represented as symbols, shows a large time variability and values are ranging from 0 to 2000 µg m -3 .

The results mainly show the high temporal variability of the measurements and the fact that the model is able to reproduce this variability. The "J78+Var oro" configuration provides PM 10 lower that the "J78". This means that over the erodible regions leading to mineral dust measured at these sites, we have an orography variance less important than the mean average value and thus a increased k parameter, leading to a more peakly Weibull distribution and thus less mineral dust emissions.

The results are quantified as statistical scores in Table 3. The first part of this Table corresponds to the model results for "J78" and the second part to "J78+Var oro". For the correlation and the RMSE, the best statistical scores are bolded. For the correlation and the RMSE, the results are better for "J78+Var oro" than for "J78" (except for the correlation in Banizoumbou, but the difference is low and only 0.01). The most important benefit with the "J78+Var oro" configuration is for the RMSE. This is mainly because the shape parameter k is increased, leading to a sharper Weibull distribution and thus less variability in the mineral dust emissions. The model error is reduced by at least ≈ 50% for each of the three stations. This benefit being true for the three stations, and these three stations being not close but representative of different regions, we can expect that this "J78+Var oro" configuration improves the results in a realistic way.

c. Aerosol Optical Depth time series

Results for the AOD comparisons with the AERONET data and the CHIMERE model outputs are displayed in Figure 8 as time series, from the 1st January to 30 April 2012.

The time series show that several high peaks of mineral dust occured during the four studied months. For all these peaks, the model is able to retrieve correctly their timing (the correct day and the duration) and their magnitude. A few peaks are missing (as in January in Banizoumbou or mid-April in Izana) but they are very peaked and probably correspond to local events, difficult to catch with the used model resolution, such as convective cold pools or sporadic wind speed extremes. In general the model is able to estimate the background values and the major dust events. As for the PM 10 time series, the "J78+Var oro" configuration simulates less AOD than the "J78" one.

In order to quantify these results of AOD, statistical values are presented in Table 4. The simulation with the orography variance improves the model results compared to the observations: the RMSE is better for all stations, showing that the model error is globally reduced. For the correlation, the differences between the two simulations are less important, even if the best results are also mainly for the "J78+Var oro" version. The bias decreases with "J78+Var oro", highlighting the results showed with the time series: for these stations, the AOD is generally lower than for the "J78".

d. Aerosol Optical Depth maps

As a complement of the time series and the statistical scores, maps of AOD values is presented in Figure 9. The modeled period lasted four months and the maps represent four days during the period: 15th of Januray, February, March and April 2012 at 12:00 UTC. Note that other maps were analyzed and the discussion and conclusion are similar. The left column represents the AOD absolute values calculated with the "J78+Var oro" simulation and the right column correponds to the difference of AOD calculated with AOD(J78+Var oro)-AOD(J78). The AOD highest values are located in different places according to the day. But for the four days, the "J78+Var oro" corresponds always to a decrease in AOD, the largest differences being colocated with the highest AOD values. It shows that the use of the orography variance leads to an increase of k and thus a decrease of emissions, thus AOD. It means that the change is mainly applied over erodible areas where the orography variance is low.

For example and for the 15 April 2012, two main plumes are visible: a first one in the center of Africa (over Niger and Chad), and a second one in Guinea and Senegal and transported towards the Atlantic sea. In these plumes, the AOD values reach 2 at the maximum when the 'background values' over Africa are between 0 and 1. For the AOD 'background' values, the differences reach 0.2. In 4. Scores for the comparisons between observations (AERONET) and model (CHIMERE) for the Aerosol Optical Depth (AOD). Results are presented with N (%) the percentage of hourly mean available measurements for the period from 1st January to 1st May 2012, the mean values over the period (called "Obs" and "Mod", for the observations and modelled values, respectively), the temporal correlation (R t ), the Root Mean Squared Error (RMSE) and the bias (model minus observations). Values are bolded for the simulation having the highest correlation or the lowest RMSE.

the plumes, where the AOD are maximum, the differences reach 0.8, corresponding to an important change in the AOD estimation.

e. Aerosol Size distribution

As presented in Figure 4, a change of the k parameter acts on the mineral dust emissions flux magnitude but also on its size distribution. In order to quantify the impact of the k changes on the size distribution, we compare here the model results to the AERONET inversion. In the AERONET database, note that all AOD measurements are not inverted and we had to find available data. Two comparisons are presented but many other size distributions were analyzed and the conclusions are similar. The results presented correspond to available data close to the results showed in Figure 9, and corresponding to the 15 March and April 2012. Comparisons are presented in Figure 10.Note that, here, the distributions are expressed in radius (to be consistent with the raw data provided by AERONET).

For the four days and the four sites, the most important mode corresponds to the "coarse" mode with r p ≈ 3.35 µm. Less important, concentrations are visible for the "fine" mode, r p =0.8 µm, with the model but not in the observations. The simulation using "J78+Var oro" is closer to the observations than the "J78" simulation. If the "coarse" mode dominates the mass, the "fine" mode explains the most important differences between the two simulations. And for the calculation of AOD as compared with AERONET (with λ =550nm), this is the most sensitive part of the size distribution. These figures show that the variability of k acts on the correct part of the size distribution and proves that the RMSE is highly improved for correct reasons.

Conclusions

This study examined the impact of the wind speed Weibull distribution on the calculation of mineral dust emissions. This distribution depends mainly on one parameter, k, parameterized as a function of the mean wind speed. In this study, we consider that for regional models, having a grid cell of tens of kilometers, the subgrid scale orography may also be an important parameter to take into account in the estimation of k since the orography has a direct impact on the value of the mean wind speed close to the surface.

Using the information of subgrid orography variance, we proposed an adjustment of the estimation of the k parameter. We first test this change with academic cases in order to estimate realistic boundaries to the potential evolution of the k value. Second, we applied this change in the framework of a realistic modeling: four months, from 1st January to 1st May 2012 and over a large domain, encompassing Western Africa and Europe. The modeled results were compared to surface PM 10 concentrations and AERONET Aerosol Optical Depth and aerosol size distribution. Correlation, Root Mean Square Error and bias were also calculated. We showed that for a large part of the comparisons between observation and model results, the correlation is slightly improved (≈ +0.01) and the RMSE is greatly improved (≈ -30%) both for mass and AOD. Over the studied period, a systematic decrease of AOD and PM 10 is modeled with the changed k shape parameter. It means that, over erodible surfaces, k was increased. Considering the proposed function, it means we are in the range of 'low' variability of orography.

Even if the RMSE is better, the correlation is not really improved. It showed that the subgrid scale variability of orography is not the main driver for the mineral dust emissions. This is surprising because a significant impact was expected: the DPM schemes are very sensitive to the mean wind speed and the mean wind speed is sensitive to the orography close to the surface. The fact to have quantified this impact remains useful even if the impact intensity is lower than expected. For processes such as mineral dust emissions, the sub-grid scale variability remains a way to link a local phenomenon and an averaged model grid cell. It is a better approach than the model tuning, sometimes efficient, but not able to explain the physics and often limited to a specific model, with its own parameterizations and resolution. The model tuning is also systematic where the change proposed in this study varies in each model grid cell, adding realism to the system. A potential future direction could be the convective cold pools, knowed to increase mineral dust emissions very locally and certainly missing in many regional models.

FIG. 1 .

 1 FIG. 1. The domain defined for the real test case. (top) Map of the erodibility (0:1). (bottom) Map of the orography variance (m 2 ).

FIG. 2 .

 2 FIG. 2. Example of Weibull wind speed distribution for a mean value of |U|=8 (in black) and 14 (in blue) m s -1 . The three curves correspond to k=3, k=4 and the parameterization in Justus et al. (1978) (here J78), respectively. With J78, the corresponding k values are k=2.66 for |U|=8 m s -1 and k=3.52 for |U|=14 m s -1 .

•

  FIG. 3. The f k factor variability as a function of the orography variance.

FIG. 4 .

 4 FIG. 4. (top) Values of the k parameter values for the J78 parameterization and including the orographic weight factor proposed in this study. The two curves represent the extrema of the change, with f k,min and f k,max corresponding to the lowest and largest orography variances, respectively. (middle) Weibull distribution with and without taking into account orography variance. (bottom) Sensitivity of the mineral dust emission vertical flux, F v , as a function of the impact of the orography.

FIG. 5 .

 5 FIG. 5. Distribution (in %) of the orography variance, SSO variance (m 2 ), over the whole modeled domain and for the erodible cells only.

FIG. 6 .

 6 FIG. 6. 'k' parameter values using the orographic weight factor proposed in this study, and for |U|=10m s -1 . Values are estimated for cells with an erodibility of 10% at least. The frames indicate the main mineral dust sources in Africa. The surface stations are indicated with their names.

FIG. 7 .

 7 FIG. 7. Hourly PM 10 surface concentrations time-series from 1st January to 30 April 2012, for the model (with the two simulations without and with the orography variance) and the Sahelian Dust Transect sites: Banizoumbou, Cinzana and Dakar.

FIG. 8 .

 8 FIG. 8. Aerosol Optical Depth (λ =675nm) time series from 1st January to 30 April 2012, for the model and the AERONET sites. Observations are represented by symbols and the two plain curves represent the model simulation without and with orographic variability of the k parameter.

FIG. 9 .

 9 FIG. 9. Aerosol Optical Depth (λ =600nm) maps for the 15 January, 15 February, 15 and 15 April 2012, 12:00 UTC. (left) Absolute valuesof AOD for the "J78+Var oro" simulation and (right) AOD differences as "J78+Var oro"-"J78".

TABLE 1 .

 1 Names and locations of the AERONET stations used for model comparisons to AOD and ASD data. The stations are ordered from South to North. The altitude ASL are not presented, the measurements being representative of the vertically integrated atmospheric column Above Ground Level (AGL).

	Site	Longitude Latitude	Country
		( o W)	( o N)
	Ilorin	4.34	8.32	Nigeria
	Cinzana	-5.93	13.27	Mali
	Banizoumbou	2.66	13.54	Niger
	ZinderAirport	8.98	13.75	Niger
	Dakar	-16.95	14.39	Senegal
	CapoVerde	-22.93	16.73 CapoVerde
	Tamanrasset	5.53	22.79	Algeria
	Izana	-16.24	28.47	Tenerife
	Saada	-8.15	31.61	Algeria

TABLE 3 .

 3 Scores for the comparisons between observations (Sahelian Transect) and model (CHIMERE) for the surface concentrations of PM 10 . Results are presented with N (%) the percentage of hourly mean available measurements for the period from 1st January to 1st May 2012, the mean values over the period (called "Obs" and "Mod", for the observations and modelled values, respectively), the temporal correlation (R t ), the Root Mean Squared Error (RMSE) and the bias (model minus observations). Stations are sorted in increasing latitude, from south to north. Values are bolded for the simulation having the highest correlation or the lowest RMSE.

	Site	N	PM 10	R t	RMSE	bias
		(%)	Obs Mod
	k=f(J78)
	Banizoumbou	90.9 269.	249.	0.46	1.46	-20.
	Cinzana	99.9 293.	238. 0.33	2.12	-54.
	Dakar	99.4 256.	321. 0.52	4.23	64.
	k=f(J78,varSSO)
	Banizoumbou	90.9 269.	183. 0.45	1.06	-85.
	Cinzana	99.9 293.	177.	0.34	1.63	-116.
	Dakar	99.4 256.	228.	0.55	2.73	-27.

Acknowledgments. For the Sahelian Dust Transect data, we thank Bernadette Chatenet, the technical PI of the Sahelian stations from 2006 to 2012, Béatrice Marticorena and Jean-Louis Rajot, the scientific co-PIs, and the African technicians who manage the stations. We thank FIG. 10. Aerosol Size Distribution compared between the AERONET inversions and the CHIMERE modeled concentrations.

the principal investigators and staff for establishing and maintaining the AERONET sites used in this study.


Thi Mai

Trang Nguyen
email: trang.nguyen@lip6.fr

Lyes Hamidouche

Fabien Mathieu

Sébastien Monnet

Syphax Iskounen

SDN-based Wi-Fi Direct Clustering for Cloud Access in Campus Networks

Keywords: Wi-Fi Direct, software-defined networking, cloud network access

Mobile cloud is changing the way to enroll teaching activities in university campus. Lectures and lab sessions can be carried out directly from tablets in a class room by accessing a server in the cloud. In this paper, we address the problem of high density cloud access with wireless devices in campus networks. We propose to use Wi-Fi Direct clustering to solve the problem of Quality of Service (QoS) degradation when a high number of wireless devices want to access a content in the cloud at the same time. A centralized softwaredefined network controller is used in our proposed architecture to capture the network state and organize the Wi-Fi Direct groups. The optimized number of clusters can be calculated in function of the number of devices in the room. By simulations, we show that we can provide a better QoS in terms of download time and application's throughput by reducing the interference in this dense wireless network environment.

Introduction

Cloud communications and networks since their invention have shown an important impact on our society. Starting by online file storage services such as Google Drive and Dropbox, today's cloud networks provide users with numerous advanced services with mobile access. For the education sector, mobile cloud services are changing the way that we teach at universities. Tablets have become a common tool for users and can be given to students at the beginning of the lecture with slides and lecture notes. Demonstrations can be easily done by video streaming through a Wi-Fi access available in the auditorium. Teaching materials can be updated with less cost and effort compared to the printed version. For numerous teachings in computer science, lab sessions can be done over tablets via an Internet connection to a distant server where the experiments are really running. Small tablets also consume less energy than fixed PC, making the campus greener.

However, providing Quality-of-Service (QoS) in dense wireless networks is a challenge. Many wireless devices in a small geographic area can drastically degrade the overall network performances and decrease the user's Qualityof-Experience (QoE). Two main reasons are the interference generated by simultaneous transmissions and the heavy traffic load requested by all the devices through the access point. With the high density of devices set at a high transmission power to reach the access point in the auditorium, the number of transmission collisions at the Medium Access Control (MAC) layer is very high. Collided frames are erroneous and must be retransmitted by the transmitter with an additional waiting Back-off time. As a consequence, data transmission delay and packet loss ratio are very high in the network.

Wi-Fi Direct [START_REF] Camps-Mur | Device-to-Device Communications with WiFi Direct: Overview and Experimentation[END_REF], a new technology widely available in tablets and smartphones, allows users to communicate directly to each other. Two nearby Wi-Fi users can activate Wi-Fi Direct communications to exchange pictures. A group of users can form a Wi-Fi Direct group to communicate without the necessity of having an access point. Recognizing the high automation and flexibility offered by Wi-Fi Direct to form user groups, we propose in this paper to use the Wi-Fi Direct clustering to reduce the interference in dense Wi-Fi campus networks. By dividing a large number of Wi-Fi users into smaller groups, combined with frequency assignment and power control techniques, we can reduce the overall interference level in the network and increase network performances. For the Wi-Fi Direct group formation to be efficient, a centralized controller following the Software Defined Networking architecture is used for network control.

While the main research motivation of this paper is focused on campus networks, other use cases such as events or conference venues can also take advantage of Wi-Fi Direct-based D2D communications. This technology can reduce data dissemination time in large scale events. There are companies (e.g. Magency [2]) that are specialized in providing mobile technology-based solutions to animate events. During an event, the attendees receive devices that they use to participate in different activities. Depending on the needs of each event, the number of devices simultaneously connected to the network can go up to thousands. The size of the data used by devices can go up to more than 1 GB. Using traditional Wi-Fi networks poses serious problems of scalability. Even if some practical solutions can be used to meet user's needs such as oversizing the network infrastructure or pre-loading data onto the devices, they only partially solve the problem. In fact, these solutions cost additional resources in terms of hardware and time required. Furthermore, when a high number of access points are physically close to each other, network throughput can be drastically reduced due to wireless interferences.

The remainder of this paper is organized as follows. In section 2, we give an overview on the Wi-Fi Direct technology and SDN architecture. Section 3 presents our SDN-based architecture for Wi-Fi Direct clustering. Section 4 presents our implementation of Wi-Fi Direct in the OMNeT++ simulator and evaluates the performances of the proposed architecture. Finally, section 5 concludes the paper.

Background

Dense Wi-Fi networks

Since the introduction of Wi-Fi networks at the beginning of the 21st century, this technology has been widely deployed in university campus, event centers and hotspots. The density of Wi-Fi hosts and access points have also been significantly increasing during the last five years. Especially, with the coming of smartphones and tablets with Wi-Fi integrated, users are moving to use Wi-Fi as the default connectivity to access cloud-based applications. Wi-Fi is also a principle wireless technology to offload mobile traffic. In most of campus networks, wireless connectivities are totally based on Wi-Fi. Following the Cisco VNI report [START_REF]Cisco, Cisco Visual Networking Index: Global Mobile Data Traffic Forecast Update[END_REF], 60 percent of mobile traffic has been offloaded over Wi-Fi and small-cell networks in 2016. Wi-Fi traffic from both mobile devices and Wi-Fi only device contributes 42 percent of the total IP traffic in 2015 and will account for 49 percent of total IP traffic by 2020. Total public Wi-Fi hotspots will grow six-fold, from 94 million in 2016 to 541.6 million by 2021. By 2021, there will be 8.3 billion handheld or personal mobile devices all over the world.

In high density Wi-Fi networks, the Quality-of-Service is degraded. Maity and al. [START_REF] Maity | TCP Download Performance in Dense WiFi Scenarios: Analysis and Solution[END_REF] have shown that TCP download performance degrades significantly in a dense network, e.g. a network with 20-30 clients per access point. Especially, in a classroom scenario, a student may need up to 229 seconds to download a quiz file of 200KB and up to 478 seconds to complete the download of a reference file of 5MB. In a residential scenario, the number of clients per access point is small but there can be up to 500 access points per square kilometer area [START_REF] Baid | Understanding Channel Selection Dynamics in Dense Wi-Fi Networks[END_REF] resulting in high interference and poor network performances.

The solutions to dense Wi-Fi network performances have been proposed mainly in three directions -Device-to-Device (D2D) communications, power control and channel selection. D2D communications [START_REF]Green PI: Device2Device wearable networks[END_REF] allow direct communi-

GO

Client

Client Client

Fig. 1 Wi-Fi Direct mode cations between devices to reduce the traffic going through the access points. Power control and channel selection techniques [START_REF] Baid | Understanding Channel Selection Dynamics in Dense Wi-Fi Networks[END_REF] enable an efficient interference management.

Wi-Fi Direct mode

In the Wi-Fi networks, two communication modes have been defined in the IEEE 802.11 standard: the infrastructure mode and the ad-hoc mode. The infrastructure mode, the most widely deployed mode, is based on a network infrastructure composed of access points deployed by a network administrator. When a user wants to send data to another, the data is transmitted to the access point which will retransmit the data to the destination, even if the two users are close to each other. In the ad-hoc mode, there is no access point. Devices in the transmission range of each other can communicate directly as a single-hop communication.

Recently, the Wi-Fi Alliance has defined the Wi-Fi Direct mode to facilitate the D2D communications among nearby Wi-Fi users. In a Wi-Fi Direct group (Figure 1), one of the nearby devices is elected as Group Owner (GO) which plays the role of access point. The other devices are clients and communicate to each other via the GO. Wi-Fi Direct can be seen as a combination of the infrastructure mode and the ad-hoc mode because we still have an AP (the GO) like the infrastructure mode but this AP is also a client as in the ad-hoc mode.

A standard Wi-Fi Direct group formation (Figure 2) is composed of four steps: scanning, discovery, GO negotiation and provisioning [START_REF]Wi-Fi Simple Configuration Technical Specification[END_REF]. First, a device scans to verify if there is an existing group. If there is no beacon message sent by an existing GO, the device moves on to the second step to discover a peer. The device switches between scanning and listening states over randomly chosen channels until there are two devices sending the ProbRequest message and receiving the ProbResponse message over the same channel. The third step, GO negotiation, is a three-way handshake. Each peer communicates a GO Intent value which is randomly chosen between 0 and 14. The device with the greater GO Intent value becomes GO. In the last step, the peers exchange authentication messages to secure the communications.

Wi-Fi Direct is a strong candidate to offload cellular traffic [START_REF] Pyattaev | 3GPP LTE traffic offloading onto WiFi Direct[END_REF]. The ability of quickly creating groups and disseminating information of Wi-Fi Direct can be used in mission-critical communications for public safety [START_REF] Rajadurai | Enhanced Interworking of LTE and Wi-Fi Direct for Public Safety[END_REF]. Proactive routing protocols can be used with Wi-Fi Direct to build Mobile Ad-hoc Networks (MANET) [START_REF] Liu | Development of Mobile Adhoc Networks over Wi-Fi Direct with Off-the-Shelf Android Phones[END_REF]. An open-source software has been also developed with Android devices for opportunistic data transmission with Wi-Fi Direct [START_REF] Amaral | Oi! -Opportunistic Data Transmission based on Wi-Fi Direct[END_REF].

2.3 Software-Defined Networking in Wi-Fi Networks SDN is a new concept which separates the network control from data forwarding (Figure 3). The data plane includes network devices which forward user data. The control plane is composed of one or more SDN controllers which have a consistent and global view to control the network. The SDN controllers communicate with the APs via the an SDN protocol such as OpenFlow [START_REF] Mckeown | OpenFlow: Enabling Innovation in Campus Networks[END_REF]. In Wi-Fi networks, the centralized view of the SDN controller can efficiently optimize the complex wireless environment. In [START_REF] Riggio | Interference Management in Software-Defined Mobile Networks[END_REF], SDN has been used for interference management. The SDN abstractions provide a global view of channel quality and interference in the network and allow quick actions on reconfiguration such as channel assignment. The authors in [START_REF] Stiti | Virtual OpenFlow-Based SDN Wi-Fi Access Point[END_REF] proposed an architecture using SDN for the management of virtual Wi-Fi access points. Flow tables in virtual access points are reprogrammed when users move in the network to provide seamless handovers. In [START_REF] Seyedebrahimi | SDN-Based Channel Assignment Algorithm for Intererence Mangement in Dense Wi-Fi Networks[END_REF], SDN has been used for channel assignment in dense Wi-Fi networks. Channel assignment problem is modeled as an optimization problem and implemented in the SDN controller to react when the interference level in the network gets higher than a predefined threshold. The centralized wireless management enabled by SDN can improve up to 12 Mbps of channel capacity measured at physical layer compared with the uncoordinated one.

When applied to wireless networks, the SDN concept is not limited to flow configuration function initially proposed by the OpenFlow protocol. SDN in wireless networks takes advantage of the centralized view of the controller which is very useful for efficient channel assignment, power control, interference management and resource allocation. Some proposed solutions are compatible with OpenFlow while others are not OpenFlow compliant [START_REF] Haque | Wireless Software Defined Networking: A Survey and Taxonomy[END_REF]. Many research contributions are related to the control of virtualized and programmable radio networks.

SDN-based Wi-Fi Direct Clustering

Network architecture

In this paper, we propose an SDN-based architecture for Wi-Fi Direct clustering as illustrated in Figure 4.

SDN controller(s) AP1 AP2 AP3 GO1 GO2 GO3 GO4 GO5 Group1 Group2 Group3 Group4 Group5
Fig. 4 Network architecture

In the control plane, one or more SDN controllers are connected to the network of APs. When a new user wants to join the Wi-Fi network, SDN controllers and the overlay authentication service ask the user for information not only for user authentication such as user name and password but also the information related to the position of the user such as the seat number in an auditorium. Based on the number of users and their positions, the SDN controller can decide to form Wi-Fi Direct groups to optimize the network performances. For this purpose, the SDN controller can activate certain users to be GOs, configure their transmission powers and assign them the frequencies to be used. The SDN controllers also tell the other devices the Wi-Fi Direct group that they should join.

In the data plane, a device acting as a GO plays the role of an AP for the clients within its group and the role of client vis-a-vis the APs in the campus network infrastructure. In other words, it should be able to forward data between a client in the group and the infrastructure AP. From a network architecture point of view, the GO in our architecture is equivalent to an ad-hoc node which plays both roles of host and router with a routing policy controlled by the SDN controller. An OpenFlow compliant protocol can be used for the communication between the GO and the SDN controller for this purpose. For other devices which are not GOs, the data transmission is simpler. They select the specified GO in the same way that they select an AP among the available APs in the surrounding and use the GO as their single point of attachment to the network. All data from a client to the network will be sent to its GO. With this architecture, network administrators can develop advanced control algorithms for the SDN controller such as clustering, power control, chan-nel assignment or data transmission scheduling. The centralized view on network topology at the SDN controller can enable efficient clustering and interference management algorithms especially for dense wireless networks. The control in the above-presented network architecture can be summarized by a five-step control loop as shown in Figure 5. The control loop repeats periodically at a frequency configurable by network administrator or every time there is a new user coming in the room. When a student changes the seat, the seat number can be updated through the application developed by the event planner. Indoor positioning techniques can also be used to detect the position of the user. However, in the case of campus networks, requesting users to update the seat number through an application is feasible. If the moving device is a client and the user moves out of the range of its GO, it will automatically detect the presence of the new GO close to the new seat thanks to beacon messages and perform the group joining procedure. If the moving device is a GO, the clients of this group will detect its absence due to the absence of beacon messages. As a consequence, the clients will try to connect to the access point of the auditorium using the infrastructure mode, similar to the case of new arrivals. The SDN controller can intelligently detect this case and assign a new GO node for the group. The details of these control steps are presented in the following subsections.

Clustering Algorithm

It is possible for the SDN controller to use different clustering algorithms. In this paper, we use the K-means algorithm [START_REF] Cabria | Potential-K-Means for Load Balancing and Cost Minimization in Mobile Recycling Network[END_REF] which divides a set of N nodes in the network into K groups called clusters by minimizing the sum of the distances from the nodes belonging to a given group to the center of gravity of that group, called the centroid. Applied to our network architecture presented in Figure 4, the node that is closest to a centroid is designed as GO for the corresponding cluster. Each cluster corresponds to a Wi-Fi Direct group.

K -means algorithm is an unsupervised learning algorithm which is used to solve clustering problems. It is used to compute K clusters out of a given set of N elements based on their positions. The aim of this algorithm is to find K centroids, one for each cluster.

Our network N is composed of a set of N = |N | devices. For ease of notation, we map the devices in N to the integers 1, . . . , N . Device positions are stored in a set P = {p 1 , p 2 , . . . , p N }, where p j = (x j , y j ) is the position of device j.

C = {C 1 , C 2 , . . . , C K } is the partition of N into K clusters computed with the K-means algorithm. O = {o 1 , o 2 , . . . , o K } is the associated set of centroids (o k = (x o k , y o k ) denotes the coordinates of the centroid of cluster C k ).
Informally, the K-means works as follows: initially, K centroids are randomly chosen among the coordinates of the N devices. After some iterations, the K-means algorithm returns a set C of K clusters, with the associated cen-troids, which forms a Voronoi diagram: each node belongs to the cluster whose centroid is closest to it. This property is summarized by Equation (1):

∀C k , C l ∈ C, C k = C l , ∀j ∈ C k , dist(j, o k ) < dist(j, o l ),
(1) where dist is the Euclidean distance:

dist(i, j) = (x i -x j ) 2 + (y i -y j ) 2 .
Algorithm 1 K-means algorithm end while 12: end function Algorithm 1 details the procedure which computes the cluster partition C and the centroids O. The inputs of these algorithms are the description of the set N , (including the set of positions P), and the number of the clusters to calculate, K. As explained above, the first centroids are chosen with random values among the locations of the devices. Then the function iterates until we converge to a set of centroids that doesn't change. This means that the centroids of the i th and the (i + 1) th iterations are the same. At every loop, we first save the values of the centroids of O in order to compare them with the new centroids. Then, the assignClusters(O, N ) function builds a cluster C by computing, for each point j ∈ N , the minimal distance between j and every cluster centroid o k , and assigning j to its closest centroid.

Algorithm 2 details the instructions of this function.

Algorithm 2 Assign Clusters

1: function AssignClusters(O,N ) 2: C ← {C 1 = ∅, C 2 = ∅, . . . , C |O| = ∅} 3: for i ∈ N do 4: k = arg min l∈[1,K] dist(i, o l ) Find the closest centroid 5: C k = C k ∪ {i}
Add device to cluster of closest centroid 6: end for 7:

return C 8: end function Back to Algorithm 1, the cluster set computed by the function assignClusters(O, N ) will then be used next to compute new centroids in line 7. For each cluster C k , the distance between its centroid and the points is minimal in comparison with the other centroids of O. After this step, function calculateCentroids(C) computes the new centroids of the new clusters. The new centroid set O contains K centroids. The coordinates of the centroid o k of a cluster C k ∈ C are obtained by averaging the coordinates of its points (see (2)).

x o k = i∈C k x i |C k | , y o k = i∈C k y i |C k | (2)
If the centroids did not change after an iteration (O = O ), the algorithm stops iterating and returns the optimal k-partition for the network. For each cluster, the device that is located in the closest position to the corresponding centroid will be elected as GO for data dissemination.

Interference management and transmission scheduling

As indicated in Figure 5, after the Wi-Fi Direct clustering step, the SDN controller calculates the transmission powers of GO and client in other to minimize the interference in the network. Based on the positions and the role of the devices, two types of transmission power are calculated as illustrated in Figure 6. Each GO is assigned with a high transmission power P GO so that it can communicate with the access point in the auditorium. Non-GO devices are only assigned with a small transmission power P c which is sufficient to reach the GO of the group.

After the power control procedure (Figure 5), the SDN controller verifies if there may be interference between groups. Two clients at the border of two clusters can interfere each other if they use the same frequency. To avoid interference as much as possible, the SDN controller builds an interference graph. Each node in the interference graph represents a cluster. An edge connecting two nodes in the graph represents the possible interference between members of the two corresponding groups. With this interference graph, the SDN controller can apply a coloring algorithm [START_REF] Brelaz | New Methods to Color the Vertices of a Graph[END_REF] for frequency assignment.

Finally, the data transmission is scheduled into two steps : data transmission between an access point and GO, and data transmission between a GO and their clients. For this purpose, the Notification of Absence (NoA) [START_REF] Camps-Mur | Device-to-Device Communications with WiFi Direct: Overview and Experimentation[END_REF] protocol defined in Wi-Fi Direct is used by the GO to switch between the Wi-Fi Direct mode and the infrastructure mode. This protocol allows the GO to notify the clients in the group that it will be absent for a moment then switch into the infrastructure mode to get the data via the access point.

To some extent, the interference management problem in Wi-Fi Direct is similar to Femtocell networks. In Femtocell networks [START_REF] Saquib | Interference Management in OFDMA Femtocell Networks: Issues and Approaches[END_REF], a set of femto base stations are deployed in the same coverage area of macro base stations. A user can be associated to a femtocell or a macrocell. The femto base stations in cellular networks are equivalent to the GO nodes in our architecture. The macro base stations in cellular networks are equivalent to the Wi-Fi access points in campus networks. Algorithms for frequency assignment and power control designed for Femtocell networks may be applied to Wi-Fi Direct. However, the main differences between Femtocells and Wi-Fi Direct reside in the medium access control mechanism and the orthogonality of the wireless channels. The access to radio medium in Wi-Fi networks is based on random access while cellular networks rely on radio resource allocation. While the frequencies in Femtocell networks are orthogonal, the channels in Wi-Fi networks are almost partially overlapping. These differences make the interference management in Wi-Fi Direct more difficult and the algorithms designed for Femtocells less applicable to the Wi-Fi Direct environment.

By dividing the whole network into smaller groups and reducing the transmission power within the groups, our approach can reduce the overall interference in the network and increase the network performances. However, the transmission power of the GO devices is higher than the rest of the devices. Residual energy state can be included as a parameter in the GO selection algorithm run by the SDN controller to select the GO nodes in addition to the position-based clustering method. It is possible to change the GO role to another device at a reasonable frequency to balance the battery lifetime among devices in the network.

Optimized number of clusters

The K-means algorithm can only be applied when the number of clusters has been determined. It is necessary for the SDN controller to determine the number of clusters K in the network before applying the K-means algorithm.

While a complete, exact determination of an optimal value of K is beyond the scope of this paper, we propose as a rule of thumb to choose K to be of the same order as √ N . This choice is motivated by considering and solving the following simplified performance model:

-We suppose to begin with that all devices share the same frequency (the more realistic case of multiple frequencies will be deduced afterwards). -We assume that the communications can be described by a disk graph:

each device p has a given communication range ρ p . p can send information to any device in the disk of center p and radius ρ p , but will also create interferences with all other devices inside that disk. p has no impact outside its communication range. -We make the assumption that the overall performance of the system is nonincreasing with the number of interferences. In other words, to maximize the performance, we need to minimize the total number of interferences.

Each of the K GO needs to communicate with the access point. As the distance between a GO and the access point is of the same order as the dimensions of the room, we expect that a typical device will receive interferences from a non-negligible fraction of the GOs, i.e. αK for some constant K.

Then, each device needs to communicate with its GO, so we can expect that the transmission range of a typical non-GO device will have the same order of magnitude as the radius of a cluster. Assuming that the devices are uniformly distributed in the room, this means that the number of devices within the interference range of a typical non-GO device should be of the same order as the size of a typical cluster, i.e. β N K for some constant β > 0. In the end, the typical number of interferences perceived by a device should be like αK + β N K , which is minimized for K = β α N . Intuitively, α and β should be similar, as they both measure the same type of geometrical ratio, only at different scales (room and cluster respectively). This leads to K ≈ √ N . Note that the same type of reasoning should apply if one uses multiple frequencies: both α and β would be roughly divided by the number of orthogonal frequencies available1 , leaving the result unchanged. Similarly, the transmission modeling could be refined, for example by considering some AWGN channel instead of a simple disk range, but this should not change the main reason of the √ N rule. That means that the interference comes:

from one part from the communications by the K GOs to the AP, from another part from the communications by devices to their GO, which impact a number of devices proportional to the typical size of a cluster, N K . Note that in a dynamic context, the value of N may fluctuate with time, changing the optimal number of clusters. As modifying the cluster architecture may be costly, it would be best to delay the modification with some hysteresis mechanism. For example, one could change the number of clusters only when the difference between current and actual K is greater than a given threshold.

4 Performance evaluation

Clustering performance evaluation

In this section, we evaluate the data rate that a client can receive in function of the number of clusters that we have in the network. Let's consider a network of N devices divided into K clusters. If we suppose that medium access control mechanism is fair to all the device sharing the same medium, each device will have the same probability to access the channel and share the channel capacity in an equitable manner. Figure 7 shows the average data rate that a non-GO device should receive following the assumptions presented in Section 3.4. All devices share the same frequency. The capacity of the channel is shared between the N -K K devices within a cluster and the K GOs. As we can see in Figure 7, for an 802.11g network working at 54 Mbps, the optimized average data rate can be obtained when the number of clusters K ≈ √ N . For a class room of 20 students, the optimized number of clusters is 4 or 5 offering a data rate of 6.75 Mbps per device. For an auditorium of 200 students, the optimized number of clusters is 10 or 11 with the offered data rate per device of 1.86 Mbps or 1.91 Mbps respectively. It is worth noting that these results on data rate are theoretical and do not take into account the overhead introduced by protocol headers at different layers and the reduction of data rate due to the Back-off mechanism in 802.11 medium access control protocol.

For recent 802.11 standards offering higher data rates such as 802.11n or 802.11ac, the optimized data rate that a user can obtain is also higher. As illustrated in Figure 7, in the case of 802.11n network working at 150 Mbps, each user in a classroom of 20 students can receive a data rate of 18.75 Mbps and each user in an auditorium of 200 students can receive a data rate of 5.32 Mbps. Standards supported by a device can also be used as a parameter in the GO selection algorithm run by the SDN controller. A device supporting among others a more recent 802.11 standard with higher data rate should be a better candidate for a GO node.

The clustering computational time depends on the number of the devices in the network N and the number of clusters K. Figure 8 presents the measurements of computational time of the K -means algorithm in function of number of clusters to be formed. We can see that the computation time of the Kmeans algorithm linearly increases with the number of cluster K and the size of the network N. However, for a configuration of 200 devices and 200 clusters, the computation time does not exceed 0.3 seconds which is acceptable for a campus network.

Simulation with OMNeT++

In order to evaluate the data transmission performances in a more realistic conditions, we have implemented Wi-Fi Direct in OMNeT++ [START_REF] Iskounen | Wi-Fi Direct Simulation for INET in OMNeT++[END_REF]. The simulated network topology is presented in Figure 9. In this network, 50 tablets are arranged in a grid topology composed of five rows and ten columns. The horizontal and vertical distances between two adjacent devices are one meter which is approximately equivalent to the distance between two adjacent seats in an auditorium. The data rate of the Ethernet connection between the server and the access point is 100 Mbps. An SDN controller is used to control the Wi-Fi Direct clustering and data transmission in the network. The radio parameters configured in the simulation are presented in Table 1. To evaluate the performances, we have defined two test scenarios. In the first scenario, all the 50 nodes download a file directly from the access point using the infrastructure mode. In the second scenario, the 50 nodes are clus-tered into K Wi-Fi Direct groups. Transmission power of non-GO devices is limited to only reach their GO. Data are transmitted between a device within a cluster and the access point via the GO of the cluster.

Figure 10 shows the average download time in function of the file size when the 50 nodes are clustered into 4 Wi-Fi Direct groups. The file size has been varied between 50 KB, which can correspond to an excel file containing the time table of the class, and 2 MB which can be the size of a heavy PDF file of 50 slides for a lecture. These simulation results show that Wi-Fi Direct clustering can reduce up to 34% the average download time compared to the Wi-Fi infrastructure mode. To understand the the origin of performance gain in terms of download time, we consider the relationship between the packet error rate and the download time. Figure 11 shows the download time of a 50 KB file and the packet error rates in three cases of clustering, K = 1 (i.e. the infrastructure mode), K = 4 and K = 8. The packet error rates explain the gain in download time that we have obtained. When all user devices use a high transmission power to communicate with the access point to download the file simultaneously, the interference level is high leading to a high packet error rate of 21% and a high average download time of 4.68 seconds. When using Wi-Fi Direct clustering and reducing the transmission power of the non-GO devices, the overall interference in the network is also reduced, resulting in much lower packet error rates, 14.1% in the case of 4 clusters and 13.08 % in the case of 8 clusters.

To analyze the impact of the number of clusters on the overall performance, we measure the download time, the group formation time and the packet error rate in the cases of 4, 6, 8 and 10 clusters. Figure 12 shows that the configuration of 8 clusters can offer better performances in terms of download time and packet error rate than the configuration of 4 clusters, 6 clusters or 10 clusters. This observation also confirms our calculation of the optimized number of clusters should be K ≈ √ N . In this paper, we have proposed the use of SDN to control the Wi-Fi Direct clustering in campus networks. This approach is particularly useful in dense wireless networks such as auditorium and class rooms. When diving a dense wireless networks into smaller Wi-Fi Direct groups, we can reduce the transmission power of the clients in each group hence reduce the overall interference in the network. We also propose a simple method to calculate the optimized number of clusters. The simulation results have shown that our proposed architecture gets better performances in terms of download time and packet error rate compared to the Wi-Fi infrastructure mode. In the future work, we will consider data caching, multicast transmission and develop a prototype with Android tablets to be used in university campus.

Fig. 2

 2 Fig. 2 Standard Wi-Fi Direct group formation

Fig. 3

 3 Fig. 3 SDN in Wi-Fi networks

Fig. 5

 5 Fig. 5 SDN control loop

Fig. 6

 6 Fig. 6 Power control combined with Wi-Fi Direct clusetring

Fig. 8 K

 8 Fig. 8 K-Means algorithm computation time measurement

Fig. 9

 9 Fig. 9 Simulation topology

Fig. 10

 10 Fig. 10 Average download time at K = 4

Fig. 11

 11 Fig. 11 Impact of Clustering on packet error rate and download time

Table 1

 1 Simulation parameters

	Parameter	Value
	Operation mode	802.11g
	Carrier frequency	2.4 GHz
	Channel bandwidth	20 MHz
	Radio transmission rate	48 Mbps
	Receiver sensitivity	-85 dBm
	SINR threshold	4 dB
	Background noise power	-90 dBm
	GO transmission power	1.47 dBm
	Client transmission power	-23 dBm

IEEE 802.11 has 11 frequency bands in the US, 13 in the UE; however, the number of orthogonal channels is about 3 in both cases.