Malik Khalfallah

Nicolas Figay
email: nicolas.figay@eads.netghodous

Parisa Ghodous

Catarina Ferreira
email: catarina.ferreira-da-silva@liris.cnrs.fr

Da Silva

Cross-organizational Business Processes Modeling Using Design-by-Contract Approach

Keywords: design-by-contract, interoperability, product model, business process, model driven architecture, UML

Reaching agreements between organizations in a collaborative environment is a way to ensure interoperability between these organizations at all levels. For business processes interoperability this agreement can be reached by well defining the cross-organizational process. However most BPM frameworks have used meta-models centered on flows of activities, with the data manipulated by these activities seen as second-class citizens. For business analysts (for example in complex product design collaborative environments) data plays a major role. In this paper, we propose a methodology backed by a conceptual framework to model the cross-organizational process relying on the product model. This framework defines the evolution of the product model through a finite number of states, and then automatically generates executable artifacts to support the collaboration during run-time phase. This approach is being implemented in the living laboratory provided by EADS in the context of the European project IMAGINE.

Introduction

The aeronautic and aerospace industries are seeing a rapid shift to the extended enterprise strategy. Companies specialized in this domain, including EADS, are increasingly externalizing activities for subcontractors while focusing on their core activities that are specific to their disciplines.

Ensuring interoperability between organizations in the extended enterprise is a complex problem and several studies have been conducted and reported in the literature. While many of these proposals focus on the static extended enterprise, where partners do not leave the network and are not replaced by new ones, ensuring sustainable interoperability for a dynamic extended enterprise has not been addressed well. We have recently seen the events of the volcano's eruption in Iceland and the nuclear disaster in Fukushima and their impact on enterprises. These events have reaffirmed the need for greater flexibility and maintaining organizations' interoperability and to cope with the dynamic nature of collaborative environments. To address this problem, the concept of Dynamic Manufacturing Network (DMN) has emerged [START_REF]Innovative end-to-end Management of Dynamic Manufacturing Networks (IMAG-INE Project)[END_REF].

In the context of product design collaborative environments, standards such as ISO-10303 (STEP) have been proposed in order to ensure partners' interoperability. The overall objective of STEP is to provide a mechanism that describes a complete and unambiguous product definition throughout the life cycle of a product. STEP provides both broadly useful data modeling methods and data models focused on specific industrial uses [START_REF] Kramer | STEP in a Nutshell, Advanced Design and Manufacturing Based on STEP[END_REF]. Even though STEP is an accepted standard in today's industry, it is still not sufficient to ensure a satisfactory level of interoperability in the extended enterprise. Indeed, relying on the classification of interoperability levels provided by Lewis et al. [START_REF] Lewis | Why Standards Are Not Enough to Guarantee End-to-End Interoperability[END_REF], we can say that STEP is limited to data interoperability of the first three levels (machine, syntactic, and semantic) but does not address data interoperability at the organizational level, as illustrated in figure 1, which remains the most complex one [START_REF] Lewis | Why Standards Are Not Enough to Guarantee End-to-End Interoperability[END_REF]. To fill this gap, EADS has participated to the European project CRESCENDO [START_REF]Collaborative and Robust Engineering using Simulation Capability Enabling Next Design Optimization[END_REF] that defined the concept of Behavioral Digital Aircraft Business Object Model (BOM). The BOM is a data model built with UML and defines high level concepts to be used by engineers pertaining to different organizations involved in the design and development of aircrafts. The BOM and its mapping to STEP concepts are the building blocks of the agreement necessary to ensure data interoperability at all levels.

In addition to data interoperability, another issue in the extended enterprise is the interoperability of organizations' business processes [START_REF] Papazoglou | Service oriented architectures: approaches, technologies and research issues[END_REF]. Cross-organizational processes aim to achieve an agreement between organizations at the process level. They specify what messages each organization should send and receive as well as their sequences. To model the cross-organizational process languages such as BPMN/XPDL (Business Process Modeling Notation/XML Process Definition Language) can be used. As illustrated in figure 1, XPDL covers the first three levels but not the fourth one. The reason is that when modeling the cross organizational process using XPDL, the business expert of each organization defines the exchanged messages based on their representation using STEP concepts (serialized in XML). Consequently, the high level concepts defined by the BOM are not usable and the business expert needs to understand the concepts defined by STEP. Moreover, the constructs provided by process modeling languages do not fit well for DMN processes due to the dynamic nature of DMN. They do not support well the changes on-the-fly that could occur [START_REF] Van Der Aalst | Dealing with workflow change: identification of issues and solutions[END_REF]. This analysis raised the following challenges: (i) how to define an agreement between participants that will ensure their business processes interoperability at the organizational level while using the concepts defined by the BOM? (ii) How to make this agreement flexible enough to maintain a sustainable interoperability of the whole DMN?

Contributions. Since the BOM defines the concepts used to design the whole aircraft from different point of views (i.e. physical, functional, system), we start from this data model and then, define the evolution of the product model properties along the collaboration until reaching the final configuration. The evolution is modeled using temporal logic relations between product model states' speci-fications. Each state describes what constraints should the product model fulfill and which organizations should ensure these constraints.

Using the proposed framework to model the cross-organizational process is promising because it decreases the coupling between organizations' processes.

In addition it provides the business expert the means to naturally define the cross-organizational process based on high-level concepts, since existing business process modeling approaches fail to capture the business intent underlying the interactions [START_REF] Telang | Business Modeling via Commitments[END_REF].

The rest of the paper is organized as follow: in section 2, we provide motivations regarding the proposed approach. In section 3 we give an overview of the methodology accompanied with an example. Section 4 elaborates on the conceptual framework that supports this methodology. Section 5 gives an overview of the implementation. Section 6 discusses the related work. Finally, section 7 concludes the paper and out-lines the perspectives of the current work.

Background

In this section we give a brief description of the building blocks of the proposed framework.

Model Driven Architecture

At the conceptual level, MDA is a holistic approach to improving the entire IT life cycle specification, architecture, design, development, maintenance and integration based on formal modeling [START_REF] Guttman | Real-Life MDA: Solving Business Problems with Model Driven Architecture[END_REF]. Today's MDA is less about generating code per se and much more about precisely capturing requirements, enforcing architectural standards, maintaining traceability, and facilitating effective communication between the business and IT.

From the analysis above, we figured out that we are facing the same issues that software developers faced in the domain of software engineering. One major issue is the lack of abstraction. Since MDA has proved to be successful as illustrated by industrial case studies demonstrated by [START_REF] Guttman | Real-Life MDA: Solving Business Problems with Model Driven Architecture[END_REF], or even internally in EADS, where it has been used in multiple projects, we carry on in this direction by using it. Our aim is not limited to generate executable code but also to generate more complex artifacts as explained below. Such an approach provides engineers with high level concepts when modeling the cross-organizational process without caring how this process will be implemented using workflow engines.

Design By Contract

Design by Contract (DbC) is an approach to building reliable software that focuses on making the contract of a software module explicit and formal [START_REF] Frankel | Model Driven Architecture: Applying MDA to Enterprise Computing[END_REF]. DbC is not limited to software design but can be used for high-level modeling [START_REF] Bocchi | A Theory of Design-by-Contract for Distributed Multiparty Interactions[END_REF]. DbC involves writing two kinds of formal constraints:

-Preconditions and post-conditions that are assertions about operations; -Invariants that are assertions about the system state that must be true, except during the execution of an operation.

In this paper we extend the core idea of DbC -software development through elaboration of type signature with logical predicates -to the design and modeling of cross-organizational process by linking constraints' specifications on product model properties using temporal logic relationships. This approach allows us to formalize an agreement between all organizations involved in the collaboration based on the product model. This is an important element to ensure interoperability at organizational level as noticed by Lewis et al. [START_REF] Lewis | Why Standards Are Not Enough to Guarantee End-to-End Interoperability[END_REF].

Business Rules

A business rule is "a statement that defines or constrains some aspect of the business. It is attended to assert business structure or to control or influence the behavior of the business" [START_REF] Van Eijndhoven | Achieving Business Process Flexibility with Business Rules[END_REF]. Business rules are executable by rules engine, thus they bridge the gap between contract constraints expressed at the organizational level and their implementation in the execution platform. A typical pattern of business rules is the ECA (Event Condition Action) notation [START_REF] Knolmayer | Modeling Processes and Workflows by Business Rules, Business Process Management, Models, Techniques, and Empirical Studies[END_REF]:

-The event component specifies when the rule has to be executed; -The condition component indicates a condition to be checked before any action is triggered; -The action states what has to be done depending on the result of the evaluation of the condition component. In general an action terminates by raising one or more relevant events.

