Nicolas Bozon
email: nicolas.bozon@gmail.com

Carole Sinfort

Bijan Mohammadi
email: bijan.mohammadi@universite-montpellier2.eu

A GIS-based atmospheric dispersion model

Keywords: Atmospheric Dispersion, risk analysis, reduced-order modelling, Gaussian plume model, multi-scaled, digital elevation models, open source GIS development, Quantum GIS MOTS-CLÉS : Dispersion Atmosphérique, Analyse de Risque, Modélisation d'ordre réduit, modéle Gaussien, multi-échelle, modèle d'élévation numérique, SIG libre, Quantum GIS

Atmospheric pollution due to the use of agricultural pesticide is a major concern today, regarding both public health, sustainable agriculture and ecosystems quality monitoring. This work aims to model pesticide atmospheric dispersion and to propose an useful air pollution prediction tool, using fluid mechanics equations and open-source GIS programming capabilities. This work focuses on pesticide spray drift modelling for viticulture applications from the plot to the watershed scales. Our goal is to provide a GIS based atmospheric transport simulation platform at very low calculation cost using solution space reduction and reduced order modelling. This is necessary to make affordable simulation and very helpful to quickly perform georeferenced pollution predictions, which can then be used as a basis for risk analysis. The coupling of fluid mechanic equations with GIS is first described, regarding the enhancement of the mathematical model according to GIS parameters like topography and scale variations. Indeed, ground variations are implemented in the model using digital elevation model (DEM). The use of DEM is an asset to replace former Cartesian outputs by georeferenced raster pesticide cloud and also, it allows to use the model on any area. Furthermore, the implementation of topography in the spatialized Gaussian plume model allows scale changes. The plateform has been developed as a plugin for Quantum GIS open source software, which offers many spatial modelling programmation capabilities. This paper describes the platform and proposes vector and raster based simulation results.

RÉSUMÉ. La pollution atmosphérique liée à l'usage agricole des produits phytosanitaires représente aujourd'hui un enjeu d'actualité pour la santé publique, pour l'agriculture durable et pour le contrôle de la qualité des écosystèmes. Un modèle de dispersion atmosphérique des pesticides a été développé à partir d'équations de la mécanique des fluides intégrées dans une programmation SIG pour fournir un outil de prediction de la pollution de l'air . Ce travail concerne la modélisation de la dérive des pulvérisations lors des applications sur vigne, de l'échelle de la parcelle à celle du bassin versant. Notre objectif est de fournir une plateforme de simulation du transport à très faible coût de calcul en utilisant des techniques de réduction de l'espace. Cette méthode permet de fournir rapidement des prédictions de pollution géoréférencées qui peuvent par la suite servir de base pour une analyse de risque. Le couplage entre équations de mécanique des fluides et le SIG est tout d'abord décrit en insistant sur les modifications apporées au modèle mathématique pour l'intégration de données géographiques comme la topographie ou les variations d'échelle. En effet les variation de niveau du sol ont été intégrées dans le modèle à partir d'un modèle d'élévation numérique (DEM). L'utilisation du DEM permet aussi de remplacer les descriptions du nuage en coordonnées cartésiennes par des données raster et permet de travailler sur n'importe quelle zone. De plus l'intégration de la topographie dans la modélisation d'un nuage Gaussien spatialisé permet de travailler à plusieurs échelles. La plateforme a été développée sous la forme d'un plugin du logiciel libre Quantum GIS (QGIS) qui offre de nombreuses possibilités de programmation avancée. L'article présente les développements réalisés et les résultats des simulations vecteur et raster.

Introduction

Let us consider environmental phenomenons presented by the atmospheric dispersion process. These are common in agricultural pesticide applications. The latter is part of our research interests, and this paper focuses on pesticide spray drift modelling over vineyards, using mathematical modelling and Geographical Information Systems (GIS). Our aim is to provide a platform dedicated to georeferenced atmospheric transport simulations at very low calculation cost, using reduced order modelling and open source GIS.

Atmospheric dispersion equations are based on parameters which are strongly affected by spatial and temporal variations [1]. These variations significantly affect the spray drift's behaviour and trajectory. The spatial dimension is thus essential in atmospheric dispersion modelling, but also represents the Geographical Information System's paradigm [2]. GIS have so become an adequate tool to analyze and visualize spatial based environmental models [3]. The coupling of the model with such technologies is presented here, regarding both the use of DEM and the role that GIS plays in the enhancement of the former model.

The optimization of the mathematical model regarding the spatial dimension first appears as necessary as we wish to couple a low-complexity model with the landscape representation provided by GIS. This work implies to establish links between the spray drift model and GIS algorithms, by selecting the most significative geospatial parameters that can improve the calculation. Topography effects and scale variations have been naturally chosen as one needed to get an optimized height dimension, in order to implement more reality in the drift phenomenon calcultation. Then, a second side of the coupling consists in integrating the model within a GIS environment in the aim to manage and automate the model inputs and to promote its outputs using cartographic rendering.

In a previous work we have treated the problem of agricultural pesticide emission and dispersion, by coupling local (i.e. emission and near-field distribution) and global (i.e. transport over large distance) models [4,6]. In this approach, the local model provides the inlet conditions for the levels above. The dependency between levels is a major asset to avoid the solution of partial differential equations using model reduction, as detailed in previous publications. [4,5,6,7] The goal of the present work is first to improve the existing long range transport model, and then to account for non uniform topographies to get more realistic simulations. We also want to be able to locally optimize the calculation accuracy according to scales variations, by using nested numerical zoom and modifying the wind flow field construction algorithm. Thus, two main aspects are presented in this paper. The optimisation of the dispersion model according to topography and scale changes is first presented using mathematical modelling. The coupling with GIS is then detailed in a second part, regarding the model integration and the main computational tasks.

Reduced-order modelling

Transport and non-symmetric geometry

We consider the situation of a source releasing a time dependent quantity c inj (t) in the atmosphere at a given location, and aim to develop a low-complexity model to represent the dispersion of this quantity. The primary factors influencing the dispersion of a neutral plume are advection by the wind and turbulent mixing. The simplest model of this process is to assume that the plume advects downwind and spreads out in the horizontal and vertical directions. Hence, the distribution of a passive scalar c, emitted from a given point and transported by a uniform plane flow filed U along x coordinate, can be represented by:

c(x, y, z) = c c (x)f ( y 2 + z 2 , δ(x)) [1] where c c (x) ∼ exp(-a(U )x) and f ( y 2 + z 2 , δ(x)) ∼ exp(-b(U, δ(x)) y 2 + z 2 )
c c is the behavior along the central axis of the distribution and δ(x) characterizes the thickness of the distribution at a given x coordinate. An analogy exists with plane or axisymmetric mixing layers and neutral plumes where δ is parabolic for a laminar jet and linear in turbulent cases [12,8]. a(.) is a positive monotonic decreasing function and b(., .) is positive, monotonic increasing in U and decreasing in δ. In a uniform atmospheric flow field, this solution can be used for the transport of c + above. This solution has been then generalized in a non symmetric metric defined by migration times calculation based on the flow field as explained in previous works [4,5,6,7] . For ground variations modelling, the flow is locally rotated to remain parallel to the ground (see figure 1) In realistic configurations, simulation needs be carried out over several hundreds square kilometers domains. At the same time, we must be able to account for local topography variations with details provided every few meters. As wind measurements are most of the time available on very coarse grids, it is unrealistic and inefficient to perform the whole simulation with a metric topographic accuracy. Rather we would like to account large scale variations of topography on a coarse level simulation and include gradually the details of the ground variations near the main points of interest. To perform this task, one recursively applies the modelling described above on a cascade of embedded rectangular homothetic domains ω i , i = 0, ..., with ω 0 = Ω the full domain1

No information is transferred at this time from fine to coarse. Indeed, we emphasize that the grids correspond to the locations where topography is available. As we mentioned, our approach is mesh free in the sense that no meshes is used for calculation. Only, evaluated information on wind and concentration are stored at these locations. The total wind field is expressed as:

U H = n level i=0 u i χ(ω i ) [2]
where u 0 = u H is calculated in equation for the coarser level and χ(ω i ) is the characteristic function for the subdomain on which level i is defined. In other words, the correction is equal to zero outside ω i . n level is the total number of levels used. For i > 1, velocity restriction from level i -1 to i is evaluated using the information from the observation point being the information at the four corners q j of a rectangle:

u i = -∇φ i , -∆φ i = 0, φ i (q j ) = φ i-1 (q j ), j = 1, .., 4
Once again we can take advantage of the linearity of the operator, to use a similar decomposition for u i than for u 0 where the observation quantities are becoming the values at the corners of the homothetic restriction:

u i (x) = 4 j=1 λ j (x) u i-1 (q j ) [3] W here 4 j=1 λ j (x) = 1, and λ j (x i ) = δ ij , i, j = 1, .., 4
If u i-1 is divergence free, then this construction guarantees that

∂ω i u i . n i dS = ∂ω i u i-1 | ∂ω i . n i dS = 0
Hence, the velocity restriction in ω i remains divergence free and is compatible with the overall field. In the simulation presented here, three levels has been used to link

Ω = ω 0 ∼ 10km 2 to the ω 2 ∼ 10m 2 .
Once the velocity restriction is defined, the concentration restriction is calculated as follow:

c i = 4 j=1 λ j (x)c i-1 (q j ) [4]
This construction guarantees that the total mass in ω i fits the entering quantity 2 :

ω i c i dV = ∂ω i c i ( u i u i . n i )dS
In practice, we use the same sampling for each level, typically a 20×20 grid which means topography needs to be defined every 500 meters at coarse level and every 50 cm at the finest level.This is interesting as it permits to locally provide a more accurate Digital Elevation Model.

Integral data

Once the species distribution c(x, y, z) is found and that the total injected quantity in time interval [0, T ] is known, one can assumes:

K = T 0 c inj (t)dt
Various quantities can be computed. For instance, we can have an estimation of the amount of species which has reached the ground using:

C g (x, y) = z≤z0 c(x, y, z)dz
or estimate the quantity still in the atmosphere beyond a distance R 0 from the source, using:

C a = K - R≥R0 C g (x, y)dV R = x 2 + y 2 corresponds
to the radius from source.

Multilevel correction for ground variations

Multi-level correction

At this point we would like to account for the topography or ground variations ((x, y) → ψ(x, y)) in the prediction model above. These are available from Digital Elevation Models (DEM) [26], [25]. Despite this plays an important role in the 2. The four points trapezoidal rule is exact for numerical integration of bilinear functions

ω f (x, y)dxdy = 1 4 4 j=1 f (qj).
dispersion process, it is obviously hopeless to launch direct simulations using a Computational Fluid Dynamic (CFD) model, based on a detailed ground description . We should mention that ground variations effects are implicitly present in observation data for wind. However, as we said, wind observations are quite incomplete. In particular, wind measurements are available every few kilometers while topographic data are available on a metric basis. At each level i of the construction (see figure 2) we introduce a correction u t i to the restriction

U H = n level i=0 ( u i + u t i )χ(ω i ) [5]
Various local modelling can be considered for u t i going from simple algebraic expres- sions to more sophisticated local CFD models. We propose the following correction 3

u t i = - 1 ρ sgn(U t )∇p i , p i = p i-1 r U 2 t [6]
where

U t = u i-1 u i-1 . n t i and p i-1 r = 1 2 ρ( u i-1 . n i ) 2 - 3
. This is the Bernoulli-Newton formula widely used in aeronautics and reproducing well the pressure distribution over a cylinder for a potential flow.

ρ is the density of the fluid. p r is a local pressure reference based on averaged entering velocity into subdomain i:

( u i-1 . n i ) -= 1 n - 4 j=0 min(0, u i-1 (q j ). n i j )
1 ≤ n -< 4 being the number of entering flow corners. The normal to the ground evaluated from the digital terrain model restriction at level i is denoted by n t i . This is different from the normal n i to subdomain i. In absence of ground variations the two normals are orthogonal (see figure 3).

Figure 3. Sketch of topography variation and normals definitions

In case the topography is not constant, we have