We elaborate on the usage of business rules in the implementation overview in section 5.

Problem Statement

The general problem addressed in this paper is the need to develop a productbased technique to improve the change management and maintainability of crossorganizational business processes for DMN. Such a method should allow the product architect to predict and manage the impact of change on the processes used to develop the product, while decreasing the time needed to implement changes.

Methodology of Design by Contract for DMN

In this section, we present a methodology that defines how our approach can be used to model the cross-organizational process. This is illustrated in figure 2. Step 1. The Business Object Model is an ontology that defines the concepts used by all organizations in the DMN. It ensures semantic interoperability between all participants. It defines the concepts used to develop all product components.

Business
Step 2. The product model breakdown is a static model that defines the structure of the product. It structures the concepts defined in the BOM and establishes the relationships between these concepts. The product model is the building block of the contract, because all participants agree on this structure. However it remains not sufficient because it does not cover the behavior of the collaboration.

Step 3. Adding generic constraints to the classes' attributes of the product break-down allows the architects to formally define what values the properties should take. The objective of the collaboration is to iterate over the crossorganizational process until finding the optimal values of all product and sub products properties and defining the final configuration of the product. Relying on the classical relationships between constraints (e.g. and, or) is not sufficient to define the cross-organizational process. We need to add temporal relationships between constraints in order to define the succession of the product model states until reaching its final state.

Step 4. The logic-based model of the product model evolution is stable and it is in-dependent from: the organizations involved in the DMN, the processes enacted by these organizations and the IT systems that support the collaboration. Therefore, the last step adds organizations to the model and specifies what constraints each organization should fulfill.

Example

The following is an example of applying the proposed methodology using the results of the project CRESCENDO [START_REF]Collaborative and Robust Engineering using Simulation Capability Enabling Next Design Optimization[END_REF]:

Step 1. The concepts defined in the real BOM are complicated for non-aircraft engineers to understand. In this example we use usual concepts and simple attributes in order to demonstrate our solution.

Step 2. The product model of an aircraft is a tree-like structure. For a single aircraft, we can have several views that generate multiple decompositions (physical, functional, system etc.). Figure 3 illustrates this decomposition.

Step 3. Relying on the product model, the aircraft architect develops the crossorganizational process by defining the constraints on the aircraft functions (i.e. what functions the aircraft should fulfill). Additionally, he adds relationships between these constraints to indicate their dependencies as illustrated in figure 4. For this specific example, since a function is realized by physical components, the aircraft architect sets constraints on these components as well. This is illustrated in figure 5.

Step 4. Finally, each constraint is attached to two actors the requester, who sets the values of the attributes (the aircraft architect in our case) and the supplier (e.g. the wing designer) who tries to design the sub-product using these values.

At this stage the design of the contract is complete. The next step consists in iteratively running the collaboration between the involved partners until satisfying all constraints. Product architects define the desired evolution of the product along the collaboration by defining constraints and their relationships. Consequently, the focus is put on the product model which remains stable even if some organizations quit the network. This decreases the coupling between organizations' processes which ensures a sustainable interoperability of the whole collaborative environment.

For example, in figure 5, if the wing designer quits the DMN and he should be replaced by two separate designers: one to design the spoiler and another to design the flaps. This will have no impact neither on the evolution of the product nor on the remaining organizations in the DMN.

A Conceptual Framework for DbC

The objective of this framework is to back the presented methodology by formal basis. For this purpose, we define the concept of Obligation.

Definition 1. An obligation defines both what one party is obliged to guarantee and, dually, what other parties can rely on.

More specifically, each product component designed by an external organization needs to respect predicates over attributes representing its properties. These predicates are obligations for the component designer and a guarantee for all other organizations (e.g. the aircraft architect decides that the length of the wing shall be 20m. This is a guarantee for the aircraft architect and an obligation for the wing designer).

Contract Formalization

It is essential to have a tractable and rigorous representation of obligations. Formal representation of obligations will provide accurate and unambiguous specifications. For this reason we use logic formulae to formalize obligations. Logic formulae are a general, rigorous and flexible tool to describe constraints. Indeed using logic formulae, the product architect specifies the properties of the subcomponents composing the product and the relationships between these properties.

Obligations are the building block of the cross-organizational process. Relationships between obligations are important as well. Using temporal logic provides the means to define how the collaboration should evolve through different product model states.

Definition 2. A product model state is a set of obligations specifying constraints on a subset of objects defined in this model.

Two states are interconnected by a temporal logic relationship LeadsTo, as illustrated in figure 4, to express the precedence relationship between these states.

In this framework we purposely use temporal logic relationships patterns to inter-connect obligations and states. For instance, the knowledge required to use formal models and their complexity remain a significant obstacle for their widely adoption by business experts (i.e. product architects) [START_REF] Kharbili | Bringing Agility to Business Process Management: Rules Deployment in an SOA[END_REF]. Using patterns supports shielding the complexity of formalisms from business experts and facilitate their specification in the abstract. [START_REF] Dwyer | Property specification patterns for finite-state verification[END_REF][START_REF] Tretken | Enforcing compliance on business processes through the use of patterns[END_REF] already identified temporal logic patterns that we use to formalize the relationships between obligations and states. Given two obligations:O 1 , O 2 and two states: S 1 , S 2 :

Pattern Description O1inclusiveO2
The fulfillment of O1 mandates the fulfillment of O2 O1exclusiveO2

The fulfillment of O1 mandates the non-fulfillment of O2 O1prerequisiteO2

The non-fulfillment of O1 mandates the non-fulfillment of O2 O1mutexchoiceO2 Either O1 or O2 are fulfilled but not any of them or both of them

S1LeadsT oS2

State S1 must always be followed by state S2 Table 1. A set of patterns to interconnect obligations and states

Besides states and obligations, we add the concept of role. There are two types of roles:

-The requester: this is the partner who instantiate the obligations and waits to see whether the physical component can be designed under these obligations.

-The supplier: this is the partner who fulfills the obligations set by the requester.

Local Projections of the Contract

Modeling the cross-organizational process using obligations and product model states makes the local projections easy to generate. There are three kinds of artifacts to generate:

The requester process. The process of the organization that relies on obligations is illustrated in the right hand side of figure 6. For example, the aircraft architect and the wing designer collaborate in order to build an optimal model of the wing. The aircraft architect starts by giving a value to the flaps dimension.

After that the wing designer executes its business process, the aircraft architect receives the value of the new flaps dimension. This value can be equal to the already set value in the constraint which means that the wing designer has fulfilled its obligation in this case the aircraft architect carries on the execution of his internal process. Otherwise, the aircraft architect may give a new value to the constraint and repeat the process.

The supplier process. The process of the organization that fulfills a particular obligation has the pattern illustrated in the left hand side of figure 6. For example, the wing designer and the aircraft architect collaborate in order to ensure state 1 in figure 5. The wing designer receives the value of the flaps dimension, carries out an internal process to design the flaps following the set conditions, and then he notifies the aircraft architect of the new value.

The business rules that manage the update of the product shared model. Since in our implementation we are using a shared product model of the aircraft that is updated during the collaboration until reaching the final state, we decided to generate business rules that handle the updates of attributes' values and that generate notifications. Notice that even if this shared product model was absent, the generated artifacts would have been different but this remains an implementation detail and it has no impact on the way the crossorganizational process is modeled. processes interoperability issues in this case has been addressed in our previous paper [START_REF] Khalfallah | Product-based Business Processes Interoperability[END_REF]. Figure 7 gives an excerpt of the generated business rules of the given example. We rely on the business rules of type ECA as defined in [START_REF] Knolmayer | Modeling Processes and Workflows by Business Rules, Business Process Management, Models, Techniques, and Empirical Studies[END_REF].

Implementation

To illustrate our proposal, we have started to develop a prototype to test it with the results of the CRESCENDO project. The following steps illustrate how the proposed architecture (figure 8) implements the presented methodology with the conceptual framework. We should notice that there are two different phases:

ON Ensure_state1_event IF State1 = true RAISE Ensure_state2_event ON Ensure_state2_event IF State1 Leads_To State2 = true RAISE Ensure_state3_event
Fig. 7. Generated business rules the design-time phase where the aircraft architect defines generic constraints and assigns them to organizations. This defines the cross-organizational process (steps: 1, 2 and 3) and the run-time phase, where the cross-organizational process is executed until reaching the final configuration of the product model (step4):

Step 1. We assume that a product model is built and is being shared through a collaborative platform.

Step 2. The aircraft architect sets the (temporal) constraints on the shared product model through the collaborative platform. He uses an extension of the language OCL (Object Constraints Language) with temporal properties to define constraints on the classes in the UML model of the product. Additionally, he assigns corresponding organizations to each constraint. This step can be further enhanced by using advanced graphical user interfaces. At this stage we use an extension of OCL that supports temporal logic constructs and organizations assignment to constraints.

Step 3. Once the cross-organizational process is ready, the workflow generation module generates the workflows of the organizations and deploys them in their work-flow engines. Actually, a generated workflow is the public view on the complete organization's process. More specifically, it contains activities involved in the exchange of data with other organizations or with the collaborative platform.

The detailed design sub-process (illustrated in figure 5) is the internal process of each organization that contains private activities and is not shared with the external world. Additionally, in this step, the workflow generation module generates the business rules that handle update events that occur on the shared product model and verify the fulfillment of every obligation of every state.

Step 4. The aircraft architect collaborates with the sub-products designer through an iterative process until reaching a satisfactory configuration of the product.

We can notice that in the cross-organizational process model illustrated in figure 4, there is no exception handling (i.e. it models only the nominal process). This is purposeful, because we noticed that it is very difficult for an aircraft architect to identify all possible exceptions that can occur during a product design at the design time phase of the cross-organizational process. These exceptions are handled during the run-time phase. Our objective in this paper is to provide a framework supported by a formal theory that allows product engineers to model the cross-organizational process (i.e. the con-tract) using business concepts that are more abstract than low-level concepts such as data and control flows modeling constructs.

Modeling cross-organizational processes and defining their local projections is a widely studied subject: standardized languages WS-CDL, BPMN as well as Decserflow for choreography [START_REF] Montali | Declarative specification and verification of service choreographiess[END_REF], Scribble [START_REF] Honda | Scribbling Interactions with a Formal Foundation[END_REF], and Bocchi et al. proposal [START_REF] Bocchi | A Theory of Design-by-Contract for Distributed Multiparty Interactions[END_REF] focus on formally defining the cross-organizational process and propose a formal definition of the projection to generate participating organizations processes' public and private views. The common limits of these proposals regarding our problematic are twofold:

-They still use constructs such as messages, activities, flows that were demonstrated to be low-level concepts Telang et al. [START_REF] Telang | Business Modeling via Commitments[END_REF]. -The coupling maintained between organizations that could be difficult to handle when an organization quits the network and be replaced by a new organization.

Business rules and business processes. Eijndhoven et al, and Charfi et al. [START_REF] Van Eijndhoven | Achieving Business Process Flexibility with Business Rules[END_REF], [START_REF] Charfi | Mira: Hybrid web service composition: business processes meet business rules[END_REF] Work aims to increase the process flexibility by separating between the procedural flow of the process and the business rules. Their proposal is limited to a single enterprise process, because they assume that business rules are internal to an enterprise and not shared with the external world.

Product model and business processes. Van der aalst et al. [START_REF] Reijers | Product-Based Workflow Design[END_REF] proposed a method to design the optimal workflow based on the product model. The proposed methodology is limited to decision-making processes. They gave the example of a product decomposition of a candidate that wants to become a pilot. Relying on multiple criteria, they generated the optimal workflow that can decide whether a given candidate is able to become a pilot or not.

Declarative modeling of business processes. Van der aalst et al. [START_REF] Van Der Aalst | DecSerFlow: Towards a Truly Declarative Service Flow Language[END_REF] already pointed out the limits and the rigidity of procedural modeling languages such as BPMN and proposed a declarative language to model business processes. This work was used by Montali et al. [START_REF] Montali | Declarative specification and verification of service choreographiess[END_REF] in order to model the crossorganizational process. The limit of this proposal is the possible ambiguity when interpreting these models. While Van der aalst et al. [START_REF] Van Der Aalst | Business process management: a survey[END_REF] pointed out the mandatory nature of a process model to be unambiguous, these studies lack this property. Indeed, the modeled process is a temporal logic formula interconnecting activities. Inherently, they assume that the activity is true when it finishes, while there is a possibility to give a different interpretation: an activity is true when it starts. These two possible interpretations generate two different processes, the first one where two successive activities are executed in sequence and the second interpretation is that two successive activities are executed in parallel.

Business processes and web services. The interoperability of services' processes is also a widely studied subject and several researches are conducted in the domain of web services. Basically, these approaches generate adapters either automatically or semi automatically that could resolve both structural and behavioral mismatches between two communicating business processes. We have already defined an ontology-based approach to resolve structural mismatches that could occur due to the difference of STEP application protocols used by two communicating organizations [START_REF] Khalfallah | A Novel Approach to Ensure Interoperability Based on a Cloud Infrastructure[END_REF], [START_REF] Figay | Interoperability framework for dynamic manufacturing networks[END_REF]. For the behavioral mismatches we have decided to use the concept of DbC because mediation generation is an organization-dependent approach. If an organization quits the network we are obliged to regenerate a new mediator and redeploy it in the collaborative platform, which can be time consuming and expensive. Munusamy et al. [START_REF] Munusamy | A comparative study of process mediator components that support behavioral incompatibility[END_REF] provide a thorough state of the art on the techniques that generate adapters to resolve mismatches between two communicating business processes.

Conclusion

Imperative languages have proven to be non-trivial to support changes of the pro-cesses on-the-fly [START_REF] Van Der Aalst | Dealing with workflow change: identification of issues and solutions[END_REF]. This is a major obstacle to ensure business process interoperability in Dynamic Manufacturing Networks (DMN). In this paper we have presented a methodology backed by a conceptual framework to ensure business processes interoperability for DMN. We have used the Design by Contract (DbC) approach. DbC provides a declarative specification of what each organization in the DMN should accomplish (obligations) while maintaining this specification independent from organizations' specificities and the IT platform that supports the collaboration. Using the MDA architecture, we were able to project the DbC specification on the execution platforms and run the collaboration in order to fulfill the contract.

The results presented in this paper are being implemented in the living lab of EADS provided in the context of the European project IMAGINE [START_REF]Innovative end-to-end Management of Dynamic Manufacturing Networks (IMAG-INE Project)[END_REF].

Fig. 1 .

 1 Fig. 1. Proposal overview

Fig. 2 .

 2 Fig. 2. Proposal overview

Fig. 3 . 2 Fig. 5 .

 325 Fig. 3. Proposal overview

Fig. 6 .

 6 Fig. 6. Process Pattern for both sides of collaboration

Fig. 8 .

 8 Fig. 8. Implementation overview


ABSTRACT:

Today several lidar networks around the world provide large data sets that are extremely valuable for aerosol and cloud research. Retrieval of atmospheric constituent properties from lidar profiles requires detailed analysis of spatial and temporal variations of the signal. This paper presents an algorithm called STRAT (STRucture of the ATmosphere) designed to retrieve the vertical distribution of cloud and aerosol layers in the boundary layer and through the free troposphere and to identify near particle free regions of the vertical profile and the range at which the lidar signal becomes too attenuated for exploitation, from a single lidar channel. The paper describes each detection method used in the STRAT algorithm and its application to a tropospheric backscatter lidar operated at the SIRTA observatory, in Palaiseau, 20 km south of Paris, France. STRAT retrievals are compared to other means of layer detection and classification; retrieval performances and uncertainties are discussed.

Introduction

Remote sensing of the atmosphere by Lidar dates from the 1960's. Lidars have since become very popular remote sensing tools. They are used extensively both in ground-based and airborne configurations to get detailed vertical profiles of cloud or aerosol properties, such as the extinction coefficient to learn more about the nature of particles present aloft. Lidars are commonly deployed during intensive observation field experiments well as for routine long-term monitoring.