u i-1 . n t i-1 = 0, but u i-1 . n t i = 0
This multi-level correction improves the predictive capacity of the model introducing a dependency between ground variations and migration time. However, this is not sufficient to correctly account for ground variations in dispersion. For instance, it is indeed clear that, even in a uniform flow, cross diffusion is not symmetric on a sloppy ground when dispersion is performed parallel to the iso-level contours (see figure (see figure 4).

We also need to correct the functions a and b appearing in the dispersion modelling. As we have assumed the construction is only coarse to fine without feedback from fine to coarse levels, we assume the correction conservative in the sense that the incoming mass into subdomain i:

K i = c i-1 u i-1 u i-1 . n i - = 1 n - 4 j=0 min 0, c i-1 (q j ) u i-1 (q j ) u i-1 (q j ) . n i j
K i = ω i c i dV = ω i c i t dV
where c t is the modified expression for the concentration to account for topography changes and n i j is the normal to face j = 1, .., 4 of subdomain i. This implies a constraint on the modified expressions of a and b (e.g. correction on b can be deduced from a) through our analytical dispersion model.

K i = ω i c i t (a, b)dV [7]
The correction in a is a scaling by a positive monotonic decreasing function worthing one in absence of topography changes. For instance, one can assume:

a i t = a i u i + u t i u i
Hence, in case a change in topography increases the local velocity the dispersion goes further downstream with less cross-diffusion due to decreasing b through constraint.

GIS based coupling

Linking Dispersion model and GIS

It should first be reminded that we want to couple the above physical model which involves unsteadiness and uncertainties, with GIS that are tending to provide an accurate numerical copy of the study area's surface. Thus, GIS can be used to apply the model in a richer georeferenced numerical environment. GIS capabilities regarding DEM generation and exploitation are significally improving the former dispersion equations, as they allow the model to be run on any local topography [26]. Furthermore, GIS permits to map directly the drift process and to get standard atmopsheric concentrations at given geographical coordinates. Then it becomes rather easy to make the pesticide cloud interacts with other relevant geodata, and so to proceed to advanced risk analysis [32] regarding for example bystanders exposure, agricultural plots or water courses contamination after treatments. Although GIS allows to gain some more precision regarding topographic impact on calculations and to georeference the model outputs, we should keep in mind that the main objective of the reduced order modelling approach is to provide mean tendancies of the spray drift with very low calculation cost, and that potential errors are duplicated into the GIS. In addition to this, automated geoprocessing tasks that are carried, like DEM generation or point-based data interpolation can also add some more spatial incoherence and introduce a new level of uncertainties. As those limits regarding precision and application on real situation have been raised, it is now interesting to explain more precisly how the model and GIS comunicate, and how we can get the best of geospatial techniques to improve the model's efficiency.

Loose coupling method

Several way to couple GIS and environmental models are known in litterature, mainly "'tight"' or "'loose"' coupling and more recently "'integrating"' systems [20,21]. Each technique presents assets and limitations, and is more or less adapted depending on the complexity of the model. In our case, the loose coupling has been chosen for several reasons that have to be explained. As the latest describes an approach where interfaces are developed with minimal assumptions between the sending/receiving parties, therefore the risk that a change in one application will force a change in another application is reduced [22]. Loose coupling has also multiple assets regarding development costs, as we want to couple the model with existing GIS, and not coding a entire GIS software able to implement natively the dispersion model, as the integrated approach would suggests. As more and more GIS programs are being made available by open-source communities, we opted for Quantum GIS (QGIS) software to achieve the coupling, as it is one of the most highly capable open-source tool which offers advanced programming possibilities [27]. Indeed, QGIS is based on a robust C ++ API that presents plenty of spatial algorithms and native GIS functions. Those have been made acessible through Python bindings, which allow a simpler programming environment for developping specific QGIS plugins that directly interact with the core source code [28]. We opted for this technical solution to propose a friendly-user atmospheric pesticide spray drift plugin. It is dedicated to agricultural atmospheric pollution prediction, and has been designed to be fast and extensible, mixing reduced-order modelling and GIS development.

GIS as input data provider

The first roles of GIS deals with the automatic DEM extraction, needed by the model to compute the effects of ground variations on the windfield and so on the pesticide cloud movements. As the multi-leveled approach has been conceptualized to gain in topographic accuracy, one has to work with several DEM resolutions and must be able to extract pixel values from any loaded DEM in the GIS. Using the Python bindings, this can simply being done with some common GDAL (Geospatial Data Abstraction Library)4 commands. In our case we use two successive gdal translate commands [30], as described bellow:

gdal_translate -ot Float32 -projwin " str ( xmin ) + str ( ymax ) + str ( xmax ) + str ( ymin ) + " input_dem . tif output_dem . tif

This first command is done to clip the loaded DEM according the user defined extent in which the calculation must be launched.

gdal_translate -of AAIGrid clip . tif clip . asc

And this one to extract the elevation value of each pixel of the extent to an ESRI grid file [31]. The obtained grid is then converted into x,y,z triplets [34] needed by the model as topography inputs, using the grd2xyz python class5 Theses succesive commands permits to get the topographic input data for the dispersion model, overriding the user's DEM resolution and spatial projection as the Fortran program is then able to convert cartesian metric into the same values it reads in input.

Spatialisation of the the model

Once those input parameters are being made available for QGIS, we must use them to provide a georeferenced environment for the model's ouput data. This step deals both with some basic file formats conversion, the multi-leveled equation implementation and some advanced geodata processing. The georeferencing technique is first presented, and then we show how the GIS deals with the multi-scaled dispersion model. Some cartographic ramblings are finally presented, regarding the best way to map pesticide atmospheric drift.

Georeferencing the model's topographic input data

The mathematical model works on a cartesian metric basis, which is not readable as is by QGIS. As we want the plugin to be able to read any resolution in any geographic projection, the spatial properties of the image DEM have to be read and understood by the model. This is done by sending the resulting file of the gdal translate commands to the Fortran program, which one reads the given tabular x,y,z file by accessing the standard Comma Separated Values (CSV) format [30]. The generated DEM is sent to Fortran using simple Fortran open, do and read commands: Each triplets (i.e each line of the former raster matrix) is then understood by Fortran and provides the elevation data on which the calcultaion have to be computed, fo every point of the domain.

Introducing the multi-leveled algorithm

The multi-leveled correction for ground variations let the user choose the number of levels wanted (i.e nlevels in equation ( 10)), as well as their spatial extent (see figure 2). This permits to define the local area where the DEM resolution must be finer in order to compute ground variations more precisly. This "micro-scale" area can be defined for example just around the considered source plot or any other area that presents particular topography or significative obstacle (like local depression, small hill or other interesting rock-formings) to the spray drift. This is done in QGIS using an adaptation of the Region Tool 6 algorithm [33] applied in a recursive way:

def doneRectangle ( self ) : level = self . iface . getMapCanvas () . setMapTool ( self . saveTool ) self . updateBounds ( self . r . bb )
This first function permits to draw and save a rectangle on QGIS map view (map canvas) that defines the new extent. Then, the previous code allows to update the four corners of the extent and so to determine a new level for calculation. This way Region Tool can be used as many times as wanted, in order to set up the right number of levels for the calculation.

Cartography of pesticide clouds

The last step of the plugin development concerns the conversion of CSV outputs into standard GIS formats, but also the way one can enhance the cartographic rendering of the pesticide cloud. Using the QGIS API once again , we can first easily generates the model output results as ESRI shapefile (.shp) or any other OGR supported GIS vector format. This is done using the QGIS QgsVectorFileWriter class as presented bellow: uri = " plume . csv ? delimiter =% s & xField =% s & yField =% s " %( " ; " ," longitude " ," latitude " )

6. This code is inspired from one of Barry Stephen Rowlingson from the Mathematics and Statistics department of Lancaster University (UK).

v = QgsVectorLayer ( uri , " vectorial ␣ plume " ) Q gs V ec to r Fi l eW ri t er . writeAsShapefile (v , " vectorial -plume . shp " ')

Where plume.csv is the input CSV file that include longitude, latitude and atmospheric concentrations fields, and vectorial-plume.shp is the created point shapefile. Once this done, one can instantaneously apply some styling options to the created layer, in order to emphasize the concentrations values. This can be done using the QGIS QgsContinuousColorRenderer class, by alloting a symbol type to the geometries and a couple of minimum and maximum colors for the continuous color rendering (see figure 5). r . smin = QgsSymbol ( v . vectorType () ," 0 " ," " ," " ) r . smax = QgsSymbol ( v . vectorType () ," 1 " ," " ," " ) r . smin . setPen ( QPen ( Qt . green ,1.0) ) r . smax . setPen ( QPen ( Qt . red ,1.0) ) Thus, the resulted vectorial pesticide cloud is readable by any standard GIS program, and can be used in a simpler way for spatial analysis and atmospheric pollution prediction. Another point of interest for mapping pesticide clouds is the raster generation, as the spray drift is a diffuse phenomenon and that a surfacic representation is much more readable than points in this case. The raster creation can significally improve the cartographic meassage. In order to interpolate point-based values, one opted for the inverse distance algorithm, assumuming that the nearer a point to be interpolated is located to a point with known value, the more similar is the value of the point to be interpolated to the known value in close distance. This can be done using the gdal grid interpolation capabilities, using the GDAL virtual format (i.e VRT driver) [30] and playing on the power and smooting values: gdal_grid -a invdist : power =1.0: smoothing =50.0 -txe " + str ( xmin ) + str ( xmax ) + " -tye " + str ( ymin ) + str ( ymax ) + " -of GTiff -ot Float64 -l driftx driftx . vrt output . tif " ) .

readlines ()

Where -txe is the spatial extent in which to interpolate (i.e the user defined extent via the Region Tool class), -of is the desired output format and -ot the raster type. As the point-based values are interpolated over the whole domain, one has to apply a vectorial mask, in order to account only for points with values and so to kill the raster nodata. This can be done using the clipping functions of GDAL, using a gdal -clip command line. Finally and as for the vector ouputs, one can apply coloring schemes and transparency values, using the QGIS QgsRasterLayer optionnal arguments (see figure 6), as suggested bellow:

CONCLUSION AND PERSPECTIVES

A low-complexity model for the prediction of passive scalar dispersion in atmospheric flows has been presented and coupled with open source GIS. The solution search space has been reduced using a priori physical information and a non symmetric metric based on migration times has been used to generalize injection and plume similitude solutions in the context of variable flow fields.

The pesticide spray drift model has been applied on realistic topographic input data through coupling the inputs reading method with digital terrain models. Some of the GIS capabilities regarding spatial data storage and managment have thus been exploited in order to fit better to the vineyards landscape reality and so to tend to a true to life calculation. Furthermore, the multi-level algorithm and its correction for ground variations provides more accurracy. Current work now deals with the optimisation of the reading and the interpretion of the topographic data, but also with the integration and the use of several DEM into a single domain and in its implementation in the Python plugin. This is one more link between fluid mechanics equations and GIS algorithms to be established. In order to fully validate the topography effect, one will also have to realize simulations over longer time-series and different kind of slopes and microrelief, in order to compare the resulting values with real atmospheric concentration values. To achieve the terrain validation, a agricultural watershed will be monitored with air sampling devices that have to be positionned according to the major wind flows. Comparison between numerical results and chemical air analysis are planned for the future.

Finally, a Quantum GIS python plugin for atmospheric pollution prediction has been presented and detailed by illustrating the coupling concepts and explaining some of the functionnal code snippets. One of the major asset of the reduced-order modelling approach is to simplify the programming aspects of the coupling, and the same logic has been used regarding the open source GIS development. The result is that both dispersion model and GIS can communicate each other but stay independent. Although the plugin is already an usable and friendly-user tool for pollution prediction, some more improvments have to be developed in order to provide optimised wind flow and pesticide clouds calculation and also additional automated GIS functionnalities.

Figure 1 .

 1 Figure 1. Left: A typical Digital Elevation Model (x and y coordinates range over 2 km). Atmospheric dispersion in a uniform north wind with (middle) and without (right) the DEM

Figure 2 .

 2 Figure 2. Example of three-level construction of a flow field, genereated by the model and loaded into GIS as a vectorial layer.

Figure 4 .

 4 Figure 4. Sketch of topography variation and non symmetry in cross-definition for a constant velocity field

  r = Q g s C o n t i n u o u s C o l o r R e n d e r e r ( v . vectorType () )

Figure 5 .

 5 Figure 5. Example of vectorial pesticide cloud generated as ESRI Shapefile (.shp) with applied gratued color on the concentration field.

Figure 6 .

 6 Figure 6. Example of interpolated and masked raster pesticide coud.

For the sake of simplicity, and also because this is enough rich for environmental applications, we deliberately limit ourself to rectangular configurations.

The GDAL/OGR open source library was developed by Franck Warmerdam, independent developer from Egganville (CA) and president of OSGeo since 2006.

This code is inspired from one of David Finlayson from the United States Geological Survey, Santa Cruz (USA).


Lionel M€nard

Lucien Wald

Philippe Blanc

Thierry Ranchin

Sitting of a solar power plant: Development of Web service based on GEOSS data and guidance

Keywords: Renewable Energy, Solar Radiation, Solar Power Plant, Web Service, GEO/GEOSS, Interoperability, Standard, AIP-2 Pilot 1

Renewable energy sources such as solar and wind energy offer a large untapped potential for electricity production. The exploitation of these energies requires accurate knowledge of the resources and of their availability (in space and time) as well as accurate forecasts in the different phases of an energy system life cycle. For instance, the site selection process for development of large solar systems, such as photovoltaic on open land, require data on time-averaged values of solar irradiance from which basic economic assessments of a plant concept can be made. The paper illustrates the exploitation of Earth Observation data in this context. It describes the approach of setting-up a series of Web services that implement key features in Earth Observation data exploitation and illustrate their use through a complex application in the sitting of a solar power plant. The scenario is built on GEOSS interoperability and standard guidance.

INTRODUCTION

Renewable energy (RE) sources such as solar and wind energy offer a large untapped potential for electricity production. In order to establish appropriate instruments and strategies for the introduction of RES in the market, well-founded information on demand, resources, technologies, and applications, is essential. Figure 1 shows the necessary steps of a cascade for the successful development of investments in RE, that are themselves the seed for an efficient introduction of RE sources in market. The cascade starts with the analysis of the available resources. These are e.g. maps of wind speeds or annual solar irradiance. These results can be combined with data on available areas for the different technologies -e.g., roofs and land area for PV, land area with suitable irradiance levels or wind speed for concentrating solar power and wind power-to determine the technical feasible potential of the different technologies. The level of the available resources can then be used to calculate where these technologies become economically viable, in order to determine the economic potential. It is beyond the scope of this communication to discuss further this figure; this may be found in [START_REF] Hoyer-Klick | Getting renewable energy to work: developing a database for solar energy investment decision in the Mediterranean Region[END_REF]. What is underlined in this figure is the crucial role of the knowledge of resources and potentials. For example, the site selection process for development of large solar systems, such as photovoltaic (PV) on open land, require data on time-averaged values of solar irradiance from which basic economic assessments of a plant concept can be made. Accurate, easy accessible and affordable wind information is needed at different levels of detail before deciding the sitting of a wind park as well as during the life-cycle of a wind farm. An error of a few percent in wind resource evaluation may drastically affect the profitability of an offshore wind park. The importance of resources is supported by other articles [START_REF] Cros | Survey of the main databases providing solar radiation data at ground level[END_REF][START_REF] Dunlop | Solar Energy Resource Management for Electricity Generation from Local Level to Global Scale[END_REF][START_REF] Hammer | Solar energy assessment using remote sensing technologies[END_REF]. Resource for solar energy is hardly known in an accurate manner.

There are a number of reasons [START_REF] Cros | The Availability of Irradiation Data[END_REF][START_REF] Dunlop | Solar Energy Resource Management for Electricity Generation from Local Level to Global Scale[END_REF][START_REF] Gschwind | Converting a successful research project into a sustainable service: the case of the SoDa Web service[END_REF]. Among them are: Improved access. The access to the relevant information is poor for many reasons. It is complicated by the various types of data, various storage standards, various units, various ways of expressing time (coordinated universal time, mean solar time, true solar time, local time) and diversity in information properties stored in databases: sampling support (pixel size versus or pinpoint measurement), observational period of observation, frequency of individual observations and averaging time intervals, etc.

Improved knowledge on space and time structures. Space and time characteristics of presently available data are unsatisfactory. Interpolation/extrapolation techniques can be employed to gain knowledge at any geographical location and time, but the present techniques approaches lead to poor quality estimates.

Improved matching to actual users needs. Raw measurements are stored in and supplied by the current databases. Most often, they are not what are is really needed by the user. More advanced information should be supplied to the user. The GEOSS initiative reveals its usefulness in bringing solutions to these three main concerns expressed by users as it will be detailed in subsequent sections. Standards for exchange and supply of data, interoperability are a means to improve access to data as well as extensive use of Internet and geographicallydistributed services that can be exploited by a collaborative information system [START_REF] Gschwind | Converting a successful research project into a sustainable service: the case of the SoDa Web service[END_REF]. Images from meteorological satellites together with appropriate processing methods are a means to improve knowledge on space and time structures as they offer every hour a synoptic view of the clouds and more generally of the optical state of the atmosphere [START_REF] Hammer | Solar energy assessment using remote sensing technologies[END_REF][START_REF] Perez | A new operational model for satellitederived irradiances: description and validation[END_REF][START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation data from satellite images[END_REF]. Finally, availability of data of various types from various sources is enabled by the GEOSS, especially interoperability capabilities, which improve the matching of service supply to actual users needs.

GEOSS APPROACH

GEO, the Group on Earth Observation

The Group on Earth Observations (GEO) is an intergovernmental organization with 76 members plus the European Commission and 56 participating organization. GEO is coordinating efforts to build a Global Earth Observation System of Systems called GEOSS.

GEOSS, the Global Earth Observation System of Systems

The Global Earth Observation System of Systems (GEOSS) aims at promoting a global synergy of Earth observation. It intends to provide a more easy and open access to Earth observation data, to coordinate and sustain observation systems and to foster the use of such systems in all applicable domains. Nine areas of critical importance to people and society called Societal Benefit Area (SBA) are being addressed: disasters, health, energy, climate, water, weather, ecosystems, agriculture and biodiversity as shown in Figure 2.

Figure 2: The Global Earth Observation System of Systems.

GCI, the GEOSS Common Infrastructure

The GEOSS Common Infrastructure (GCI) comprises four main elements (Fig. 3):

1. the GEO Portal allows the user to directly search for information and services available in GEOSS; 2. the GEOSS Clearinghouse is the element that collects search and presents the various existing GEOSS components to the users via the GEO Portal; 3. the GEOSS Component and Service Registry allows GEOSS organizations to contribute components and services to the community; 4. the Standards Registry enables contributors to GEOSS to configure their own systems to be compatible and interoperable with others systems. A Standards and Interoperability Forum (SIF) has been created to allow contributors to share experiences, ideas and initiate new standards proposals.

These four elements are intended to provide the users of Earth observation a mean to search, discover and access information, tools and services available in GEOSS.

Figure 3: The GEOSS Common Infrastructure (GCI)

The GEOSS Interoperability Process Pilot Projects

GEOSS aims at integrating existing systems keeping in mind synergies among GEOSS components for the benefit of the nine SBAs. Interoperability is the key concept of GEOSS. The GEOSS Architecture & Data Committee (ADC) is in charge of developing and implementing this concept in several phases. The GEOSS Architecture Implementation Pilot task aims to incorporate contributed components consistent with the GEOSS Architecture, using a GEOSS Web Portal and a GEOSS Clearinghouse search facility, in order to access services (Fig. 3). Within the current Architecture Implementation Pilot Phase 2 (AIP-2), initiatives, called Process Pilot Projects, have been issued not only to address conceptual design but to demonstrate interoperability through real and persistent use cases. In that respect, a call for participation was released in June 2008 in close collaboration between the ADC and the GEOSS User Interface Committee (UIC).

Renewable Energy Scenario: an answer to the AIP-2 Call

Mines ParisTech and its associates: DLR (Deutsches Zentrum f'r Luft-und Raumfahrt), Joint Research Center -European Commission, Meteotest, NASA (Langley Research Center), Rutherford Appleton Laboratory, answered theAIP-2 call. Thereby, they express the global concern of a large community of users covering a broad spectrum: researchers, investors, producers, utilities, consulting companies, in the field of renewable energy looking mainly for solar radiation data but having also needs to access other various data sets to fulfill their objectives. The needs for Earth observation data in RE has been addressed by the Community of Practice "Energy" (GEOSS-ECP). A scenario divided in two main consecutive phases has been proposed. It intends to federate the RE community towards GEOSS interoperability concepts in order to provide as much as possible standard access to energy related catalogues and resources. Then it plans to demonstrate interoperability concepts through a real use case for "Sitting Solar Power Plants".

FEDERATE THE RE COMMUNITY TOWARDS GEOSS INTEROPERABILITY CONCEPTS

Identify, check and validate, candidate resources for the RE SBA & Scenario

The first action undertaken has been to provide the most large and thorough list of already existing available resources addressing concerns of the RE. These resources must implement interoperability arrangements in order to comply with the AIP-2. This living list is available at the dedicated section for the RE scenario within the collaborative development site for the GEOSS AIP-2. Emphasis has been put on identifying the interoperable arrangement proposed by each resource. Various interoperable standards stand together (OGC/WMS, W3C/WSDL, JSR6168,…) in this list; each of them addresses specific concerns of the RE community. An excerpt of this list is given in Table 1.

Provider

Resource Desc. 1. However, several of them necessary for the scenario were not practically implemented as such. For the purpose of the scenario and also as a strategic move to support interoperability for the sake of the RE community, Mines ParisTech has developed a community portal (www.webservice-energy.org) for deploying, hosting, promoting and allowing access to W3C Web services useful in this community. It acts a repository where providers may deploy their services that can be launched in an automatic manner by interested, and possibly authorized, users. Mines ParisTech has deployed several of such services, providing a standard access to several well-used data sets in solar radiation. For example, the services for databases HelioClim3 and SOLEMI (Table 1) are hosted by this portal. This dedicated platform comprises in the back-end a Red-Hat JBoss application server for deploying and hosting W3C Web services. In the front-end an Apache Web server has been installed for human description of the Web services, to host the INSPIRE metadata of these Web Services and to act as a catalogue component thanks to the WAF (see later). The WAF URL has been provided to the GEOSS Registry as the end point for harvesting purpose.

Interop. Arrangement

Metadata and catalogue for resources

Implementing interoperability arrangement implies the provision of key components to ease mechanisms for search and discovery of resources. These two key components are metadata and catalogue. Various standards exists but selecting the one that best suits specific need can sometimes be difficult.

"SITTING A SOLAR POWER PLANT" SCENARIO USE CASE IMPLEMENTATION

Scenario summary

The purpose of the second part of the response of Mines ParisTech and associates is to demonstrate interoperability concepts through a real use case for "Sitting Solar Power Plants". The scenario summary comes from a real business case encountered by various solar energy stakeholders. It is presented as the following steps:  investors and electricity producers willing to invest in solar plants need precise and thorough information to support decision-making;  on their behalf, consulting companies perform feasibility studies in order to decide where to sit power plants and which technology to use ensuring a profitable return on investment;  to reach that goal, consultants need an easy and unified access to data sets. Such data sets include meteorological, geographical and environmental parameters.

Defining inputs and outputs

An overview of the inputs and outputs implementing this scenario can be summarized as follows:

As inputs it is mandatory to enter the AOI (Area Of Interest), the time period and the output scale. This will trigger the request for "core" layers. "Optional" layers can be selected on purpose. The output is a downloadable ZIP archive containing "core" and "optional" layers.

Figure 4: "Core" layers Figure 5: "Optional" layers

Resources description

The scenario addresses two kinds of resources. 1. The "core" layers systematically provided in the output archive.

The "core" layers are built and added to the archive by internal processing of Mines ParisTech resources. These resources are made of meteorological data such as map of annual irradiation (kWh m - †) extracted from the HelioClim database. Other "core" layers include shadows (%), terrain elevation (m) and local max slope (in degree) derived by an appropriate processing of terrain data from the SRTM database. 2. "Optional" layers which are optional but valuable for the renewable energy planner may be provided upon user selection. These layers will be triggered and added to the archive from already known and identified interoperable resources of GEOSS partners. A first list of such layers has been identified. It includes geographical data layers providing information on hydrological features (rivers, lake, channels), gazetteer and land use, environmental data like protected area, risks and hazards. As previously discussed in this document, a work to identify operational and persistent resources to ease the scenario rolls out has been conducted (Table 1). These resources include the Web services deployed on the RE community portal and the CEISIN Columbia University resources on hazard, risks, population density resources, implementing OGC/WMS-compliant Web map services.

Web services chaining

In order to provide users with the map layers archives, specific software developments need to be carried out. Beside having all the resources available in an GEOSS compatible interoperable format, an orchestration of request, access, retrieve actions need to be performed. Some platforms such as SSE from ESA provide a workflow framework for Web service orchestration. The use of this tool has not been considered given the time frame of the scenario with respect to the efforts that it would have requested. Our approach uses a simple daisy chain mechanism by sequentially request the optional selected resources beside the always included core ones. As stated previously, it is admitted that the resources are already known and are permanently and persistently operating. Classical time-out access parameters will be use in the programming approach to encompass possible data failure. As this approach might sound simplistic, it nevertheless illustrates data sets interoperable approach; it will underline the benefit of bringing resources alive complying to GEOSS standards and at least it will provide a valuable outcomes to renewable energy practitioners.

Client development

Bringing the above mentioned concept alive from the user perspective implies a design of a graphical user interface (GUI). By taking advantage of previous experience gained in an on-going project called MESoR and funded by the European Commission, we are working towards the provision of the GUI in a form of a JSR-168 compatible Portlet. According to JBoss terminology a Portlet is a Java technology based web component, managed by a Portlet container, that processes requests and generates dynamic content. This Portlet is deployed on a Liferay portal on the top of a JBoss platform. It is worth noticing that the Liferay portal that we are using within this scenario has been chosen by two of the current portal candidates in GEOSS. The GUI elements within this Portlet will include Google Map API for Area Of Iterest (AOI) selection and a set of HTML forms for the other parameters. It will provide the coverage of the data resources as a KML file. Similarly to what has been developed within the MESoR project we want to provide the users with meta information about the application. Such meta information can include Web service overall description, reference to IPR (Intellectual Property Rights) and credits, and description of inputs and outputs. To provide a more interactive GUI, the meta information is embedded in tabs developed with AJAX technology for asynchronous browsing and querying approach. Figure 6 provides an example of a Portlet, by which the above mentioned GUI elements trigger the HelioClim database through the corresponding Web service deployed at the RE community portal. This example fully addresses the GEOSS interoperable approach promoted in the renewable energy scenario and will be used as a template to implement the client.

CONCLUSIONS

Though the work has not yet been fully completed, it illustrates the exploitation of Earth observation data for the benefit of the development and integration of renewable energy in energy production. Within the framework of the AIP-2 initiative, a series of services has been developed that are compliant to the standards and interoperability requirements. As the RE community has more than 10-years experience in Web services, these services were developed using W3C standards such as WSDL and SOAP. The present study demonstrates the capability of the ADC to develop interoperability arrangements that accommodates various standards: W3C, ISO, INSPIRE and OGC. One of the benefits of this collaboration between the architects (ADC) and the users (UIC) is the development of the Web Accessible Folder (WAF), which is a new type of catalog.

As practitioners in RE, we have been able to achieve a more easy and open access to Earth observation data and services and to propose user-oriented services specific to RE that are compliant to GEOSS requirements.

Figure 1 :

 1 Figure 1: The cascade of market introduction of renewables (in Hoyer-Klick et al. 2009)

Figure 6 :

 6 Figure 6: HelioClim3 Portlet at MESoR prototype

Table 1 :

 1 Resources listing details excerpt

				Spatial and Time
				Coverage
	Mines	HelioClim3:	W3C WSDL Europe -Africa.
	ParisTech	Database of		Time series -Minute
		Solar		to Month
		Radiation
	DLR	SOLEMI:	JSR-168	Europe, Africa,
		Database of	Compatible	Western Asia -Time
		Solar	Portlet	series of hourly
		Radiation		values
	CIESIN	Hazards,	OGC WMS	World
		Population	and WFS
		Density,	Compliant
		Natural	Mapping
		Disaster,	Services
		Human
		Footprint,...
		layers
	Mines	RE	ISO-19119	N/A
	ParisTech	Community	Metadata
		portal	WAF

3.2 Renewable energy community portalCommunity portal

  as a whole is a recognized component of the architecture of the GEOSS Common Infrastructure. The scenario intends to provide access to interoperable resources as shown in Table

  As the RE community has developed already several W3C Web services, this implies the ability of the GEOSS to build metadata and catalogue upon W3C WSDL Web services. Though UDDI catalogue exists, its implementation requires a considerable infrastructure to operate and is not worth the effort when only tens of services are considered. Existing UDDI hosting facilities have not revealed appropriate to our case. We have therefore considered a new approach that is by far simplest but is perfectly compliant with GEOSS interoperable arrangements. This approach exploits INSPIRE ISO-19119 metadata for resources description and Web Accessible Folder (WAF) for cataloguing purpose. The INSPIRE metadata have been generated thanks to the INSPIRE metadata editor available on-line at the INSPIRE web site. This easy-to-use editor allows to generate ISO 19119 compliant metadata that describe W3C Web services through simple forms displayed in various tabs. Once the metadata are ready, they can be saved on the computer of the user. In order to provide a repository for those metadata in the community portal, a new component was designed: the Web Accessible Folder (WAF). The WAF is a basic directory where all the metadata XML files are stored. It is registered as a GEOSS standard component within the GEOSS component registry in order to offer search and discovery capabilities. This registry in then harvested by the GEOSS crawlers of the GEOSS Clearinghouse to allow the search and discovery actions of Web services in the GEOSS portal. Valid requests for searching and discovering W3C Web Service parameters have been successfully achieved and demonstrate the usefulness of such an approach.


A Oumbe

Ph Blanc

T Ranchin

M Schroedter-Homscheidt

L Wald

A new method for estimating solar energy resource

Keywords: solar irradiance, atmospheric optics, radiative transfer, cloud, Envisat, MetOp, MSG

This paper presents the first prototype of a new method -called Heliosat-4-for assessing solar surface irradiance. The Heliosat-4 method is based on the radiative transfer model libRadtran and benefits from advanced products derived from recent Earth Observation missions. It will provide direct and diffuse components and spectral distribution of solar surface irradiance every 3 km and • h over Europe and Africa. The outcomes of this first prototype of the method are compared to ground measurements made at 8 stations in Europe and Northern Africa. The results show that RMSE attained by the Heliosat-4 method is fairly similar to that attained by current methods. A large overestimation is observed. Nevertheless these tests are very encouraging.

aerosol optical depth at 550 nm (unitless)  Angstrom coefficient (unitless). It characterizes spectral variation of aerosol optical depth vis visibility at ground level (km)  s solar zenith angle (€)

INTRODUCTION

World energy consumption is sharply rising (Martin-Amouroux, 2004) when we are witnessing a gradual depletion of fossil resources. Other sources of energy should be better exploited to meet demand. Among them, solar energy offers a great potential. Solar energy resources have been for long estimated from measurements made within meteorological stations networks (WMO, 1981). However, such networks do not offer a dense nor sufficient coverage to provide an accurate description of spatial variations of the resource. As an alternative, several methods have been developed in the past years to assess the solar surface irradiance (SSI) from images taken by satellites [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF]. Most current methods are inverse, i.e. the inputs are satellite images whose digital counts result from the ensemble of interactions of radiation with the atmosphere and the ground during the downward and upward paths of the radiation. The authors believe from experience that limits of such methods are presently reached in terms of accuracy, except for those are specialized to a specific region using empirical fittings [START_REF] Polo | Iterative filtering of ground data for qualifying statistical models for solar irradiance estimation from satellite data[END_REF]. Consequently, a new paradigm is studied based on a direct modeling. This would permit to deliver knowledge on direct, diffuse components and spectral distribution, which is seldom offered by current methods.

Recent satellite missions such as Meteosat Second Generation, Envisat and MetOp combined with recent data assimilation techniques into atmospheric modeling offer a favorable context for the design and exploitation of a method based on a direct modeling. In this approach, the various optical processes occurring along the path of the light from the outer space towards the ground are modeled by the means of a radiative transfer model (RTM).

In this work, we propose a prototype of a direct method -called Heliosat-4-for assessing the SSI; it exploits advanced products derived from recent satellite data and is based on a radiative transfer model. Our final goal is to integrate it into an operational system for providing operationally direct and diffuse components and spectral distribution of the SSI for the usable Meteosat image pixels (9 million pixels, 3 km resolution at nadir of the satellite and each quarter-hour). This SSI data will be disseminated by the means of the SoDa Service (www.soda-is.com, [START_REF] Gschwind | Converting a successful research project into a sustainable service[END_REF]. This communication deals with the SSI integrated over the spectral interval [0.3 'm, 4 'm].

METHOD

Background

The Heliosat-4 method is under joint development by the DLR (German Aerospace Center) and MINES ParisTech. The libRadtran RTM has been selected for this development. libRadtran was originally developed for modelling UV irradiance. Its accuracy has also been demonstrated for actinic fluxes and total irradiance (www.libradtran.org, [START_REF] Mueller | Rethinking satellite based solar irradiance modelling -The SOLIS clear sky module[END_REF].

There are many inputs to RTMs. [START_REF] Oumbe | Exploitation of radiative transfer model for assessing solar radiation: the relative importance of atmospheric constituents[END_REF] performed an inventory of the variables (e.g., cloud) and their attributes (e.g., optical depth) available in an operational mode and assessed to which degree the uncertainty on an attribute of a variableincluding the absence of value-leads to a variation on the SSI. They found a number of significant inputs:  solar zenith angle ( s ) and number of the day in the year,  cloud optical depth (τ c ),  cloud type,  water vapor amount,  aerosol optical depth and its spectral variation (τ aer550nm , ),  aerosol type,  ground albedo () and its spectral variation,  atmospheric profile,  ground altitude (z).

Another result of this sensitivity analysis is that the influences of vertical position and the geometrical thickness of clouds in the atmosphere are negligible. Thus, the SSI I for a cloudy atmosphere can be considered as equal to the product of the irradiance obtained under a clear sky (I clearsky ) and a function of the cloud extinction and ground albedo contribution (T cloud+albedo ):

I  I clearsky T cloud+albedo (1)

Scheme

For our objectives, look-up tables can be used [START_REF] Gimeno-Ferrer | Visiting Scientist Scientist Reports, Eumetsat "RTM calculations for surface incoming shortwave radiation[END_REF]. The look-up tables are cumbersome to compute and then to implement in the routine operations. If new descriptions of atmospheric properties are available, the look-up tables should be recomputed for the new conditions and reimplemented in the processing software. They have the advantage of making the subsequent operations much faster. Here, we have opted for another approach by using libRadtran and a combination of simpler approached models that run fast enough for our purposes.

Considering Eq. 1 and taking into account that i) the variations with the altitude z of the term T cloud+albedo are negligible and ii) the cloud parameters may be known in an operational mode every ƒ h and 3 km and the clear-sky parameters only every day and 50 km, the concept of the Heliosat-4 method is that the SSI I can be computed by the products of three models as follows:

I = I clearsky f(z) T cloud+albedo
(2) daily and 50 km ƒ h and 3 km where each model differs from the others with respect to space and time scales. The first model computes the clear-sky irradiance I clearsky ; it takes into account all atmospheric parameters but clouds and ground albedo and focuses the bulk of computation resources.

The second model f(z) corrects for ground elevation for a given site and the third model take into account the clouds extinction and the contribution of the ground albedo to the SSI (T cloud+albedo ). The second and third models are much simpler and run faster than the first one. A sketch of the scheme is given in Fig. 1.

The clear-sky part of the scheme

The Heliosat-4 method exploits a clear-sky model which is in this case the RTM libRadtran combined with the modified Lambert-Beer (MLB) relation proposed by [START_REF] Mueller | Rethinking satellite based solar irradiance modelling -The SOLIS clear sky module[END_REF]. The MLB relation is a fitting function modelling the influence of the solar zenith angle in order to save computing time:

I clearsky ( s ) = I 0 exp(-τ 0 / cos a  s ) (3) I B clearsky (s) = I 0 exp(-τ 0B / cos b  s )
where τ 0 , τ 0B and a, b are determined from the knowledge of the global and direct irradiances at solar zenith angles equal to 0€ and 60€. This knowledge is obtained by two runs of libRadtran.

Inputs to libRadtran are the optical depth and type of the aerosols, the water vapour and ozone content, which are themselves assessed from data provided by the satellites Envisat and MetOp [START_REF] Holzer-Popp | Retrieving aerosol optical depth and type in the boundary layer over land and ocean from simultaneous GOME spectrometer and ATSR-2 radiometer measurements 2. Case study application and validation[END_REF][START_REF] Schroedter-Homscheidt | Total water vapor column retrieval from MSG-SEVIRI split window measurements exploiting the daily cycle of land surface temperatures[END_REF]). These inputs are provided for cells of 50 km in size.

Double-z fitting function for vertical profile of irradiance

The mean terrain altitude for each cell, z 0 , is also an input to libRadtran. The altitude of the terrain has a large influence on the irradiance and it is necessary to take into account the actual altitude of the site of interest in the assessment of the clear-sky irradiance and consequently of the SSI. It would be too timeconsuming to run libRadtran for each altitude. Given the fact that libRadtran can run with two altitudes (z 0 and z H ) as inputs in one call, we have designed the use of the double-z fitting function described below to perform the altitude correction. Knowing values of clear-sky irradiances at two different altitudes I clearsky (z H ) at z H and I clearsky (z 0 ) at z 0 , this fitting function allows to determine clear-sky irradiances at other altitudes. It is called "double-z fitting function" because it needs as inputs the clear-sky irradiances at two altitudes. The form of this function has been defined from the typical vertical profile of direct irradiance but was then found suitable to the global irradiance, i.e., the sum of the direct and diffuse components. Two attenuation coefficients, A(z), for global and, A B (z), for direct irradiances are defined:

I clearsky (z) = I 0 (1-A(z)) (4) I B clearsky (z) = I 0 (1-A B (z))
We assume that A(z) and A B (z) have the following forms:

A(z) = A(z 0 ) exp[-(z -z 0 )] (5) A B (z) = A B (z 0 ) exp[- B (z -z 0 )]
A(z 0 ) and  can then be computed:

A(z 0 ) = 1 -(I clearsky (z 0 ) / I 0 ) (6)  = -ln[(I 0 -I clearsky (z H ))/(I 0 -I clearsky (z 0 ))] / (z H -z 0 )
z H has been empirically defined to obtain the best results for elevations between 0 and 8 km and is set to, in km:

z H = max(3, z 0 +2) (7)
The same equations hold for A B (z 0 ) and  B , where I B clearsky (z H ) and I B clearsky (z 0 ) are the direct clear-sky irradiances at z H and z 0 :

A B (z 0 ) = 1 -(I B clearsky (z 0 ) / I 0 ) (8)  B = -ln[(I 0 -I B clearsky (z H ))/(I 0 -I B clearsky (z 0 ))] / (z H -z 0 )
The outcomes of the double-z fitting function are compared for several altitudes to the outcomes of libRadtran and 6S [START_REF] Vermote | Second simulation of the satellite signal in the solar spectrum (6S), 6S: An overview[END_REF] which serve here as references. libRadtran was run for several input conditions, e.g., different solar zenith angles (0,10,20,40,50,60,65,70,75,80) €, four different aerosol types, different visibilities (15,20,25,30,40,50,60,80,100) km, different amounts of precipitable water (10,15,30,60,80) kg m -2 , different standard atmospheres and different ground albedos (0.0, 0.2, 0.3, 0.6, 0.9). Similar runs were made with 6S. Similar conclusions and values were obtained, ensuring that the comparison is not biased by the use of a single RTM. Deviations between the results of the RTM and the double-z fitting function were computed. The relative value of the deviation is less than 1 % (global) or 2 % (direct) for current visibilities, i.e. greater than 30 km (WMO 1981) (Fig. 2). It may amount to respectively 3 % and 5 % for very low visibilities of 15 km. For larger visibilities, the relative deviation is almost zero for altitudes from sea level up to 6-7 km. There is little influence of the solar zenith angle on the deviation. We conclude that this function is accurate enough to be used in the Heliosat-4 method.

Clouds extinction and ground albedo contribution

Since the geometrical thickness of clouds has a negligible effect on the SSI, it is set up to a fixed value of 1 km. Then, we assume that the irradiance passing through the clouds is divided into downward and upward irradiances. Therefore, we use the wellknown two-stream and delta-eddington approximations [START_REF] Paris | A cloudy-sky radiative transfer model suitable for calibration of satellite sensors[END_REF] to model the extinction due to clouds:

I B = I B clearsky exp(-τ c / cos s ) (9) I D = t D (I D clearsky + R c  sph ) + (t B -t D ) I B clearsky
where R c is the irradiance reflected by the cloud and is given by:

R c = (r D I D clearsky + r B I B clearsky ) / (1 -r D  sph ) (10)
The effective diffuse and the direct-to-diffuse transmissivities and reflectivities -respectively r D, t D and r B, t B -are determined from the radiation boundary conditions [START_REF] Paris | A cloudy-sky radiative transfer model suitable for calibration of satellite sensors[END_REF]. The spherical albedo  sph is computed previously in the clear-sky part of the scheme.

If we consider that an infinite series of reflection and scattering takes place between the ground and the atmosphere, the SSI computed for a null ground albedo can be corrected for the actual ground albedo :

I() = I(=0) / (1 -  sph ') (11)
where  sph ' is the spherical albedo for the cloudy atmosphere. The input to the cloudy part of the scheme is the cloud optical thickness derived from an appropriate processing of the Meteosat images made at DLR. Presently, the ground albedo is set arbitrarily; a database should be found for operation.

COMPARISON BETWEEN RESULTS AND STATION MEASUREMENTS

We have applied the Heliosat-4 method to data provided by the DLR for the year 2004. Values of irradiance for every quarterhour were averaged to yield hourly means of irradiance. These assessed SSI are compared to measurements performed by pyranometers at 8 meteorological stations (Table 1). We only use measurements with irradiance greater than 50 W m -2 . For lower values, the measured irradiance is mainly of diffuse nature and is influenced by local conditions, including orography and the presence of nearby obstacles. By removing these values, we ensure better conditions for understanding the results at this preliminary stage. For the same reason, we avoid measurements taken at solar zenith angle greater than 70€. Above this angle, atmospheric parameters are not well estimated. The discrepancies between the SSI measured and assessed are computed. Table 2 displays the mean value of the SSI measured at ground level, the bias, i.e., the mean difference, the root mean square error (RMSE) and the correlation coefficient between the two series of values. The bias is positive for all sites. This means that the method over-estimates the SSI in general. The bias may be small as in the case of Tamanrasset (station 7) or large as in the case of Vaulx-en-Velin (station 8). The correlation coefficient is large in all cases: the changes in radiation are well reproduced by the Heliosat-4 method. The root mean square error ranges from 70 W m -2 to 133 W m -2 ; in relative values, it is close to 20 % in all cases, except for Tamanrasset for which it is 10 %. The solar zenith angle has an important impact on quality of the estimation: the greater the solar zenith angle, the greater the relative bias and RMSE. A similar comparison is under way for direct irradiance.

CONCLUSION

Compared to current methods [START_REF] Rigollier | The method Heliosat-2 for deriving shortwave solar radiation from satellite images[END_REF], the Heliosat-4 method should not only produce total global irradiance, but also direct and diffuse components, and spectral distribution. First results show that the RMSE attained by the Heliosat-4 method is fairly similar to that attained by current methods. These tests are very encouraging. A strong drawback is the overestimation. Several causes may be invoked: the ground albedo was set to an arbitrary value that may not be the correct one, the spectral variation of the aerosol optical depth is not taken into account as we used only τ aer550nm , as well as the aerosol type. It is also believed that a better modeling of the cloud extinction can be done if cloud effective radii are known with subsequent more accurate assessments.

This first prototype of the Heliosat-4 method shows promises. It demonstrates that principles are sound. There is space for improvements and the authors are confident that the Heliosat-4 method will reach a satisfactory quality.

  irradiance on horizontal plane under clear sky (W m -2 ) T cloud+albedo cloud extinction and ground albedo contribution (unitless) A(z), A B (z) attenuation coefficients of the atmosphere for global, direct irradiance (unitless) ,  B parameters of the double-z fitting function (m -1 )  ground albedo (unitless)  sph ,  sph ' spherical albedo for clear and cloudy sky (unitless) a, b empirical parameters for modified Lambert-Beer relation (unitless) τ 0, τ 0B optical depths for modified Lambert-Beer relation (unitless) τ c cloud optical depth (unitless) τ aer550nm

  Figure 1. Heliosat-4 scheme

Figure 2 .

 2 Figure 2. Relative value of the deviation between the RTM libRadtran and the double-z fitting function. Main inputs to the RTM are solar zenith angle (30 €), precipitable water (15 kg m -2 ), number of the day in the year (170), and ground albedo (0.2).

Table 1 .

 1 Stations used for validation

	Site	Station	Country	Latitude / Longitude(°)	Altitude (m)
	1	Carpentras	France	44.05 / 5.05	110
	2	Dresden	Germany	51.13 / 13.78	214
	3	Freiburg	Germany	48.04 / 7.87	241
	4	Geneva	Switzerland 46.25 / 6.13	407
	5	Ispra	Italia	45.82 / 8.60	191
	6	Payerne	Switzerland 46.81 / 6.94	492

Table 2 .

 2 Performance of Heliosat-4 for global SSI. RMSE means root mean square error, Num is the number of observations and Corr means correlation coefficient.

	Site Num	Mean (W m -2 )	Bias (W m -2 )	RMSE (W m -2 )	Corr
	1	2627	566.5	37.7 (7 %)	87.9 (15 %) 0.948
	2	1916	451.0	58.2 (13 %)	111.1 (25 %) 0.901
	3	936	528.7	81.5 (15 %)	121.0 (23 %) 0.932
	4	2427	491.8	23.0 (5 %)	106.5 (22 %) 0.917
	5	2405	499.9	59.6 (12 %)	101.4 (20 %) 0.942
	6	1985	556.0	61.8 (11 %)	110.0 (20 %) 0.918
	7	2747	714.9	15.5 (2 %)	69.4 (10 %) 0.966
	8	1605	542.4	87.1 (16 %)	132.5 (24 %) 0.917

ACKNOWLEDGEMENTS

This work was partly done with a support of the ADEME, the French agency for energy.


S Baron

J P Cachemiche

F Marin
email: marin@cppm.in2p3.fr

P Moreira
email: paulo.moreira@cern.ch

C Soos

Implementing the GBT data transmission protocol in FPGAs

 is a radiation tolerant ASIC that can be used to implement bidirectional multipurpose 4.8 Gb/s optical links for high-energy physics experiments. It will be proposed to the LHC experiments for combined transmission of physics data, trigger, timing, fast and slow control and monitoring. Although radiation hardness is required on detectors, it is not necessary for the electronics located in the counting rooms, where the GBT functionality can be realized using Commercial Off-The-Shelf (COTS) components. This paper describes efficient physical implementation of the GBT protocol achieved for FPGA devices on Altera and Xilinx devices with source codes developed in Verilog and VHDL. The current platforms are based on Altera StratixIIGX and Xilinx Virtex5.

We will start by describing the GBT protocol implementation in detail. We will then focus on practical solutions to make Stratix and Virtex transceivers match the custom encoding scheme chosen for the GBT.

Results will be presented on single channel occupancy, resource optimization when using several channels in a chip and bit error rate measurements, with the only aim to demonstrate the ability of both Altera and Xilinx FPGAs to host such a protocol with excellent performances. Finally, information will be given on how to use the available source code and how to integrate GBT functionality into custom FPGA applications.

I. GBT PROTOCOL PRESENTATION

A. Introduction

The general architecture of a high-speed optical link implemented using the GBT chipset and FPGA is represented in Figure 1.

Figure 1: GBT optical link implementation scheme Logically the link provides three "distinct" data paths for: Timing and Trigger, Data Acquisition and the Slow Control. In practice, the three logical paths do not need to be physically different and are merged. The aim of such architecture is to allow a single link to be used simultaneously for data readout, timing and trigger distribution, readout and experiment control. The link establishes a point-to-point optical bidirectional connection (using two optical fibers).

The GBT chipset [START_REF] Moreira | GBTx specifications[END_REF] is under development to match such architecture. It targets high-speed (3.36Gb/s) data transmission between the detectors and the counting room.

As illustrated in Figure 1, such a link is implemented by a combination of custom and Commercial Off-The-Shelf (COTS) components. In the counting room, the receivers and transmitters will be implemented using COTS components and FPGAs while, embedded on the detectors, the receivers and transmitters will be implemented by the GBT chipset and Versatile Link Components [START_REF] Vasey | Versatile Link[END_REF]. This architecture clearly distinguishes between the counting room and front-end electronics specificities: that is, the on-detector front-end electronics works in a hostile radiation environment requiring custom made components while the counting room electronics operates in a radiation free environment allowing the use of COTS components. Moreover, the availability of FPGAs with up to 48 Hard-IP serializer blocks would allow concentrating data from several front-end sources into a single module in the counting room facilitating data merging and leading to compact systems.

The study presented below will focus on proving the usability of COTS components and FPGAs to implement the GBT protocol in counting rooms [START_REF]GBT-FPGA project web site[END_REF].

B. GBT Protocol

Due to the beam luminosity planned for SLHC, the high speed data transmission link will be exposed to high Single Event Upset rates. SEUs are a major impairment to error free data transmission. To deal with this, the GBT line coding adopts a robust error correction scheme that will allow correction of bursts of errors caused by SEUs. A significant fraction of the channel bandwidth must therefore be assigned to the transmission of a Forward Error Correction (FEC) code.

The code to be used must provide a high level of protection, since errors occurring during transmission can also occur as burst errors and not only as isolated events. Because of this, a double interleaved Reed-Solomon correcting code was chosen. The code is built by first scrambling the input data to provide DC-balancing of the frame, and then interleaving two Reed-Solomon encoded words (using 4-bit symbols), each capable of correcting a double symbol error (Figure 2). The interleaving operation allows increasing the correction capability of errors up to 4 symbols. This in practice means that a sequence of up to 16 consecutive incorrectly-received bits can be corrected. This correction technique requires an extra field of 32 bits in the frame to protect the 88 transmitted bits (including data, header and slow control), resulting in a code efficiency of 73%.

The frame (sketched in Figure 3) is composed of 120 bits that are transmitted during a single SLHC bunch crossing interval (25 ns) resulting in a line data rate of 4.8 Gb/s. Of these, 4 bits are used for the frame Header (H) and 32 used for Forward Error Correction (FEC). This leaves a total of 84 bits free for data transmission corresponding to a user bandwidth of 3.36 Gb/s. In these 84-bits, 4 are always reserved for the Slow Control (SC) field (see "Slow control channel") and 80bits are reserved for data transmission. Among the 4-bit of slow control, 2 are reserved for GBT control and 2 are user defined. The "D" field use is not pre-assigned and can be used indistinguishably for Data Acquisition (DAQ), Timing Trigger & Control (TTC) or Experiment Control (EC) applications [5][6].

This makes 2+80 = 82 bits of data available to the user for a frame of 120bits, giving a payload of 68%.

II. GBT PROTOCOL IMPLEMENTATION IN FPGAS

A. FPGAs constraints

In the same way it is done in the GBT ASIC, the DCbalance of data transmitted over the optical fiber is ensured by the FPGA by scrambling the data contained in the SC and D fields. For forward error correction the scrambled data and header are Reed-Solomon encoded before nibble interleaving, and serialization. The line encoding/decoding process is represented in Figure 4.

No problem was encountered to configure the hard-IP transceivers of the FPGAs, as the portability of this protocol was carefully checked during the specification phase of the GBT. In particular, the ability of Stratix and Virtex transceivers to transmit 120 bits at a frequency of 40 MHz was ensured at that time.

Figure 4: Block diagram of a full GBT link in an FPGA However, these transceivers provide neither specific encoding schemes like the one we selected nor flexible word alignment functions. This is mainly due to the fact that they target the most common telecommunication protocols. We had thus to implement in user logic all the encoding and decoding blocks, as well as a customized pattern detection and word alignment block (see Figure 5). At power on or after a loss of synchronization, the receiver starts a frame-lock acquisition cycle to find the frame boundaries, that is, to acquire frame synchronization.

The frame-lock acquisition mode operates as follows. In the StratixIIGX, the transceiver hard-IP word aligner block cannot be bypassed. It is thus configured to lock on an arbitrary pattern. Once completed the process is not repeated, except at power on or upon a command from the pattern detection state machine. For all the other devices, we bypass the word aligner inside the transceiver.

The parallel output of the receiver feeds the custom pattern detection and word aligner blocks, which take control of the frame alignment process: for each received frame the four bits in the header position are checked for header validity. Because the header pattern can be found in the data, 23 consecutive frames must contain a valid header before the frame is considered locked (the probability of false boundary detection is then reduced below 10 -20 as demonstrated in [START_REF] Papotti | Architectural studies of a radiationhard transceiver ASIC in 0.13 mm CMOS for digital optical links in high energy physics applications[END_REF]). Otherwise, the frame is shifted by one bit and the valid header checking procedure is repeated. After frame-lock is achieved, the receiver switches to the frame-tracking mode, which maintains frame synchronization even in the presence of headers corrupted by noise or single event upsets.

The phase tracking mode must thus be tolerant to a low rate of detection of invalid headers. Provided that frame synchronization is maintained, the detection of a corrupted header will not introduce a transmission error since the header field is also protected by the forward error correction code transmitted with the frame. A corrupted header will thus be corrected and properly identified by the Reed-Solomon decoder. The frame tracking mode operates as follows: after a successful frame-lock acquisition cycle has been executed the receiver enters the frame-tracking mode. In this mode the receiver strives to maintain frame synchronization. It checks the validity of the headers and counts the number of invalid headers received in 64 consecutive frames after the first invalid header has been detected. If the number of invalid headers received in 64 consecutive frames is bigger than 4 then the receiver re-enters the frame-lock acquisition mode. Otherwise the receiver resets the count of invalid frames and remains in the frame-tracking mode.

B. Resource Usage

The full serializer-deserializer, as described above, was implemented both in a StratixIIGX and in a Virtex5FXT. Besides the transceivers and PLLs, which do not consume any resources as they are hard-coded, a single link consumes 1542 ALMs (Adaptative Logic Modules) for the StratixII and 1481 Slices for the Virtex5.

The table 1 shows the number of links which can be implemented in a selection of StratixIIGX and of Virtex5FXT devices, taking into account the available transceiver blocks and logic elements. Differences of occupancy between Table 1 and Table 2 emphasize the different policies used by Altera and Xilinx in term of ratio between the number of logic cells and the number of transceivers. However, these numbers should be used with care. It is obvious that the occupancy of logic cells is too high if one tries to use all the available transceivers of a chip for GBT protocol implementation. This is tempered by the fact that a design using GBT links will not dedicate all its links to GBT transceivers: some links must be left to output processed data and therefore occupancy will be lower. However, as a back-end FPGA has to dedicate a significant part of its logic to other tasks, optimization of the resources used by the decoding block is a must.

C. Optimization

An analysis of the resource usage per block for a single link (see Figure 6) quickly shows that more than half of the logic elements are used by the Reed-Solomon decoder. It was thus natural to study optimization schemes, particularly for designs hosting several GBT links in one device. The first possibility is to share one decoder block between several links, multiplying its operating frequency by the same factor. The Reed-Solomon decoding algorithm is a large combinatorial circuit, and the maximum operating frequency achieved was 134MHz for the StratixIIGX, applying all the timing optimization constraints available. This allowed to share one decoder block between 3 links. An analysis of the resources used for 12 links implemented in a StratixIIGX type EP2SGX90 was carried out with and without optimization. As shown in the Figure 7, the device occupancy dropped from 51% of ALMs to 40% thanks to the optimization. Indeed, the fraction of the resources used by the decoder blocks dropped from 28% down to 10%. However, 7% of new logic elements were added due to the resource consuming multiplexers and de-multiplexers required to share the decoder.

This implementation was tested on a PCIe SIIGX development kit with three optimized links using loopback cables mounted on the HSMC connectors. It ran several days without a single error being detected. The next step for optimization could be to pipeline the decoder algorithm to increase the clock frequency. The drawback of this implementation, beside its complexity, is that it increases the decoding latency.

III. MEASUREMENTS

A. Setups and equipment

Two evaluation boards were used to implement the GBT protocol on FPGAs. The ML523 (hosting a Virtex5FXT type XC5VFX100T) for Xilinx [START_REF]Documentation on Xilinx Virtex5 FXT ML523 RocketIO GTX characterization Platform[END_REF], the PCIe SIIGX Development Kit (hosting a StratixIIGX type EP2SGX90) for Altera [START_REF]Altera PCI express Development Kit[END_REF], both powered by the power supply given in the kit (See Figure 8). The reference clock was generated by the J-BERT 4903A from Agilent on differential SMA cables.

For all the qualitative measurements, the very same SFP+ 1300nm optical transceiver module from MergeOptics was used (mounted and dismounted from one board to another). The optical patch cords were 50cm long.

The jitter measurements were made at the optical receiver level with the Lecroy SDA100G sampling scope equipped with 10 GHz optical sampling head.

B. Platform testing

Various platforms and technologies were tested by implementing the GBT protocol in both Altera and Xilinx chips presented above. As described on the Figure 9, a generator instantiated in the Virtex5 was sending parallel data (80 bits @ 40 MHz, either constant words or flying bits) to the encoder and serializer. The signal (that looks like a PRBS due to the scrambling) was transmitted by an SFP+ to the receiver in the StratixII over a short optical fibre (A). After full decoding (and remote monitoring of the decoded values), the data were encoded back, serialized again and transmitted using another SFP+ module and an optical fibre (B) back to the Virtex5, where it was decoded and compared to the generated words. We let the system run during several hours without counting any error. Besides providing us an opportunity to implement the GBT protocol on both main technologies, this test allowed us to check the compatibility between the GBT-ASIC protocol and its VHDL translation: the Virtex5 had the Reed-Solomon encoder and decoder implemented in Verilog (the direct copy of the GBT protocol implementation in the ASIC), whereas the StratixII encoder and decoder were implemented in VHDL.

C. Jitter performances

Using the same setup, we measured the jitter out of the two optical fibres A and B in Figure 9. For each of the results below, the SFP+ module transmitting the optical signal was the same (it was successively mounted on A and B fibres to test Xilinx and Altera devices).

As presented in Figure 10, Xilinx and Altera platforms both showed excellent performances. The eyes were widely open, and the total jitter of the order of 80ps PP and 5ps RMS.

IV. SOURCE CODE AVAILABILITY

Reference designs of the GBT protocol will be made available before the end of 2009 for both Altera and Xilinx FPGAs. They will be presented as a firmware-based starter kit, downloadable on request via the CERN SVN repository. This starter kit will include the source code for both implementations, and, as much as possible, for various types of devices (StratixII and IV GX, and Virtex5 and 6 FXT) and various flavors of optimization. It will also include documentation.

Basic support will be provided on how to use and optimize the implementation.

V. CONCLUSION

With this study, we proved that the GBT protocol can indeed be implemented with success both in Altera and Xilinx FPGA chips. The scheme proposed in the introduction where GBT ASICs are used in detector areas and FPGAs in counting rooms is thus a valid prospect, and the developed code will now be used as a basis to test the GBT serdes chip once it becomes available.

A firmware-based starter kit will be made available upon request to the users. It will be progressively completed by several implementation flavors for StratixIV and Virtex6, and new optimization techniques like a pipelined Reed-Solomon decoder are being considered.

Figure 2 :

 2 Figure 2: GBT encoding scheme

Figure 3 :

 3 Figure 3: GBT frame

Figure 5 :

 5 Figure 5: Frame alignment procedure in FPGAs

Figure 6 :

 6 Figure 6: % of ALMs/Slices of one GBT link used by each functional block

Figure 7 :

 7 Figure 7: Effect of optimization by 3 on 12 links implemented on a EP2SGX90

Figure 8 :

 8 Figure 8: Evaluation platforms. ML523 from Xilinx (left) and PCIe SIIGX from Altera (right)

Figure 9 :

 9 Figure 9: Test setup based on two platforms

Figure 10 :

 10 Figure 10: Eye diagrams for Xilinx Virtex5 FXT (left) and Altera StratixIIGX (right)

Table 1 :

 1 Maximum GBT links for StratixIIGX

Table 2 :

 2 Maximum GBT links for Virtex5FXT