Today several global and regional networks of atmospheric observatories exist that use lidars as their main monitoring instrument. The European Aerosol Research Lidar Network (EARLINET) monitors aerosol transport over Europe based on 21 Rayleigh/Mie and Raman lidar stations (Boesenberg et al. 2003). Similarly the National Institute for Environmental Studies of Japan coordinates a network of 12 automatic dual wavelength polarization lidars continuously to study Asian dust transport [START_REF] Shimizu | Continuous observations of Asian dust and other aerosols by polarization lidars in China and Japan during ACE-Asia[END_REF]. The continent between Europe and Asia is covered by the CIS-LiNet project that aims at aerosol and ozone monitoring through a network of 6 lidar stations since the beginning of 2004. Building from the research capabilities already established at a number of Eastern North America lidar facilities, a Regional East Atmospheric Lidar Mesonet (REALM) has been proposed to monitor air quality in the vertical from multiple locations in that region [START_REF] Hoff | A Regional East Atmospheric Lidar Mesonet (REALM)[END_REF]. The National Aeronautics and Space Administration funds a Micro-Pulse Lidar Network (MPLNet ; [START_REF] Welton | Global monitoring of clouds and aerosols using a network of micro-pulse lidar systems[END_REF] to monitor troposheric aerosols based 20 micropulse lidars located at climatologically diverse locations on the planet and coordinated with the Aerosol Robotic Sunphotometer Network (AERONET). The Network for the Detection of Stratospheric Change (NDSC) monitors stratospheric ozone using lidars that range through the stratosphere.

Continuous or routine operations of such systems produce upwards of 1000 vertical profiles per day for each channel of each system. All systems produce fundamental information about the structure of the atmosphere, that is the vertical distribution of the particle layers from near the ground to the top of the lidar range. Those are basic parameters, yet our current understanding of the vertical distributions of cloud and aerosol layers in the atmosphere remains limited by lack of large-scale analysis of available data sets. Furthermore, in the context of new satellite missions carrying active remote sensing payloads, such as the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO, Winker et al., 2003), ground-based lidar observatories have an increasingly important role to play in establishing regional climatologies that will tie in the temporally sparse global-scale satellite measurements.

The need for robust algorithms designed to process large lidar data sets has been established for over a decade (e.g. [START_REF] Platt | The Experimental Cloud Lidar Pilot Study (ECLIPS) for Cloud-Radiation Research[END_REF]. Several authors have exploited long-term lidar datasets to derive statistics and climatologies of cloud and aerosol macrophysical properties (e.g. Cadet et al., 2003;[START_REF] Comstock | Ground-based lidar and radar remote sensing of tropical cirrus clouds at Nauru Island: Cloud statistics and radiative impacts[END_REF]Sassen et al., 2001). In such analyses, particle layers, being clouds or aerosols, are typically detected using slope change in the lidar profile (e.g. [START_REF] Flamant | Lidar determination of the entrainment zone thickness at the top of the unstable marine atmospheric boundary layer[END_REF][START_REF] Shimizu | Continuous observations of Asian dust and other aerosols by polarization lidars in China and Japan during ACE-Asia[END_REF], or comparing the lidar power return to an expected clear-sky value (e.g. [START_REF] Clothiaux | An Automated Algorithm for Detection of Hydrometeor Returns in Micropulse Lidar Data[END_REF]. Several authors developed methods using wavelet transforms to identify particle layers in the lidar profile (e.g. [START_REF] Cohn | Boundary-layer height and entrainment zone thickness measured by lidars and wind pro-filing radars[END_REF][START_REF] Brooks | Finding Boundary Layer Top: Application of a Wavelet Covariance Transform to Lidar Backscatter Profiles[END_REF], but those studies were limited to boundary layer height analyses.

In this paper we present an end-to-end algorithm designed to retrieve an ensemble of basic parameters of the particulate atmosphere from single channel lidar profiles. The algorithm is designed with modules that retrieve the height of the boundary layer, the vertical distribution of particle layers (clouds and aerosols), and that identify layers that are near particle free. Each retrieval taken independently is not new, however, the added value comes from the use of multiple analyses on each lidar backscatter profile to derive a self-consistent classification. The specifications of this algorithm are:

• To detect all particle layers from near ground to the top of the lidar range, that is, cloud and aerosol layers in the free troposphere and the vertical extent of the boundary layer from the raw uncalibrated lidar backscatter profile.

• To identify pristine layers predominantly populated by molecules, and hence considered as near particle free, so that the lidar signal can be calibrated (with respect to a calculated molecular backscatter) in an automated manner for further inversion of the signal.

• To be widely applicable to any single-or multi-wavelength ground-based lidar system (micro-to milli-joule pulse energy, high or low pulse frequencies, U.V., visible, near-infrared wavelength) by using the signal-to-noise ratio as only threshold reference.

• To be automatic and robust so that long lidar time series can be processed in an operational environment.

The objective of this algorithm development is to provide a tool that can be used to process large ground-based lidar datasets and hence give access to a very detailed macrophysical classification of the particulate atmosphere. This is the first and necessary step towards extending regional statistics of cloud and aerosol macrophysical properties.

In Section 2, we present the lidar measurement data set used to develop the STRAT algorithm. Section 3 is dedicated to the description of the algorithm and Section 4 to the evaluation of the algorithm. Conclusions are provided in Section 5.

Data description

To develop the STRAT algorithm we use a database produced by the SIRTA 1 observatory, a facility dedicated to observing the atmosphere in support of cloud and aerosol research. SIRTA is located in Palaiseau, 20 km south of Paris, France; the data can be retrieved through its web site (www.sirta.fr). SIRTA gathers active and passive remote sensing instruments to retrieve optical, radiative and dynamic properties of the atmosphere and its constituents. This includes a dual wavelength polarisation lidar LNA Lidar for Lidar Nuages Aérosols, that is Cloud Aerosol Lidar), a millimetre wave Doppler radar, a near-IR ceilometer, a surface broadband flux station, and standard weather measurements [START_REF] Haeffelin | SIRTA, a ground-based atmospheric observatory for cloud and aerosol research[END_REF]. The SIRTA database includes also radiosonde profile data produced by Météo-France 15 km from the site as part of their national operational network (00 and 12 UT).

The main dataset used to develop STRAT is that of the LNA lidar. Table 1 lists its technical characteristics. LNA lidar profiles corrected for electronic noise, atmospheric background signal and range divergence are available in the SIRTA data base in netcdf format for the period starting March 2002 through today (3000 hours).

1 Site Instrumental de Recherche par Télédetection Atmosphérique These level-1 files contain lidar backscatter profiles for 6 different channels (see Table 1) along with a quality flag that indicates when signal-to-noise ratio drops below 3 (see Section 3.1). Several studies [START_REF] Hodzic | Comparison of aerosol chemistry transport model simulations with lidar and Sun photometer observations at a site near Paris[END_REF][START_REF] Naud | Assessment of MISR and MODIS cloud top heights through comparison with a back-scattering lidar at SIRTA[END_REF][START_REF] Haeffelin | SIRTA, a ground-based atmospheric observatory for cloud and aerosol research[END_REF] are based on this lidar dataset and the output of the STRAT algorithm described in this paper. The lidar signal analysis by STRAT makes use of density profiles derived from radiosounding data (temperature and pressure) to simulate the theoretical lidar backscatter solely due to molecules. Two radiosonde profiles per day are enough to capture most density profile variations.

In section 4 we use ceilometers measurements collocated with the LNA lidar to perform comparisons of cloud detection and cloud base height retrievals. The ceilometer uses a laser diode with an emitted wavelength of 855nm, its pulse repetition frequency is 6494 hz with maximum emitted light power of 50 µW.m -2 and its optical divergence is 1.2 mrad. Cloud information is retrieved through a Vaisala proprietary algorithm that looks for high backscatter values in the profile.

In section 4 we also compare STRAT retrievals of the boundary layer height with estimations based on potential temperature profiles. For that we derive highresolution (10 m) temperature and pressure profiles from the standard significant points provided in the Météo-France file.

Algorithm description

The STRAT algorithm includes four successive detections carried out on individual lidar profiles:

• noise detection, to identify where the lidar signal becomes too weak (signal-tonoise ratio) for further exploitation, so that data is not misinterpreted.

• molecular layer detection, to identify regions that are near particle free, i.e. populated by molecules only or with a negligible concentration of particles, because lidar calibration is carried out in those regions by comparison with computed molecular backscatter profiles. This will enable lidar data users to automate processes that require calibration such as inversion of extinction profiles or calculation of layer optical depth or investigation of cloud thermodynamic phase using depolarization ratio.

• particle layer detection with separate cloud and aerosol layer identification, so cloud and aerosol processes can be studied

• boundary layer detection, so boundary layer process studies can be addressed.

Finally, a flag variable is derived with value 0 for noise, 1 for molecular layers, 2 for boundary layer, 3 for aerosol layers, 4 for cloud layers and 10 for unidentified layers.

Figure 1 shows a diagram of the STRAT algorithm data processing.

a. Noise detection 1) Method

A simple signal-to-noise ratio threshold allows us to determine where the signal is too noisy to extract information from the lidar measurement. For each profile, it has been shown [START_REF] Durieux | Measurement of the lidar signal fluctuation with a shot-per-shot instrument[END_REF]) that the noise level σ(r,t) of the measured signal, assuming all its components are independent, can be written as σ r , t=  σ sig r , t 2 σ b r ,t  2 σ d r , t 2 , 1 where σ sig (r,t) 2 is the shot noise induced by and proportional to the backscattered lidar signal, σ b (r,t) 2 is the shot noise resulting from background light and σ d (r,t) 2 is the shot noise resulting from dark current.

Assuming that σ sig (r,t) 2 , σ b (r,t) 2 and σ d (r,t) 2 are proportional to the backscattered lidar power P(r,t), to the background noise power B d (r,t) and to the dark current noise power B p (r,t), respectively, equation 1 becomes

σ r , t=C  P r , tB b r , t B d r ,t  , 2
where C is a proportionality coefficient to be determined. Assuming that B(r,t), defined as the sum of B d (r,t) and B p (r,t), is constant along the profile, the noise level can be estimated at the altitude range where there is no lidar return (P(r,t)=0) by computing the signal standard deviation σ P=0 (r,t) of P(r,t) at this range.

The proportionality coefficient C can thus be estimated by averaging the ratio derived from equation 2 over a given number of points, as

C= σ P=0 r , t   B r , t  = 1 N ∑ 1 N σ P=0 r , t   B r , t  . 3
The noise level can then be expressed as σ r , t= σ P=0 r , t  B r , t   P r , t B r , t .

4 So the signal-to-noise ratio SNR(r,t) can be written as

SNR r , t = P r , t  σ P=0 r , t   B r ,t   P r , tB r , t  . 5
Note that C can also be derived from pre-trigger data if available. This retrieval can be checked against the result of equation 3 for consistency.

2) Threshold determination

In the STRAT algorithm, the signal is considered to be too noisy for further analysis when the SNR falls below a threshold of T 1 =3. Indeed, for a gaussian noise 99% of values are contained in the interval +/-3σ P=0 (r,t). Boundary layers, molecular layers, and cloud and aerosol layers will be detected on the part of the signal that is above that threshold. For systems with very low signal to noise ratio, the algorithm must be applied to time averaged profiles.

Note that if B(r,t) is large (e.g. daytime) σ P=0 (r,t) can be used as the signal noise for convenience. This simplification does not introduce significant error in determining if SNR(r,t) is less or greater than 3, as for small SNR values σ (r,t) is quite close to σ P=0 (r,t) because P(r,t)/B(r,t) tends toward 0.

b. Molecular layer detection

1) Method

Identification of particle-free or near particle-free layers is of particular importance as they are often used in lidar calibration algorithms (e.g. [START_REF] Sassen | Optical Scattering and Microphysical Properties of Subvisual Cirrus Clouds, and Climatic Implications[END_REF][START_REF] Platt | Remote Sounding of High Clouds: I. Calculation of Visible and Infrared Optical Properties from Lidar and Radiometer Measurements[END_REF]. For simplicity, these layers will be labelled as molecular layers although they may contain aerosols in small quantities (see Section 4.a. for discussion). The identification algorithm for molecular layers is based on analysis of the variability of the lidar signal around a theoretical molecular backscatter profile computed from pressure and temperature profiles. Thermodynamic profiles can be obtained from collocated atmospheric sounding measurements or extracted from model analysis data.

A normalization coefficient K(r,t) is estimated at each range of the lidar signal as

K r ,t = ∑ r '=r-N r N β mol r ', t  β lidar r ', t  2N1 6
where

β lidar r , t=P r , t r 2 7
is the non-normalized attenuated backscatter coefficient, β mol (r,t) is the computed theoretical molecular backscatter coefficient and 2N+1 is the size of the averaging window.

The variability, V(r,t), of the normalized lidar signal around the molecular backscatter profile at a range r is determined as

V r , t = ∑ r '=r-N r N  1 r ' 2 [ β lidar r ',t -1 K  r , t  β mol r ', t  ]  2 2N1 , 8
where K is the coefficient used to normalize the lidar profile to the molecular profile in the averaging window. A sensitivity study showed 20 gates (e.g. 300m for the LNA lidar) to be a wide enough averaging window to detect only atmospheric variations.

In molecular layers, the lidar backscatter signal β lidar (r,t) can be expressed as the sum of molecular backscatter P' mol_lidar (r,t) and an additional zero mean noise Μ(r,t) :

β lidar r , t= P ' mol lidar  r , t M  r ,t  r 2 . 9
Hence, in molecular layers Eq. 8 becomes

V r , t = ∑ r '=r-N r N M r ', t 2 2N1 . 10
The variability V(r,t) is only due to the noise variability and so can be compared to the noise variance σ(r,t) 2 . So, a threshold value V thr (r,t) can be defined with respect to the noise variance σ (r,t) 2 as

V thr  r , t =T 2 σ  r ,t  2 11
where Τ 2 is the molecular layer threshold coefficient. Thus, if V(r,t) is below the threshold value V thr (r,t), we consider that the lidar backscatter power is characteristic of a molecular layer.

2) Threshold determination

Here again we use σ P=0 (r,t) as a substitute for σ(r,t) because in molecular layers the two values are quite close. Indeed, in our case P(r,t) is typically lower than the background backscatter B(r,t) and hence σ(r,t) tends toward σ P=0 (r,t).

To determine this threshold, values of V(r,t)/σ P=0 (r,t) 2 (Eq. 8 and 11) have been computed from one profile per day over the available database (October 2002 -October 2005) only on signal values with SNR greater than 3.

Figure 2 illustrates the probability density function (PDF) of V(r,t)/σ P=0 (r,t) at two wavelengths (532 and 1064 nm) with 2N+1=21 gates = 315m. The PDF of the function V(r,t) expressed in Eq. 10 is also represented in Figure 2 with a solid line.

The simulated noise used to estimate this PDF is a Gaussian noise similar to the real one.

Distributions are divided in two separate regions. The first one is a narrow Gaussian distribution ranging between 0 and 3. It can be associated with molecular layers where the variability is smaller than in particle layers. The second one is a very broad distribution of V(r,t)/σ P=0 (r,t) 2 with values greater than 3. Because the objective is to use molecular layers for calibration, it is important not to falsely detect stable particle layers as molecular layers. For our application, a threshold value V thr (r,t) of 3σ P=0 (r,t) 2 is a good compromise to separate molecular layers from other layers. The difference in distribution width between the 532 and 1064nm curves of Figure 2 can be attributed to the signal quality of the two different channels.

c. Cloud and aerosol layer detection

1) Method

The majority of particle layer detection techniques described in the literature use thresholding tests on the first derivative of the backscatter intensity (e.g. [START_REF] Pal | Automated method for lidar determination of cloud base height and vertical extent[END_REF]. Such methods give satisfying results as long as the signal-to-noise ratio remains high. Others techniques use algorithms that depend on cloud type (e.g. [START_REF] Chazette | Determination by Spaceborne Backscatter Lidar of the Structural Parameters of Atmospheric Scattering Layers[END_REF]). While they are suited for case studies, they cannot be used for automated detection. In the STRAT algorithm, a combination of wavelet transform and Pr 2 ratio thresholding is used to identify particle regions in lidar profiles.

The Continuous Wavelet Transform (CWT) is used to detect discontinuities in the lidar signal as the base, the top and the peak backscatter of individual particle layers.

This method, based on seeking high correlation between the lidar signal and the wavelet characterized by the "Mexican hat" shape for each range and for each scale, is inspired on studies by [START_REF] Mallat | Singularity Detection And Processing With Wavelets[END_REF] and an algorithm developed by [START_REF] Brooks | Finding Boundary Layer Top: Application of a Wavelet Covariance Transform to Lidar Backscatter Profiles[END_REF]. The "Mexican hat" wavelet Ψ(r), shown in Figure 3, is the second derivative of a gaussian. It is used because its shape is very similar to the shape of the lidar signal backscattered by cloud or aerosol layers. Additionally, "Derivates of Gaussians are most often used to guarantee that all maxima lines propagate up to the finest scales", which is not the case of the Haar wavelet [START_REF] Mallat | Singularity Detection And Processing With Wavelets[END_REF].

Firstly, the CWT is computed for each P(r,t) profile as

Ψ a , b  r = 1 a Ψ  r-b a  , 13
where a is the wavelet dilation (or scale) and b is the location of its center. CWT coefficients can be interpreted as a correlation coefficient between the wavelet (centred on b and scaled by a) and the signal P(r,t).

Secondly, the modulus of CWT coefficients are determined to extract the lines of modulus maxima of the CWT a,b (r,t) that are lines (or ridges) formed by all maxima found at all dilations. This skeleton of the CWT, formed by all ridges, represents the highest correlation and anti-correlation between the signal P(r,t) and the wavelet from the largest to the finest scale.

Figure 4 illustrates the performance of this method for an ideal cloud or aerosol backscatter profile (Figure 4a) and a real one (Figure 4d). Figure 4b and Figure 4e show CWT coefficients for the ideal profile and the real one, respectively. Figure 4c and Figure 4f show the corresponding maxima lines. Ridges of highest correlation and anti-correlation coefficients propagate to the finest scale at the base and top of each particle layer, as well as at the location of the maximum backscatter.

So, each ridge shown in Figure 4c is associated with a discontinuity of the P(r,t)

signal. The value M CWT (i ridges ) of the average CWT coefficients along this ridge allows us to discriminate a backscatter peak from a layer base or top as M CWT  i ridg =CWT a , b  r ,t | a , b∈i ridg { ¿ 0 : layer peak ¿ 0 : layer base or top 14

For each identified backscatter peak, the base (top) of the same layer can be found by looking for the first base (top) detected below (above). If the top of one layer is the base of the next one, the STRAT algorithm is designed to link these two layers into a single one with a peak defined as the maximum P value of the two original peaks.

Finally, we apply a threshold value R thr on the difference of backscatter power between peak height and base height defined as

R=P  r peak -P r base R thr .
15

This threshold removes over-detections that are due to noise variations such as the discontinuities detected between points 200 and 350 shown in Figure 4f. The R thr threshold implemented in the STRAT algorithm is derived with respect to the noise level σ(t) as

R thr =T 3 σ t  16
where T 3 is the particle layer threshold coefficient.

2) Threshold determination

As this threshold is used to identify false peak/base detections in layers with low backscatter signal (i.e. molecular layers) we use σ P=0 (r,t) as a substitute for σ(r,t).

A PDF of R(r,t)/ σ(t) values is derived from the LNA lidar database (10 profiles per day) to determine the R thr . The PDF is shown in Figure 5. Because of noise-related signal variations, discontinuities can also be detected in molecular layers but corresponding R(r,t) values are smaller than for particle layers. The PDF of R(r,t)/ σ(t)

values for discontinuities identified in molecular layers are also represented in Figure 5 (solid line). This curve is obtained by processing simulated noisy molecular profiles derived from radiosonde data. The distribution is divided in two separate regions. The gaussian distribution between 0 and 10 that contains 86% of the detections is due to noise variations. This effect corresponds to the many short CWT ridge lines shown in Figure 4f. Values R(r,t)/ σ(t) for true particle layers are logically greater. Picking a threshold value R thr of 10σ(t) allows us to remove over-detections.

d. Cloud and aerosol distinction

1) Method

Aerosol and cloud layers can have similar signatures in lidar backscatter profiles.

However, for near I.R, I.R, visible and U.V. wavelengths, the lidar backscatter power is generally greater for liquid-water and optically thick ice for clouds than for aerosols.

The cloud and aerosol distinction algorithm is based on the study by [START_REF] Wang | Cloud Type and Macrophysical Property Retrieval Using Multiple Remote Sensors[END_REF] who applied a threshold on the peak Pr 2 to base Pr 2 ratio. The ratio is expressed as

d Pr 2 = P r peak , t r peak 2
P r base , t r base 2

. 17 Ratios greater (less) than a threshold T 4 classify a layer as cloud (aerosols). Figure 6a shows a 7-hour time series of 532nm backscatter power profiles. The measurements show significant backscatter between ground and 2500m. After 11.30

UT one can see several occurrences of very strong extinction that are characteristic of dense water clouds. Figure 6b shows profile-by-profile dPr 2 ratios for the areas identified as particle layers. It reveals a large profile-by-profile variability of dPr 2 values.

To improve this method, we derive average dPr 2 values for a given object (cloud or aerosol layer). To obtain this averaged value on a vertically and temporally consistent particle layer, range-time processing is required.

An average dPr 2 layer value is computed for each identified particle layer and the T 4 threshold is applied to this average value to separate cloud from aerosol layers as

{ d

e. Boundary layer height detection

1) Method

The atmospheric boundary layer (ABL) is the lowest part of the troposphere which is directly influenced by the Earth's surface and responds on short time scales to surface forcing. This is the region which is well-mixed due to convectively driven mixing. Several boundary layer height (BLH) detection methods are described in the literature. Methods using a simple signal threshold (e.g. [START_REF] Melfi | Lidar observations of the vertically organized convection in the planetary boundary layer over the ocean[END_REF][START_REF] Boers | Lidar observations of the fine-scale variability of marine stratocumulus clouds[END_REF] are not appropriate for cases with varying aerosol extinction. Methods based on gradient properties at the top of the boundary layer (e.g. [START_REF] Flamant | Lidar determination of the entrainment zone thickness at the top of the unstable marine atmospheric boundary layer[END_REF]) need averaged or smoothed signals and so lose resolution. In the presence of boundary layer clouds, all methods, including wavelet analysis methods (e.g. [START_REF] Cohn | Boundary-layer height and entrainment zone thickness measured by lidars and wind pro-filing radars[END_REF][START_REF] Brooks | Finding Boundary Layer Top: Application of a Wavelet Covariance Transform to Lidar Backscatter Profiles[END_REF]) are likely to identify the top of the cloud as BLH because the strongest gradient (or correlation) will occur in that part of the profile. In the STRAT algorithm we use the output of the molecular layer module and particle layer module to help distinguish the low-altitude clouds from the boundary layer below them. The boundary layer height detection method used in the STRAT algorithm, is similar to the particle layers detection method described in section c. It is inspired by the work of [START_REF] Mallat | Singularity Detection And Processing With Wavelets[END_REF] and [START_REF] Brooks | Finding Boundary Layer Top: Application of a Wavelet Covariance Transform to Lidar Backscatter Profiles[END_REF].

The wavelet used here is the first derivative of a gaussian Ψ'(r), shown in Figure 8 because its shape is very similar to the negative gradient of the backscatter signal at the top of the boundary layer during daytime. A standard boundary layer backscatter signal is shown in Figure 9a. As for particle layer detections, the CWT is computed for each P(r,t) profile as

CWT ' a , b r , t= ∑ r P r ,t  ψ ' a , b r  20 where Ψ ' a , b  r = 1 a Ψ '  r -b a 
21 a is the wavelet dilation (or scale) and b is the location of its center.

Then modulus maxima lines of the CWT' a,b (r,t) are also determined to detect all gradients in the backscatter signal. Because of the wavelet shape, negative gradients can be discriminated from positive ones using average values of the CWT' coefficients along this ridge, as follow:

M ' CWT i ridg =CWT ' a , b r , t | a , b ∈i ridg { ¿ 0 : positive gradient ¿ 0 : negative gradient 22
The detection of negative gradients combined with the altitude of the lowest molecular range, H min_mol (shown in Figure 9), and the base height of the lowest particle layer (H min_part ) allows us to estimate the boundary layer height. Four different cases, illustrated in Figure 9a,c,e, must be considered:

if H min_mol < H min_part (molecular layer below the lowest identified particle layer)

o there exists a ridge with M' CWT <0 that propagates up to a range r<H min_mol : BLH is the range r (if there is more than one ridge, only the ridge with the minimum M' CWT value is kept). This case is illustrated in Figure 9a with an example of a standard lidar backscatter signal and in Figure 9b with the corresponding wavelet coefficients.

o there does not exist a ridge with M' CWT <0 that propagates up to a range r<H min_mol : BLH is undefined if H min_part < H min_mol (molecular layer above the lowest identified particle layer)

o there exists a ridge with M' CWT <0 that propagates up to a range r<H min_part : BLH is the range r (if there is more than one ridge, only the ridge with the minimum M' CWT value is kept). A cloud or aerosol layer is located near the top of the BL. This case is illustrated in Figure 9c with an example of a backscatter signal where the molecular layer is above the lowest identified particle layer. Figure 9d shows with the corresponding wavelet coefficients.

o there does not exist a ridge with M' CWT <0 that propagates up to a range r<H min_part : BLH is H min_part . A cloud or aerosol layer is located at the top of the BL. This case is illustrated in Figure 9e with an example of a backscatter signal where a cloud layer is present at the top of the BL.

Figure 9f shows with the corresponding wavelet coefficients.

After daytime convection ceases, aerosol layers become stratified and multiple layers can form near the surface (boundary and residual layers). In such situations, the STRAT algorithm is not able to distinguish the top of the boundary layer and the top of the residual layer.

Evaluation of the STRAT algorithm a. Evaluation of the molecular layer detection

Here we evaluate if layers identified by STRAT as molecular layers contain any additional extinction due to the presence of some quantity of aerosols. To do so we apply a classic approach of optical thickness estimation (Platt, 1973) that is based on the ratio of the lidar power attenuation from the base to the top of the molecular layer to a theoretical molecular attenuation. Analysis of 4 years of SIRTA lidar profiles containing molecular layers extending more than 1 km reveals that these layers exhibit attenuation uncertainties of +/-2.10 -5 m -1 in terms of equivalent extinction.

Hence the parameters used in the molecular layer detection module (see Table 2)

imply that STRAT will allow layers whose attenuation is somewhat different from that of theoretical molecular layers to be identified as particle free. As a result such layers could contain up to 2.10 -5 m -1 particle extinction, equivalent to a 0.02 optical depth for a 1 km deep layer. This uncertainty can be reduced by increasing the test range (e.g. from 300 to 500 m) or by reducing the variability threshold (T 2 ).

b. Cloud and aerosol layer detection

1)

Performance evaluation based on simulated data

Figure 10 shows results obtained by the STRAT algorithm cloud and aerosol layer detection with a simulated backscatter profile containing a cloud (Figure 10a). Two slopes S 1 and S 2 are used to describe the majority of cases, where S 1 is the molecular slope and S 2 is the slope of the backscatter profile in the cloud between the base and the peak. Figure 10b and c illustrate results obtained on CBH and CTH detection, respectively. We describe results obtained for slopes between -0.5.10 -10 and -2.10 -10 m -1 .sr -1 /m for S 1 and between 1.10 -8 and 7.10 -8 for S 2 . The CBH detection is sensitive to both slopes but the maximum resulting error is -3 gates (-45m for LNA lidar profiles). CTH detection depends on the slope of backscatter in the particle layer. CTH errors are biased high between 0 and 5 gates (0 to 75m for LNA lidar data) for the largest S 2 values.

2) Comparison with retrievals from ceilometer

We compare cloud base height retrievals derived by applying STRAT to LNA lidar data to those derived by a VAISALA ceilometer located nearby (100m). Figure 11a shows the PDF of cloud base height for the two systems based on 12 months of observation. We limit our comparisons to situations where each retrieval is consistent for 10 minutes. We note that the LNA misses clouds below 1300m, due to the large overlap function. Above 5000m, the ceilometer data becomes unreliable due to limited power. Comparisons of cloud and cloud free occurrence detections by the two systems for the 1300-5000m vertical range are shown in Table 3. In situations labelled as cloud free by the ceilometer, we find 92% agreement, and 8% cloud detection by LNA/STRAT. In situations where the ceilometer detects a cloud between 1300 and 5000m, LNA/STRAT detects a particle layer 93% of the time with 74%

clouds and 19% aerosols. The cloud vs aerosol discrepancy can result from the simple cloud/aerosol threshold used in STRAT as well as possible aerosol detection by the 855nm ceilometer. Next we compare CBH retrieved by both systems when they agree that a cloud is present in the 1300-5000m range. Figure 11b shows a scatter plot of LNA CBH versus ceilometer CBH and Figure 11c shows the PDF of the difference between the two retrievals, based on 12 months of observation. The VAISALA CBH detection method (Vaisala proprietary algorithm) is based on the detection of high backscatter in the profile so the retrieved CBH is frequently placed between the base of the cloud and the altitude of maximum backscatter in the cloud.

The position of the ceilometer CBH in the lidar backscatter is illustrated in Figure 12.

The mean difference between LNA/STRAT and ceilometer CBH is -178 m which is consistent with the result of Figure 12. The standard deviation of the comparison is 265 m. The PDF can be divided into three zones: zone 1 contains 80% of detections and gather the situations with most consistent retrievals, zone 2 includes 5% of the situations for which ceilometer CBH are lower than corresponding STRAT retrievals, and zone 3 contains 14% of the distribution for which STRAT retrievals are lower than the ceilometer CBH.

c. Boundary layer height detection, comparison with radiosounding retrievals

Figure 13 shows a comparison of the BLH estimated from radiosoundings (launched every day at 1200 UT 15 km from SIRTA) and BLH processed by the STRAT algorithm on LNA lidar data. The method used to extract BLH from soundings is a threshold method applied on the Richardson number R ib (z) [START_REF] Menut | Urban boundary layer height determination from lidar measurements over the Paris area[END_REF] calculated as

R ib z= g  z-z 0  θ  z   θ  z -θ  z 0  u z 2 v  z  2 ,

23

where θ is the potential temperature, g is the acceleration due to gravity, z is the height, z o is height of the surface and u and v are the zonal and meridian wind components. The BLH is estimated with a threshold value of 0.21 [START_REF] Vogelezang | Evaluation and model impacts of alternative boundary-layer height formulations[END_REF].

The lidar derived BLH is the median BLH extracted between t 0 and t 0 + 5' [START_REF] Menut | Urban boundary layer height determination from lidar measurements over the Paris area[END_REF] where t 0 is the radiosonde launch time.

We study 200 temporally and spatially collocated radiosonde and lidar profiles. 125 situations correspond to clear-sky events without clouds below 5000 m and without aerosol layers above the boundary layer in a 20-min window around the RS launch.

The mean difference between lidar and RS derived BLH estimates is 99 m with a standard deviation σ of 452 m, hence the standard error in the mean is 62 m for a 95% confidence interval. We find that 83% of the population is within +/-500 m (close to 1 σ).

We assume that points beyond +/-500 m are outliers and restrict the comparison to situations when the difference is in the interval [-500 m, +500 m], the mean difference becomes 21 m, the standard deviation 200m, and the standard error is 30 m. When we further restrict the comparison to clear-sky situations, we find very similar statistics (see Table 4). The population of BLH differences contains two subgroups, one representing 83% of the situations where BLH retrievals agree within 20m +/-30m (95% confidence) and the other (17% of the population) representing cases with very large discrepancies (between 500 and 1500 m). The inconsistency between the two retrieval methods can be due to lack of mixing or entrainment of aerosols, or poor collocation of radiosonde and lidar profiles. Table 4 shows that the presence of cloud does not introduce additional discrepancies.

Conclusion

The STRAT algorithm has been developed to analyze large data sets of lidar backscatter profile and retrieve the vertical structure of particle layers in the atmosphere. The algorithm is based on four successive detections carried out on individual profiles. The signal noise level is a key parameter in the algorithm as the detection thresholds at each step of the process are determined with respect to it.

Hence the algorithm automatically adjusts to varying levels of signal noise. Molecular or (near) particle free layers are determined with a conservative approach to minimize false detections so that those layers can effectively be used for automated normalization processes. Identification of particle layers is done by using continuous wavelet transforms to identify discontinuities in the lidar profile and choosing those that effectively correspond to cloud or aerosol layer boundaries. We find good consistency between cloud-base heights retrieved by STRAT and those provided by a commercial ceilometer analysis. The height uncertainty inherent to the method is evaluated to be less than 3 times the vertical resolution (e.g. less than 45 m for the LNA lidar). Similarly the transition from the boundary layer to the free troposphere is analyzed with wavelet transforms. When compared to boundary layer heights retrieved from radiosondes, we find no significant bias in the STRAT retrievals, but the comparison reveals a large scatter due to the inconsistency between the aerosolbased and the thermodynamic-based BLH definition. Even though a few test cases have been carried out with the STRAT algorithm on 355-nm, 532-nm and 1064-nm lidar systems, with both analogue and photon counting detection systems, the true portability of the STRAT algorithm to diverse large lidar data sets is still under study.

Figure captions:

FIGURE

Figure 2

 2 Figure 2 PDF of V(r,t)/ σ 0 (t) values estimated on 1 profile per day for the entire

Figure 3

 3 Figure 3 Second derivative of a Gaussian wavelet called "Mexican hat" wavelet

Figure 5

 5 Figure 5 PDF of R(r,t)/ σ(t) values estimated on 10 profiles per day for all the

Figure 6

 6 Figure 6 (a) lidar LNA data 532 nm WFOV telescope on May 26 2003, (b) profile by

Figure 7

 7 Figure 7 PDF of dPr 2 layer values estimated on the entire database for all the range

Figure 8

 8 Figure 8 First derivative of a Gaussian wavelet.

Figure 9

 9 Figure 9 (a), (c), (e) not normalized range corrected backscattered signal, (b), (d), (f)

Figure 10

 10 Figure 10 (a) Simulated backscattered profile with an added noise, where S 1 and S 2

Figure 12

 12 Figure 12 CBH lidar estimated on the range corrected backscattered profile (dark solid

  CWT a , b r ,t =∑

r P r , t  ψ a , b r 12 where

  Pr layer 2 T 4 then layer is cloud layer

	that aerosol concentrations are not sufficient to be detected above 7.5km whereas
	cloud layers extend from 0 to 15km. So, distributions of values dPr 2 layer are due to
	different contributions: under 7.5km a combination of cloud and aerosol contributions,
	above 7.5 km only cloud contributions. Thus, the distribution of dPr 2	layer for aerosols
	layers is located between 1 and 4 where the two dashed lines are distinct whereas
	the distribution of dPr 2	layer for clouds layers is wider. To separate cloud from aerosol
	layers based on these distributions, we select a threshold value T 4 =4. So equation 18
	becomes
				{ d Pr layer 2 4 then layer is cloud layer d Pr layer 2 4 then layer is aerosol layer
				d Pr layer 2	T 4 then layer is aerosol layer	.
				18
	Figure 6c shows the dPr 2	layer values for cloud and aerosol layers observed on 26 May
	2003. Some particle layers appear with significantly stronger dPr 2 layer than others.
	2)	Threshold determination
	The PDFs of dPr 2	layer values are shown in Figure 7 for three different vertical range
	intervals. The distribution based on the complete vertical range (0 to 15 km) is
	represented by a solid line, the distribution of dPr 2	layer for layers below 7.5 km is shown
	in dashed line with square markers and the distribution of dPr 2 layer for layers above
	7.5km is drawn with a dashed line and diamond markers. Those intervals of altitude
	are used because except for exceptional events like volcanic eruptions, we assume

, 19 for particle layers below 7.5km. Above 7.5km we assume that 100% of particle layers correspond to clouds.

Table captions :

 captions 1 STRAT ALGORITHM DIAGRAM...................................................................38 315M........................................................................................................................39 FIGURE 3 SECOND DERIVATIVE OF A GAUSSIAN WAVELET CALLED "MEXICAN HAT" WAVELET.............................................................................................40 ...............................................................................................................44 FIGURE 8 FIRST DERIVATIVE OF A GAUSSIAN WAVELET.....................................45 SIMULATED BACKSCATTERED PROFILE WITH AN ADDED NOISE, WHERE S1 AND S2 ARE TWO SLOPES THAT ALLOW US TO DESCRIBE THIS PROFILE. (B) DIFFERENCE BETWEEN DETECTED CBH BY STRAT AND THE REAL CBH (C) DIFFERENCE BETWEEN DETECTED CTH BY STRAT AND THE REAL CTH....................................................................................................................48 CBHCEILOMETER......................................................................50 FIGURE 12 CBHLIDAR ESTIMATED ON THE RANGE CORRECTED BACKSCATTERED PROFILE (DARK SOLID LINE) AND THE CORRESPONDING CBHCEILOMETER ..............................................................................................................51 FIGURE 13 (A) SCATTER PLOT OF BLH RETRIEVED BY STRAT ALGORITHM ON 532 NM NFOV TELESCOPE DATA AND BLH ESTIMATED BY RADIOSONDES, FOR ALL CASES WITH CROSS MARKERS AND FOR CLEAR CASES WITH RING MARKERS (B) PDF OF BLHLIDAR-BLHRADIOSONDES FOR ALL CASES...........52

	FIGURE 2 PDF OF V(R,T)/ σ0(T) VALUES ESTIMATED ON 1 PROFILE PER DAY
	FOR THE ENTIRE DATABASE AT TWO WAVELENGTH : 532NM (GREY SOLID FIGURE 11 (A) VERTICAL DISTRIBUTION OF CLOUDS SEEN BY THE LIDAR LINE) AND 1064NM (GREY DOTTED LINE) AND PDF OF SIMULATED V(R,T)/ σ0 LNA (DARK LINE) AND THE CEILOMETER LD40 (GREY DASHED LINE) (B)
	(T) VALUES IN MOLECULAR LAYERS (BLACK SOLID LINE) WITH 2N+1 = 21 SCATTER PLOT OF CBH DETECTED BY STRAT ALGORITHM ON 532 NM NFOV
	TELESCOPE DATA AND CBH DETERMINED WITH CEILOMETER DATA, (C) GATES = FIGURE 4 (A) SIMULATED AND NOT NORMALIZED BACKSCATTERING PDF OF CBHLIDAR -
	POWER RECEIVED FOR AN IDEAL CLOUD OR AEROSOLS CASE IN
	FUNCTION OF ALTITUDE, (B) CORRESPONDING CWT COEFFICIENTS
	CALCULATED FOR DIFFERENT DILATION A (FINEST HIGH UP) AND
	DIFFERENT LOCATION OF WAVELET'S CENTER B WITH HIGHEST
	COEFFICIENTS IN WHITE AND LOWEST IN BLACK, (C) SKELETON (MAXIMA
	.42
	FIGURE 6 (A) LIDAR LNA DATA 532 NM WFOV TELESCOPE ON MAY 26 2003,
	(B) PROFILE BY PROFILE DPR22 RATIO (17), (C) LAYER BY LAYER
	DPR2LAYER RATIO, (D) FLAG OBTAINED WITH PARTICLE LAYERS
	DISTINCTION, CLOUD LAYERS ARE IN RED, AEROSOL LAYERS IN ORANGE
	(18). CEILOMETER CTH DETECTIONS ARE REPRESENTED WITH BLACK
	POINTS...43
	FIGURE 7 PDF OF DPR2LAYER VALUES ESTIMATED ON THE ENTIRE
	DATABASE FOR ALL THE RANGE (SOLID LINE), RANGE BELOW 7.5 KM
	(DASHED LINE WITH SQUARES) AND RANGE ABOVE 7.5 KM (DASHED LINE
	WITH DIAMONDS)FIGURE 9 (A), (C), (E) NOT NORMALIZED RANGE CORRECTED
	BACKSCATTERED SIGNAL, (B), (D), (F) CORRESPONDING CWT'
	COEFFICIENTS CALCULATED FOR DIFFERENT DILATION A (FINEST HIGH
	UP) AND DIFFERENT LOCATION OF WAVELET'S CENTER B WITH HIGHEST
	COEFFICIENTS IN WHITE AND LOWEST IN BLACK. (A), (B) A CLEAR CASE,

LINES) OF THE CWT MODULUS. (D,E,F) SAME AS (A,B,C) FOR A REAL BACKSCATTERING..............................................................................................................41 FIGURE 5 PDF OF R(R,T)/ σ(T) VALUES ESTIMATED ON 10 PROFILES PER DAY FOR ALL THE DATABASE AND PDF OF SIMULATED R(R,T)/ σ(T) VALUES IN MOLECULAR LAYERS (SOLID LINE)............................................................................................................................................................................................................(C), (D) A CLOUDY CASE WITH A CLOUD NEAR THE BL, (E), (F) A CLOUDY CASE WITH A CLOUD AT THE BL..................................................................................47 FIGURE 10 (A)

TABLE 1 LIDAR

 1 LNA DESCRIPTION..............................................................................53 TABLE 2 PARAMETERS OF STRAT ALGORITHM TO PROCESS LIDAR LNA DATA........................................................................................................................................54 TABLE 3 COMPARISON BETWEEN STRAT DETECTION AND CEILOMETER DETECTION...........................................................................................................................55 TABLE 4 COMPARISON BETWEEN STRAT BLH RETRIEVALS AND BLH ESTIMATED BY RADIOSONDES......................................................................................56 Figure 1 STRAT algorithm diagram.

Table 4

 4 Comparison between STRAT BLH retrievals and BLH estimated by radiosondes

				ALL	CLEAR-SKY
	Situations	ALL	CLEAR-SKY
				(+/-500m)	(+/-500m)
	Number of cases	211	125	173	105
	Mean (BLH lidar -BLH RS ) (m)	99	89	24	51
	Standard Dev. (BLH lidar -BLH RS ) (m)	452	445	202	199
	Standard Error

of the Mean (BLH lidar -BLH RS ) (m) (95% conf. Int.)

	62	80	30	38

Acknowledgments

This work has been carried out with the financial support of the French National Space Agency CNES. The authors would like to acknowledge Stéphane Mallat for stimulating discussions on our application of the CWT method. We extend our acknowledgments to the reviewers for their comprehensive reviews and useful inputs.

in Lidar Remote Sensing for Industry and Environmental Monitoring, U. N. Singh, T. Itabe, N. Sugimoto, (eds.). Proc. SPIE, 4153, 151-158. Winker D., J. Pelon and P. McCormick, 2003: The CALIPSO mission : Spaceborne Lidar for observations of aerosols and clouds. Proc. SPIE, 4893,[1][2][3][4][5][6][7][8][9][10][11]